Evolutionary and dynamic constraints converge on a conserved coupling sector in thermally adapted GH1 β-glucosidases

How protein sequences evolve to maintain functional conformational dynamics under thermal stress remains an important question in molecular evolution. Thermophilic enzymes are often considered to achieve stability primarily through enhanced structural rigidity, yet thermal adaptation may also involve the preservation or reorganization of correlated motions across the protein scaffold. Here, we compared three glycoside hydrolase family 1 β-glucosidases-BsBgl1A, 9UPT, and LQBG8-that share a conserved (β/α)8 TIM-barrel fold but display distinct optimal temperatures of approximately 45, 65, and 70 °C, respectively. By integrating comparative structural analysis, three independent 1 μs molecular dynamics simulations at 300, 325, and 350 K, neural relational inference of residue-level dynamic coupling, and family-level coevolutionary analysis, we identified distinct dynamic strategies associated with their thermal phenotypes. All three enzymes retained their global compactness and principal secondary-structure framework across the simulated temperatures, indicating that differences in RMSD and residue flexibility did not reflect large-scale unfolding. BsBgl1A showed the strongest temperature sensitivity, with increased structural deviation and residue mobility at 350 K. Its inferred interaction network exhibited a non-monotonic response, characterized by pronounced overcoupling at 325 K followed by reduced connectivity at 350 K. In contrast, 9UPT maintained relatively restrained residue dynamics and a highly reproducible coupling topology across temperatures. LQBG8 combined substantial local conformational flexibility with persistent long-range coupling and retained the highest directed network density at 350 K. Across the three enzymes, the βα5-βα7 region repeatedly emerged as a shared structural and dynamic center. Coevolutionary couplings were also enriched in TIM-barrel-associated regions that spatially converged on this sector. These findings indicate that thermal adaptation in GH1 β-glucosidases is not adequately explained by uniform structural rigidification. Instead, enzymes with the same conserved fold may adopt distinct organizations of correlated motion, while greater thermal adaptation is associated with the persistence or reorganization of a βα5-βα7-centered dynamic coupling pattern under elevated temperature. The correspondence between dynamic coupling and coevolutionary enrichment is consistent with long-term evolutionary constraint on this structural sector, although it does not establish direct selection on the network as an independent functional unit. Our results therefore provide a testable framework in which thermostability reflects a balance among local stabilization, regional flexibility, and the resilience of higher-order dynamic organization.