Convergence, stability, and thermal adaptation in the rubisco enzyme in plants

Convergence, stability, and thermal adaptation in the rubisco enzyme in plants | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Convergence, stability, and thermal adaptation in the rubisco enzyme in plants View ORCID Profile Arthur Leung , View ORCID Profile Belinda S. W. Chang , View ORCID Profile Rowan F. Sage doi: https://doi.org/10.1101/2025.10.08.681247 Arthur Leung 1 Department of Ecology and Evolutionary Biology, University of Toronto , 25 Willcocks Street, Toronto, ON M5S 3B2, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Arthur Leung For correspondence: art.leung{at}mail.utoronto.ca Belinda S. W. Chang 1 Department of Ecology and Evolutionary Biology, University of Toronto , 25 Willcocks Street, Toronto, ON M5S 3B2, Canada 2 Department of Cell and Systems Biology, University of Toronto , 25 Harbord Street, Toronto, ON M5S 3G5, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Belinda S. W. Chang Rowan F. Sage 1 Department of Ecology and Evolutionary Biology, University of Toronto , 25 Willcocks Street, Toronto, ON M5S 3B2, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Rowan F. Sage Abstract Full Text Info/History Metrics Preview PDF Abstract Enzymes are thought to be tuned to perform similarly in different thermal regimes. Whether the photosynthetic enzyme ribulose-1,5-bisphosphate carboxylase/oxygenase (rubisco) follows similar rules, especially when considering evolutionary history, is uncertain. The molecular, structural, and ecological factors of the rubisco large subunit (RbcL) were examined in four plant clades: wood ferns, pines, sea lavenders, and viburnums. Using rbcL gene sequences, codon evolutionary models were used to test for positive and divergent selection and convergent evolution. Protein structure modeling was performed to predict side chain changes and protein stability. Phylogenetic comparative methods were used to examine the relationship between protein stability to growing season temperature. All four clades showed significant evidence of positive selection, with multiple convergent substitutions predicted to alter side chain polarity and interactions with the solvent. In viburnums, biome transition rates were dependent on amino acid substitution, with positive selection was concentrated in cold temperate and cloud forest clades. Rubiscos with higher stability occurred in species from warmer environments. However, this correlation was weaker after correcting for phylogeny. These analyses support a hypothesis that RbcL evolution is influenced by both environmental tuning and evolutionary history. Introduction The diverse range of environmental conditions on the planet are thought to contribute to the evolution of enzyme performance, leading to adaptations that tune function such as catalytic rate and substrate affinity. These properties are represented by the turnover rate per active site ( k cat ) and the inverse of the Michaelis constant ( K m ). Enzyme stability and flexibility are key determinants of turnover rate and affinity. Flexible enzymes bind substrate and release product faster but are more prone to denaturing or adopting low-affinity conformations ( Fields, 2001 ; Fields and Houseman, 2004 ). Variations in temperature disrupt this stability-flexibility balance. At supraoptimal temperatures, non-covalent interactions (e.g., hydrogen bonds) that stabilize the enzyme are weakened, initially leading to elevated reaction rates but decreased affinity; in contrast, at suboptimal temperatures, non-covalent interactions are strengthened, reducing time spent in low-affinity conformations, but inhibiting catalytic motions that contribute to greater k cat ( Fields, 2001 ). The temperature response of reaction rate is typically exponential and doubles with every 10 °C rise in temperature ( Scholander et al., 1953 ). However, in a comparative context, proteins in different species do not tend to follow this relationship. For example, a fish that typically lives at 0 °C does not have a metabolic rate that is 20 times slower than one that typically lives at 40 °C ( Fields et al., 2015 ). Evolutionary selection drives distinct adaptations in organisms from different thermal habitats such that stability and flexibility is optimized for their environment ( Hochachka and Somero, 2002 ). As shown for lactate dehydrogenase (LDH) and malate hydrogenase (MDH), cold-adapted enzymes are more flexible and less stable than warm-adapted orthologs, enabling similar performance ( Fields and Somero, 1998 ; Dong and Somero, 2009 ; Liao et al., 2017 ; Dong et al., 2018 ). Thus, protein evolution allows organisms to escape the “tyranny” of the exponential temperature response ( Bullock, 1955 ). While LDH and MDH are classic examples of thermal adaptation in enzymes, ribulose-1,5-bisphosphate (RuBP) carboxylase/oxygenase (rubisco), which catalyzes CO 2 fixation in photosynthesis, presents additional challenges due to its dual-substrate specificity. The carboxylation of one RuBP molecule produces two phosphoglycerate molecules, which is used in the photosynthetic carbon reduction cycle to produce glucose. Oxygen competitively inhibits RuBP carboxylation, and each RuBP oxygenation event produces one phosphoglycerate and one phosphoglycolate molecule ( Bowes and Ogren, 1972 ). Phosphoglycolate is toxic and must be recycled into phosphoglycerate, which consumes energy and releases previously fixed CO 2 in the process. Importantly, the k cat and K m for CO 2 and O 2 ( k cat c , K c and k cat °, K o , respectively) do not scale uniformly to temperature ( Jordan and Ogren, 1984 ). However, variation in k cat c and CO 2 specificity factor ( S c/o ) across species indicates evolutionary flexibility ( Jordan and Ogren, 1981 ; Galmés et al., 2016 ). In an extreme case, CO 2 specificity varies by an order of magnitude between Rhodospirillum rubrum and land plants ( Jordan and Ogren, 1981 ). Rubisco small subunit (RbcS) isoforms can further modulate kinetics in plants, adding an another layer of complexity ( Lin et al., 2020 ; Cavanagh et al., 2022 ). Whether these patterns are consistent with broader trends in enzyme evolution, such as those in LDH and MDH ( Hochachka and Somero, 2002 ), remains unresolved. The complexity of rubisco kinetics becomes more apparent when considering how variation in CO 2 availability, shaped by environmental pressures and photosynthetic pathways, influences its evolution. Plants with CO₂-concentrating mechanisms (CCMs), such as C₄ and CAM species, operate at higher plastidial CO₂ concentrations than C₃ plants, which has driven rubisco adaptation toward higher k cat c but lower K c ( Yeoh et al., 1980 , 1981 ; Iñiguez et al., 2020 ). High CO 2 affinity is less important than enzyme velocity when the CO 2 availability at the active sites is high, and thus these kinetic differences have convergently evolved multiple times ( Christin et al., 2008 ). In plants that lack CCMs, rubisco adaptation is constrained by the relative rates of carboxylation to oxygenation and the trade-off between k cat c and S c/o ( Lorimer and Andrews, 1973 ; Tcherkez et al., 2006 ; Savir et al., 2010 ). It is thought that better differentiation of CO 2 and O 2 transition states strengthens the binding with transition intermediates, slowing subsequent steps in the reaction ( Tcherkez et al., 2006 ). For instance, C 3 plants in warm, saline, or arid environments, where stomatal closure limits diffusion of CO 2 substrate to the active site, tend to have rubiscos with lower k cat c but higher S c/o or lower K c ( Sage, 2002 ; Galmés et al., 2005 ). However, this k cat - K c trade-off is weaker in recent, larger datasets ( Flamholz et al., 2019 ), which would raise questions regarding if and how evolutionary history can shape RbcL evolution. A comparative phylogenetic approach offers a powerful means of understanding the evolution of rubisco structure and function in relation to temperature. It can be used to identify sites in the protein that were the targets of positive selection, and to locate amino acid shifts onto specific lineages of the phylogeny in distant relatives ( Kapralov and Filatov, 2007 ; Galmés et al., 2014a ; Hermida-Carrera et al., 2016 ; Orr et al., 2016 ; Sharwood et al., 2016 ) and within genera ( Kubien et al., 2008 ; Galmés et al., 2014a ; Capó-Bauçà et al., 2023 ). Although amino acid changes in rubisco can be linked to shifts in enzymatic function, the connection between amino acid variation and its broader ecological context is less frequently explored (but see Hermida-Carrera et al., 2017 ). Phylogenetic comparative analyses are thus useful, because evolutionary relatedness between proteins can confound the interpretation of rubisco evolution ( Tcherkez et al., 2006 ; Savir et al., 2010 ; Bouvier et al., 2021 , 2024 ; Tcherkez and Farquhar, 2021 ; Bouvier and Kelly, 2023 ). For instance, the relationship of k c with K may differ when compared across broad versus narrow taxonomic scales ( Galmés et al., 2014b ). Thus, evolutionary history must be explicitly accounted for when correlating amino acid sequence with different thermal niches ( Felsenstein, 1985 ; Pagel, 1994 ). This framework allows us to account for phylogenetic non-independence and to better resolve whether observed trait correlations are adaptive or simply reflect the degree of relatedness. For instance, distantly relatedly species may differ in the genetic background other than RbcL, which could confound differences in k cat c and K c ( Hermida-Carrera et al., 2016 ) On the other hand, examining closely related species without correcting for relatedness can lead to misleading conclusions about adaptation, since similarities may be due to shared ancestry rather than adaptation to environmental pressures ( Bouvier et al., 2021 ). Based on this rationale, we selected phylogenetically disparate plant genera each spanning multiple habitats: wood ferns ( Dryopteris ), sea lavenders ( Limonium ), pines ( Pinus ), and viburnums ( Viburnum ). We asked the following questions: (1) Is adaptive evolution responsible for RbcL diversity? (2) Is RbcL diversity associated with the thermal environment of vascular plant species? (3) How do these patterns hold up when correcting for phylogenetic history? By using phylogenetics to examine structural and evolutionary patterns in tandem, we aim to evaluate how molecular changes in rubisco explain its adaptation. We hypothesized that convergent amino acid substitutions would be consistent with an increased the number of non-covalent interactions in the holoenzyme that enhance stability in warm or dry climates. By integrating analyses of selection, convergence, stability, and correlated evolution, our study provides insight into the evolutionary and ecological basis of rubisco diversity. Materials and methods Taxon sampling and analytical approach We used a model clade approach ( Donoghue and Edwards, 2019 ; Mabry et al., 2024 ). Patterns of positive selection on a large phylogeny consisting of plants and algae ( Kapralov and Filatov, 2007 ) may differ from those in closely-related species ( Galmés et al., 2014a ). Sparse taxon sampling of distantly related species can lead to the misinterpretation of other evolutionary change, such as that in rubisco small subunit (RbcS) and rubisco activase (Rca), with changes associated with the focal trait ( Heyduk et al., 2019 ). We chose four phylogenetically distant genera that each were speciose and had well-resolved phylogenies, such that each can be examined as separate examples of evolution, and such that effects of saturation are minimized. The four genera, Dryopteris , Limonium , Pinus , and Viburnum, were chosen. Tests of saturation showed minimal saturation in Dryopteris , Pinus , and Viburnum ( Table S 1 ). Limonium showed evidence of saturation, which could lead to overestimation of positive selection; thus, the significance of the results may be weakened. Gene sequences and alignments rbcL gene sequences from 621 vascular plant accessions were downloaded from GenBank: 185 Dryopteris , 184 Limonium , 104 Pinus , and 148 Viburnum accessions (see Data Availability for accession numbers). Separately for each genus, sequences were aligned using MUSCLE v5 ( Edgar, 2022 ) and then manually inspected in MEGA v11 ( Tamura et al., 2021 ). To facilitate homology modelling, the sequences were aligned with a spinach rbcL sequence (GenBank: NC 002202), which was later removed from the alignment. Species-level phylogenetic trees for Pinus and Viburnum were downloaded from supplementary data of respective papers ( Jin et al., 2021 ; Landis et al., 2021 ). The Dryopteris and Limonium phylogenies were provided by the lead authors ( Sessa et al., 2017 ; Koutroumpa et al., 2021 ); both lead authors granted written permission to make these files available in the data repository (see Data Availability). Codon evolutionary analyses Random sites models were used to analyze codon evolution at the site level using CodeML in ‘PAML’ v4.10 ( Yang, 2007 ). Codon model (M) 0, M1a, M2a, M3, M7, and M8 were used with codon frequencies empirically determined from the rbcL sequences. For each genus, the best fitting codon model with the lowest Akaike Information Criterion (AIC) value was used to reconstruct of the ancestral state at the root of the phylogeny. Likelihood-ratio test results of M1 vs. M0, M2a vs. M1a, M3 vs. M0, and M8 vs. M7 are reported in Table 1 . Bayes empirical Bayes (BEB) analysis determined the probability that each codon site was under positive selection; codon sites with a BEB posterior probability ( P ) > 0.95 were considered under positive selection. View this table: View inline View popup Download powerpoint Table 1. Likelihood ratio tests of PAML codon models in each genus ( Yang, 2007 ). For models M0–M3, ω 0 – ω 2 values are shown representing d N / d S values for each site class, along with proportion of sites for each site class. For models M7 and M8, α and β describe the shape of the beta distribution, and ω P describes the d N / d S value for the positively selected site class. n.p., number of parameters; AIC, Akaike Information Criterion; ln L , log-likelihood; κ, transition/transversion ratio; χ 2 , chi-squared; d.f., degrees of freedom; p , p -value of likelihood ratio test. Protein convergence analyses Using ‘CSUBST’ v1.4.0 ( Fukushima and Pollock, 2023 ), ancestral gene sequences were reconstructed to locate branch pairs at which convergent non-synonymous and synonymous substitutions occurred (any amino acid to a specific amino acid) and calculate the rate at which they occur for each codon site (hereafter, d Nc and d Sc , respectively). ‘CSUBST’ internally uses ‘IQTREÈ to determine the best codon substitution model for ancestral state reconstruction ( Minh et al., 2020 ; Fukushima and Pollock, 2023 ). In contrast to commonly used methods to examine rbcL convergence, this method corrects for phylogenetic errors due to introgression ( Fukushima and Pollock, 2023 ). Introgression has been shown to affect rbcL evolution ( Yao et al., 2019 ) and would similarly affect both the d Nc and d Sc , so not accounting for the d Sc would overestimate the rate of convergent codon evolution. In contrast, convergent amino acid substitutions that arise from positive selection should increase the d Nc relative to the d Sc . Thus, convergent amino acid substitutions were defined as a substitution at the same site in a branch pair any amino acid (empirical Bayesian P > 0.95) to the same amino acid, with a d Nc / d Sc > 1 ( Fukushima and Pollock, 2023 ). Amino acid- and codon-level associations with biome transitions in Viburnum For Viburnum , discrete biome scorings were obtained by binarizing published data ( Landis et al., 2021 ). Discrete representations of thermal niche based on expert observations were unavailable for the other three genera at the time of this work. Viburnum species occur in warm temperate forests, cold temperate forests, cloud forests, and tropical forests. The genus first diversified in warm temperate forests and later diversified into the other biomes ( Landis et al., 2021 ). To further examine the associations of codon evolution with transitions between biomes, clade models C and D (CmC, CmD) in CodeML were used. These models test for divergence in d N / d S between clades ‘PAML’ v4.10 ( Yang, 2007 ). Clade models allow for one class of codon sites to vary in d N / d S ( ω ) between preassigned clades (e.g., clades occupying different biomes). Compared to branch-site models, which test for rapid positive selection on specific branches, clade models can test for long-term selection across clades ( Chang et al., 2012 ). CmC allows three classes of sites with 0 < ω 0 0, with ω 2 varying between clades. CmD relaxes the assumptions of the first two classes, allowing ω 0 , ω 1 > 0 ( Bielawski and Yang, 2004 ). The null models of CmC and CmD are M2a_rel ( Weadick and Chang, 2012 ) and M3, both of which constrain ω 2 to be identical across clades. Near-identical parameter estimations resulted from different initial values of transversion/transition ratio (0.3, 10) and/or ω (1.3, 10) (see Data Availability). The dependence of biome evolution on amino acid evolution across the Viburnum phylogeny ( Landis et al., 2021 ) was tested by fitting models of evolution in which they were dependent or independent. Binary character evolution models ( Pagel, 1994 ) were fitted using the ‘fitPagel’ function in ‘phytools’ v1.9-9 in R ( Revell, 2012 ). The Pagel (1994) test provides a conservative test of correlated evolution by counting independent transitions rather than inflating sample size with phylogenetically non-independent species observations, making it appropriate for detecting coevolutionary patterns. The models can accommodate cases where biome or amino acid transitions occur before or after the other. Tests evaluated the dependence of biome transition rates on amino acid transition rates and vice versa, alternately designating amino acid or biome as the independent variable; models were selected based on the lowest AIC values. Likelihood ratio tests compared the model in which biome and amino acid evolve independently versus a model in which the two characters evolve dependently. Homology modelling The effects of amino acid substitution on non-covalent interactions were visualized using ‘PyMol’ v3.0.4 ( Schrodinger, 2015 ). Point mutations that add or remove oxygen or sulfur atoms were introduced to the large subunit of a 1.6 Å spinach rubisco structure containing residues 9– 475 in the large subunit (PDB: 8RUC). The lowest energy rotamer was used. Substitutions were assessed based on whether they altered the number of non-covalent interactions within 3.7Å of the amino acid residue. Protein stability modelling The free energy of protein folding (Δ G fold ) of the spinach rubisco structure was predicted at 298 K using ‘FoldX’ v5.1 ( Delgado et al., 2019 ). ‘FoldX’ is widely used for protein stability predictions and its predictions reliably reproduce experimental data ( Khan and Vihinen, 2010 ; Buß et al., 2018 ), but it is important to interpret modelled results with caution. Following Studer et al. (2014) , heteroatoms that are not well parameterized in FoldX (carboxyarabinitol bisphosphate and Mg 2+ ) were removed from the structure, with interacting protein residues constrained in their conformation. The residues interacting with these heteroatoms according to Andersson (1996) were constrained in space as done previously ( Studer et al., 2014 ). Using ‘RepairPDB’ in ‘FoldX’, high energy residues were adjusted to its lowest energy rotamer ( Studer et al., 2014 ). Modeling mutations on a single reference structure is standard practice for comparative protein evolution studies ( Studer et al., 2014 ) and is appropriate here because most positively selected sites occur at surface-exposed positions ( Figure 1 ) where local structural context is less constrained than in the catalytic core. Download figure Open in new tab Figure 1. Sites under positive selection and/or convergent evolution in Dryopteris , Limonium , Pinus , and/or Viburnum overlayed on a published crystal structure of rubisco large subunit protein bound to transition-state analog carboxyarabinitol bisphosphate (CABP) from spinach (PDB: 8RUC) ( Andersson, 1996 ). For each substitution at sites under positively selection and/or convergent evolution, all large subunits in the rubisco holoenzyme were mutated, with one substitution at a time, using ten replicate runs of ‘BuildModel’ in ‘FoldX’. The free energies of protein folding in the wildtype (Δ G fold,wt ) and the mutant protein (Δ G fold,mut ) were predicted. The reported ΔΔ G fold = Δ G fold,wt − Δ G fold,mut values were divided by eight to standardize the values to each large subunit ( Studer et al., 2014 ). For each amino acid substitution, the wildtype residue was set to the ancestral residue in the relevant clade, while the mutant residue was set to the derived residue. For sites where the derived residue matched the residue in spinach, the ΔΔ G was assumed to be path-independent (e.g., ΔΔ G N76→S = −ΔΔ G S76→N ) ( Studer et al., 2014 ). For sites where neither ancestral nor derived residues match the residue in spinach, the ΔΔ G values were assumed to be additive and path-independent (e.g., ΔΔ G A328→G = ΔΔ G S328→G − ΔΔ G S328→A ) ( Studer et al., 2014 ). A separate analysis predicted the ΔΔ G for each species relative to the spinach homology model. For all RbcL gap-free sequences in this study ( n = 399 species; n = 174 unique protein sequences; mean sequence coverage = 92% of the full coding sequence), the spinach crystal structure was mutated to each unique sequence in ‘FoldX’. For the remaining 8% of missing data at the N- and C-termini, the spinach wild-type sequence was used. The analysis was parallelized locally using GNU Parallel ( Tange, 2023 ). Growing season temperature estimation Selection on rubisco performance was assumed to be strong during warm and wet months, when stomata are open and leaves are photosynthetically active. We used interpolated climate data from WorldClim as a standardized, reproducible metric of thermal environment, following established protocols where species-specific microclimate data are unavailable ( Fick and Hijmans, 2017 ; Wright et al., 2017 ; He et al., 2020 ). We downloaded geographical coordinates corresponding to species occurrences for each genus ( GBIF, 2023 , 2024a , b, c). Incomplete or fuzzy species name searches on GBIF were manually inspected. ‘CoordinateCleaner’ was used to exclude duplicate coordinates, (0,0) coordinates, points with identical latitude and longitude values, oceanic points, points overlapping with biodiversity institutions, points near country capitals, and points near country centroids. For species with three or more occurrences, monthly temperature and aridity index values were downloaded from WorldClim and Global Aridity Index and Potential Evapotranspiration Database ( Fick and Hijmans, 2017 ; Zomer et al., 2022 ). The growing season was defined the months with monthly mean temperature ≥ 5°C and aridity index > 0.05 ( Wright et al., 2017 ; He et al., 2020 ). The mean temperature of the warmest month in the growing season was used. Phylogenetic generalized least squares regressions To examine the evolution of RbcL stability across genera, a backbone tree was created with crown ages of Dryopteris , Limonium , Pinus , and Viburnum with divergence times from Timetree.org ; previously published species trees ( Sessa et al., 2017 ; Jin et al., 2021 ; Koutroumpa et al., 2021 ; Landis et al., 2021 ) were then grafted onto the backbone tree (see Data Availability). Statistical tests were performed using phylogenetic generalized least squares regressions under a Brownian motion model implemented in ‘apè with post-hoc tests performed using estimated marginal means implemented in ‘emmeans’ ( Searle et al., 1980 ; Paradis and Schliep, 2019 ). Results Pervasive positive selection in rbcL To test for positive selection at the codon level in the rbcL gene, we used PAML to fit random sites models of codon evolution to the alignment and phylogeny corresponding to each genus ( Yang, 2007 ). Null models in which 0 ≤ d N / d S ≤ 1 (M1a and M7) were compared to corresponding models with an additional class of codons under positive selection ( d N / d S > 1; M2a and M8, respectively). We also fitted a model that imposes no specific constraints on three categories of d N / d S (M3), as a potential model to use for ancestral state reconstruction. As summarized in Table 1 , the likelihood ratio tests of random sites models show at least one or both of the positive selection models (M2a and M8) were found to be a significantly better fit than the respective null models ( p < 0.001), providing evidence for positive selection in rbcL across all four genera. A Bayes empirical Bayes (BEB) procedure was also used to calculate the posterior probability ( P ) of each site being under positive selection. Sites with P > 0.95 were categorized as under positive selection ( Table 1 ; Figure 2 ); these sites had d N / d S values ranging from 6 to 10, indicating a high rate of amino acid substitution. The number of sites under positive selection were 5 in Dryopteris , 5 in Limonium , 14 in Pinus , and 19 in Viburnum ( Figure 2 ). As visualized on the crystal structure, many of these sites were located on the surface of the large subunit and few were around the active site ( Figure 1 ). Download figure Open in new tab Figure 2. The non-synonymous substitution rate ( d N ): synonymous substitution rate ( d S ) ratio at each codon site of rubisco large subunit in (A) Dryopteris , (A) Limonium , (C) Pinus , and (D) Viburnum . Symbols and labels indicate sites under positive selection (Bayes empirical Bayes P > 0.95). Asterisks (*) indicate sites under convergent evolution (empirical Bayes P > 0.95). High levels of protein convergence within genera While positive selection suggests adaptive evolution, convergent evolution would indicate that amino acid substitutions occurred due to similar selective pressures rather than shared ancestry. This was tested using ancestral codon reconstruction at each branch and estimating a convergent d N / d S (i.e., d Nc / d Sc ; see Materials and Methods) across all branch pairs. Using the d Sc corrects for introgression and reduces the likelihood of overestimating convergence. Using this method, most of the sites under positive selection were also identified to be under convergent evolution (empirical Bayes P > 0.95). Each convergent substitution typically occurred across two to six branches ( Table S 3 ). However, parallel evolution in rubisco at the codon level was sometimes extensive; for instance, A262→V occurred across 19 different branches in the Dryopteris phylogeny ( Table S 3 ). Additionally, at one site, we observed convergence towards two different derived amino acids, e.g., A/C449→S and C/S449→A in Pinus ( Table S 3 ). The high probability of protein convergence at multiple sites provides further support for the adaptive evolution of the rbcL gene. Codon and amino acid substitutions were correlated with biome transitions in Viburnum Clade models C and D (CmC, CmD) were used to test for divergences in d N / d S between clades in viburnums ( Bielawski and Yang, 2004 ). This genus of plants first diversified in warm temperate (lucidophyllous) forests and subsequently underwent multiple independent radiations into cold temperate and tropical forests, as well as a single radiation into cloud forests ( Landis et al., 2021 ). Based on published ancestral biome reconstructions, these derived biomes were set as foreground clades, with warm temperate forests set as the background ( Landis et al., 2021 ). Likelihood ratio tests comparing these models to their respective null models are shown in Table 2 . View this table: View inline View popup Download powerpoint Table 2. Likelihood ratio tests of PAML clade models C and D (CmC, CmD) in Viburnum compared to their null models M2a_rel and M3 ( Yang, 2007 ; Weadick and Chang, 2012 ). ω 0 – ω 2 values are shown representing d N / d S values for each site class, along with proportion of sites for each site class. In M2a_rel, ω 2 is shared across all clades, whereas in CmC and CmD, ω 2 is estimated separately for each class of clade. n.p., number of parameters; AIC, Akaike Information Criterion; ln L , log-likelihood; κ, transition/transversion ratio; χ 2 , chi-squared; d.f., degrees of freedom; p , p -value of likelihood ratio test. See Table 1 for M3 parameters in Viburnum . CmC fit significantly better than the null model M2a_rel ( p = 0.0129), indicating divergence in the strength of selection between biomes ( Table 2 ). All clades were under positive selection as indicated by a class of sites d N / d S > 1 ( Table 2 ). Clades in cold temperate and cloud forests showed stronger positive selection than those in warm temperate forests, with d N / d S values for the freely varying site class ( ω 2 ) of 20 and 161, respectively, compared to 17 in warm temperate forests ( Table 2 ). In contrast, positive selection was relaxed in tropical clades, which had a lower ω 2 of about 9 ( Table 2 ). CmD, which allows for all site classes to vary freely, did not fit the data better compared to its null model M3 ( p = 1.00), but showed a similar trend in ω 2 values ( Table 2 ). To further examine this association, binary tests for correlated evolution ( Pagel, 1994 ) were used to evaluate the dependence of transition rates of amino acid on biome and biome on amino acid at each site under positive selection or convergent evolution. Models in which amino acid and biome were interdependent had poor model fit for all sites (data not shown). Likelihood ratio tests indicated that 12 substitutions in viburnums were associated with biome transitions ( Table S 4 ). Of the twelve amino acid substitutions that exhibited significant correlations with biome, eight showed AIC values supporting dependency of biome on amino acid. It was therefore predicted that changes in protein evolution had a greater impact on biome shifts than the reverse ( Table S 4 ). Four substitutions, S76→N, V157→S, V255→I, and A328→S, were associated with faster transitions into warmer environments ( Table S 4 ). M309→I was associated with faster transitions into colder environments. Stochastic character mapping predicted a greater number and rate of independent transitions into warmer environments in those derived amino acids ( Figure 3A,C,H,K ). In contrast, the M309→I substitution was almost always associated with cold biomes ( Figure 3J ). The remaining seven substitutions had derived amino acids that exhibited symmetrical transition rates ( Figure 3B,D,E,F,G,I,L ; Table S 4 ). Together, These findings suggest that substitutions in RbcL play a role in biome transitions in viburnum. Download figure Open in new tab Figure 3. Stochastic character mapping of correlated amino acid substitutions and biome shifts in Viburnum in Viburnum based on Pagel (1994) tests. The best ‘fitMk’ model in ‘phytools’ was used in 100 stochastic mapping simulations. The reconstructed ancestral state for the amino acid from ‘PAML’ and biome from Landis et al. (2021) was fixed at the root. Black, ancestral amino acid or warm temperate, cloud, or tropical forests; yellow, derived amino acid; green, cold temperate forests. Amino acid substitutions influenced non-covalent interactions To determine whether amino acid substitutions at positively selected or convergent sites influenced non-covalent interactions, we conducted homology modeling using a spinach crystal structure ( Andersson, 1996 ) to assess potential impacts on protein rigidity. Eleven amino acid substitutions involved shifts towards increased ability to form hydrogen bonds ( Figure 4 ), whereas three substitutions resulted decreased hydrogen bonding ability ( Figure 5 ). Possible hydrogen bonding interactions were examined with other side chains within 3.7Å, a distance that encompasses weak, moderate, and strong hydrogen bonds ( Jeffrey, 1997 p. 12). Download figure Open in new tab Figure 4. Non-covalent interactions introduced by amino acid substitutions in Dryopteris , Limonium , Pinus , and Viburnum , visualized on a spinach homology structure of the rubisco large subunit (PDB: 8RUC) ( Andersson, 1996 ). The derived amino acid (bolded) is shown for sites that ancestrally do not contain O, N, or S atoms. Grey spheres, the position of water molecules; numbers next to dashed yellow lines, distance in Å between O, N, or S atoms. Download figure Open in new tab Figure 5. Non-covalent interactions removed by amino acid substitutions in Dryopteris , Limonium , Pinus , and Viburnum , visualized on a spinach homology structure of the rubisco large subunit (PDB: 8RUC) ( Andersson, 1996 ). The ancestral amino acid (bolded) is shown for sites where the derived amino acid does not contain O, N, or S atoms. Grey spheres, the position of water molecules; numbers next to dashed yellow lines, distance in Å between O, N, or S atoms. Most amino acid substitutions that increased hydrogen bonding capability involved the replacing non-polar aliphatic or aromatic residues for residues containing hydroxyl groups. About half of these replacements (S143, S157, T230, T262, Y282, and S328) were predicted to form new non-covalent interactions within the same subunit, with another large or small subunit, or with water molecules ( Figure 4C,D,F,G,I,J ). In contrast, non-covalent interactions were not predicted to be affected by other substitutions for polar residues (see Figure 4A,B,E,H,K ). Conversely, substitutions of serine, tyrosine, or cysteine residues (S95, Y226, and C449) for non-polar aliphatic or aromatic residues were predicted to remove non-covalent interactions ( Figure 5 ). Multiple water molecules were predicted to be involved in networks of non-covalent interactions ( Figure 4C,D,F,I ; Figure 5B ), indicating that amino acid substitutions may alter solvent interactions by adding or removing water bridges. Amino acid substitutions were associated with protein stability changes To quantify the extent to which each amino acid substitution affects holoenzyme stability, we modelled the ΔΔ G fold , an index of stability and ease of protein folding. The change in ΔΔ G fold from the ancestral to derived amino acid in each substitution is shown in Figure 6 . Derived residues that add oxygen, nitrogen, or sulfur atoms were predicted to increase stability contributions from hydrogen bonds by about −0.74 kcal mol -1 per large subunit ( Figure 6A ), by locally stabilizing specific chains. Overall, the most substantial effect of substitutions was on hydrophobic group interactions with the solvent, with a predicted mean destabilization of +1.7 kcal mol -1 large subunit -1 ( Figure 6B ); the addition of polar, hydrophilic atoms may therefore perturb the exposure or burial of hydrophobic groups, thereby affecting stability. Altogether, the total ΔΔ G fold , which sums all protein interactions (see Table S 1 for the remaining contributors to stability), indicates that the addition of polar atoms resulted in a predicted mean destabilization of about +1.8 kcal mol -1 large subunit -1 for each residue ( Figure 6C ). Download figure Open in new tab Figure 6. Modelled change in free energy of folding (ΔΔ G fold ) of the rubisco holoenzyme from ancestral to derived amino acid in Dryopteris , Limonium , Pinus , and Viburnum . ΔΔ G fold values were predicted for each amino acid substitution at sites under positive selection and/or convergent evolution. ΔΔ G fold values are grouped by: ancestral amino acid side chain unable to form hydrogen bonds and derived able (“gain”); both ancestral and derived either able or unable (“none”); and ancestral able and derived unable to form hydrogen bonds (“loss”). Boxplots show the distribution of ΔΔ G fold (from three replicate simulations), standardized to each of eight large subunits in the holoenzyme. Symbols represent the mean for each category of substitution. (A) and (B), contribution from hydrogen bonds and interaction of hydrophobic groups with the solvent, respectively; (C) total ΔΔ G fold . See Table S2 for other components. One sample t-tests were used to test if each type of substitution had ΔΔ G f old values different from zero; *, p < 0.05; **, p < 0.01; ***, p < 0.001. Stability was weakly correlated with growing season temperature To evaluate the additive effects of amino acid diversity on protein stability, the large subunit sequence in the spinach homology model was mutated to predict the ΔΔ G fold of each gap-free sequence in this study. The small subunit sequences were not mutated. Across wood ferns, sea lavenders, pines, and viburnums, maximum growing season temperature had a modest but statistically significant correlation with modeled ΔΔ G fold ( p = 0.0013; R 2 = 0.028; Figure 7A ). Given the low variance explained and loss of signal after phylogenetic correction, these associations should be interpreted with caution. Climatic variables interpolated from weather station data is a coarse representation of the actual growing conditions. The most stable tercile of rubiscos occurred in environments that are 2.4 °C warmer than those of the least stable tercile. ( Figure 7B ). The relationship between ΔΔ G fold and growing season temperature, when analyzed using PGLS was altered when using a binned approach ( Figure 7B ) and disappeared when using continuous data. Together, the analyses indicate a strong phylogenetic effect in rubisco thermal adaptation. Download figure Open in new tab Figure 7. (A) Correlations between mean temperature of the warmest month in the growing season and change in free energy of folding (ΔΔ G fold ) of rubisco. The solid line is a linear regression, and the grey shaded area represents the 95% confidence interval. This correlation is not significant using phylogenetic generalized least squares (PGLS). Holoenzyme ΔΔ G fold values were predicted with each gap-free rbcL sequence ( n = 399 species with 174 unique sequences) by mutating the spinach crystal structure (PDB: 8RUC). The mean values from three replicate runs in ‘FoldX’ are shown. The growing season was defined as the months with monthly mean temperature ≥ 5°C and aridity index > 0.05 ( Wright et al., 2017 ; He et al., 2020 ). See Figure S 1 for regressions on each genus. (B) ΔΔ G fold values were also analyzed with a binning approach, which is standard practice ( Studer et al. 2014 ). Each tercile is 33% of the dataset. Boxplots show the distribution of mean temperature of the warmest month in the growing season from at least three GBIF occurrences. Asterisks indicate a significant difference based on t -tests on log-transformed data. Asterisks in parentheses indicate a significant difference based on PGLS analyses. *, p < 0.05; **, p < 0.01; ***, p < 0.001 Discussion Convergent modification of the rubisco large subunit gene Focusing on rbcL , we replicated analyses in wood ferns, sea lavenders, pines, and viburnums and found pervasive positive selection and/or convergent evolution, suggesting that adaptive evolution occurred ( Table 1 ; Table 2 ). A small fraction of sites (95, 226, 251, 255, 262, 340, 375, and 449) were consistently under positive selection in multiple genera ( Table 1 ; Figure 2 ). However, the disparate patterns of sites targeted by positive selection highlight the lineage-specific effects and the significance of evolutionary history in the evolution of rubisco ( Bouvier et al., 2021 ). Viburnums accounted for over 40% of those amino acid sites identified in this study, which may reflect its repeated shifts into tropical and cold temperate biomes ( Table S 4 ; Landis et al., 2021 ). Moreover, although tests of positive selection are routinely performed for rbcL , another indicator of adaptive evolution is the occurrence convergent amino acid substitutions on multiple branches of a phylogeny. Convergent or parallel evolution in rubisco has been shown previously by identifying multiple branches with the same derived amino acid ( Kapralov and Filatov, 2007 ; Christin et al., 2008 ). However, introgression through plastid capture can be a confounding factor in rbcL ( Yao et al., 2019 ). We used a recently developed method to account for these errors ( Fukushima and Pollock, 2023 ) and show that a high level of convergence is not due to introgression. These results support the hypothesis that rbcL has been repeatedly modified in a manner that would indicate adaptive evolution. Interactions with the solvent, not hydrogen bonds, best explain protein stability changes Unlike LDH or MDH, rubisco is constrained by dual substrate specificity, such that stability, turnover rate, and CO 2 specificity must be fine-tuned. Structural changes can be subtle; a five-fold difference in S c/o can result from a free energy difference smaller than that of a hydrogen bond ( Spreitzer and Salvucci, 2002 ). In lineages examined here, 11 and 3 unique substitutions added and removed oxygen, nitrogen, or sulfur atoms in the amino acid side chains, respectively, which we interpreted as a change in hydrogen bond ability ( Figure 4 ; Figure 5 ). Six and three of those substitutions were predicted to add or remove interactions within 3.7 Å, respectively ( Figure 4A,B,E,H,K ; Figure 5 ). These interactions are hypothesized to be hydrogen bonds, water bridges, or disulfide bridges, which would stabilize the protein. An increase in polarity is generally hypothesized to stabilize metabolic enzymes in poikilotherms through hydrogen bonds, which typically increase stability by 3–9 kcal mol -1 ( Fersht, 1995 ; Dong and Somero, 2009 ). In contrast to this hypothesis, the increases in polarity were predicted to decreased holoenzyme stability on average by 1.8 kcal mol -1 per large subunit ( Figure 6C ). Interactions of hydrophobic residues with the solvent appeared to destabilize the rubisco holoenzyme moderately ( Figure 6B ). This pattern diverges from the patterns observed in fish and mollusc LDH and MDH ( Dong and Somero, 2009 ; Fields et al., 2015 ). We propose that it is essential to consider the stabilizing effect of water bridges and especially those involving waters buried in the protein ( Park and Saven, 2005 ). It is also worth noting that all modelling was performed while keeping RbcS unchanged to isolate the effects of amino acid substitutions in RbcL. An alternative explanation is, therefore, that changes to substitutions in RbcL is coordinated with RbcS. We suggest further investigation. Rubisco large subunit evolution is associated with ecology and ecological shifts Divergence in positive selection between species in ancestral and derived biomes was tested in viburnums using previously published biome scorings ( Landis et al., 2021 ). Clade model analyses with CmC revealed that transitions to cold temperate and cloud forest clades were associated with stronger positive selection, while tropical clades experienced weakened positive selection ( Table 2 ). This is interesting because transitions into cold temperate forests in Viburnum often involve the evolution of deciduousness ( Edwards et al., 2017 ), buffering differences in leaf environment. Additionally, the cloud forest clade has a strikingly high d N / d S value of 161, which may reflect the unique, island-like nature of its disjunct cloud forest habitats ( Donoghue et al., 2022 ). CmD showed similar trends but did not fit the data better than its null. Additionally, the correlation of amino acid transition rates with those of biomes was assessed using phylogenetic comparative methods ( Pagel, 1994 ). Focusing on transition rates rather than simply correlating traits allows for inferences about trait dependencies ( Pagel, 1994 p. 199). While protein modifications have been linked to ecological shifts in animal systems ( Dong and Somero, 2009 ; Schott et al., 2014 ; Dungan et al., 2016 ), fewer studies have done so for rubisco (but see Hermida-Carrera et al., 2017 ). A large fraction of amino acid shifts in viburnums were associated with symmetrical transition rates, implying no strong preference for a particular biome ( Table S 4 ). In the substitutions that do correlate with biome shifts, the analyses suggest that RbcL evolution can either facilitate or inhibit transitions into new habitats, but that these shifts are not necessary or rapid ( Figure 3 ). The data therefore suggest a role of rubisco evolution in facilitating habitat shifts in “replicated bursts” ( Maddison and FitzJohn, 2015 ), which is when one trait is sufficient but not necessary for influencing the evolution of another trait. In short, the data demonstrate that RbcL evolution influenced biome shifts in Viburnum but did not fully explain them. Using an estimation of the growing season temperature (see Materials and Methods), we found that rubiscos with higher predicted stability occurred in warmer climates. This relationship parallels documented relationships in Δ G of unfolding with melting temperature, which in turn correlates with the temperature of maximal stability ( Rees and Robertson, 2001 ). Although we did not collect rubisco kinetics data, our combined analyses support a hypothesis that proteins from different species are fine-tuned for their environment, such that they perform similarly in their respective operating conditions. The weakened correlations when using PGLS for comparisons highlights the role of phylogenetic constraint in rubisco evolution ( Bouvier et al., 2021 ). The 2.4 °C difference between the most and least stable terciles would correspond to a 21% difference in rubisco activity with a Q 10 of 2.2 ( Cen and Sage, 2005 ). Rubisco, rubisco activase, and chaperones as an evolutionary module: phylogenetic considerations What might explain phylogenetic patterns in amino acid substitutions, holoenzyme stability, and climate? As argued by ( Bouvier et al., 2021 ), complementarity between rubisco subunits, rubisco activase, and chaperones could influence rubisco evolution, a pattern we anticipate would manifest as part of the phylogenetic signal. We suggest that rubisco subunits and related proteins act as an evolutionary module: a semi-independent set of traits that coevolve due to functional integration. In modular systems, traits evolve together because they contribute to a shared function, yet remain relatively decoupled from traits outside the module ( Wagner, 1996 ). For rubisco, although catalytic activity is localized to the large subunit, its folding, assembly, activation, and overall function are dependent on interactions with the small subunit and a suite of chaperones and auxiliary factors. Modular organization can arise because characters contributing to organismal function often evolve in coordination. This functional interdependence likely constrains evolutionary trajectories, to varying degrees, producing coordinated shifts in protein-protein and protein-solvent interactions across lineages. Evidence from recombinant expression supports this view that rubisco and related proteins function as a module. The chaperones Raf1, Raf2, RbcX, BSD2, and the Cpn60/20 chaperonin complex are required for a functional holoenzyme, defining a minimal set of proteins required for assembly. RbcL and RbcS do not assemble without chaperones, although expressing the cognate Raf1 partially rescues assembly ( Whitney et al., 2015 ; Aigner et al., 2017 ). Beyond assembly, impairment of rubisco activase expression decreases the CO 2 assimilation rate, including in non-steady state conditions ( Jiang et al., 1994 ; Hammond et al., 1998 ; Taylor et al., 2022 ), and the expression of specific rubisco activase orthologs are required for proper function of rubisco ( Milward, 2018 ; Gunn et al., 2020 ). Within the holoenzyme, the small subunit and its isoforms can have a considerable effect on rubisco activity in variable environments, further underscoring the enzyme’s sensitivity to its broader molecular context ( Wostrikoff and Stern, 2007 ; Ishikawa et al., 2011 ; Lin et al., 2020 ; Cavanagh et al., 2022 ). Patterns of protein integration are not unique to rubisco; in viral proteins, for instance, host chaperones can constrain the regions of the protein that is evolutionarily accessible, thereby shaping its evolutionary trajectory ( Yoon et al., 2022 ). Such constraints likely reflect the degree of integration among molecular components. Importantly, the strength of modular integration is itself variable, and this variation can influence how modules evolve over time ( Wagner, 1996 ; Wagner et al., 2007 ). The presence of multiple isoforms (and differing copy numbers) of RbcS and activase ( Salvucci et al., 1987 ; Spreitzer, 2003 ) suggests that integration is variable in rubisco, such that RbcL may not always reflect modifications other interacting enzymes. Deep-time divergences in rubisco and likely elements beyond RbcL may be linked to deviations from the canonical k cat c and K c trade-off when looking at large data sets ( Flamholz et al., 2019 ), such as in hornworts vs. other land plants ( Oh et al., 2024 ). Given the conformational constraints that RbcS imposes on RbcL ( Amritkar et al., 2025 ), it is worth considering how RbcS could reciprocally influence RbcL. This highlights how the functional context of RbcL has evolved, such that relying solely on RbcL to interpret kinetic variation may overlook key co-evolutionary shifts ( Bouvier et al., 2021 ). As a result, evolutionary history within the green lineage can obscure the drivers of rubisco adaptation if other components of the module are not considered. We suggest that treating rubisco and its interacting partners as an evolutionary module offers a valuable framework for understanding rubisco evolution. To better examine rubisco evolution in nature, a model clade approach may minimize confounding effects of deep phylogenetic divergences ( Donoghue and Edwards, 2019 ; Mabry et al., 2024 ). This strategy is valuable because photosynthesis is the outcome of nested processes that each can become rate-limiting ( Farquhar et al., 1980 ). As such, the rubisco module cannot be meaningfully removed from the broader organismal context; factors such as photochemistry, anatomy, hydraulics, and phenology affect the operating conditions of rubisco and in turn reflect the evolutionary history of plant lineages. Thus, dense and comprehensive sampling in a phylogenetic context (e.g., small subunit, activase, and chaperones) and the temperature response of rubisco kinetics in a well-studied model clade should provide a better picture of enablers and inhibitors of rubisco evolution. Focusing on lineages that have independently shifted between climatic regimes allows for evolutionary comparisons of the thermal adaptations. For instance, the Dentata , Mollotinus , and Urceolata clades within Viburnum are groups that each transitioned from warm to cold environments in the past 15 Ma ( Landis et al., 2021 ) and may be suitable for characterization of rubisco kinetics and its thermal responses ( Orr et al., 2016 ; Sharwood et al., 2016 ; Sargent et al., 2025 ). Conclusion Our analyses across lineages indicate that patterns in RbcL evolution resemble those seen in enzymes from thermally diverse animals, suggesting shared constraints and selective pressures; however, the complexity of the rubisco system may impose deep phylogenetic constraints. Positive selection, convergent evolution, altered side chain interactions in the homology model, correlation of d N / d S and amino acid substitutions with biome in Viburnum , and weak correlation of stability with growing season temperature together suggest a role of rubisco in thermal adaptation. However, the correlation between protein stability and temperature was weakened after correcting for phylogeny, underscoring the influence of evolutionary history on these patterns. To disentangle environmental effects from phylogenetic effects, we recommend a model clade approach that enables focused comparisons within shared phylogenetic backgrounds. A more phylogenetically grounded understanding of rubisco in natural systems could then dovetail into efforts to predict and improve plant resilience. Data availability All data and code required to reproduce the analyses are deposited on GitHub (github.com/aleungplants/rubisco-evolution). Author contributions A.L. performed the research with guidance from B.S.W. and R.F.S. Supplementary material View this table: View inline View popup Download powerpoint Table S 1. Results of substitution saturation assessments in DAMBE v7 ( Xia, 2018 ). Observed indices of substitution saturation ( I ss ) are compared to critical values for symmetrical ( I ss.cSym ) and asymmetrical ( I ss.cAsym ) topologies. t , df , and p report the t-statistic, degrees of freedom, and significance for testing whether I ss is significantly lower than I ss.c . A significant I ss < I ss.c indicates little saturation. View this table: View inline View popup Table S 2. Change in free energy of folding (ΔΔ G ) of the rubisco holoenzyme for each amino acid substitutions at sites under positive selection and/or convergent evolution. ΔΔ G values are the mean from ten simulations in FoldX 5.1 for the whole holoenzyme. Contributions from hydrogen bonds shown in Fig. 5 were summed from the backbone and sidechain hydrogen bond contributions. SD, standard deviation of the total ΔΔ G from ten simulations. View this table: View inline View popup Table S 3. Error-corrected convergent amino acid substitutions in each genus ( Fukushima and Pollock, 2023 ). Convergent substitutions are defined as a transition from any amino acid to the same amino acid in multiple branches of the same phylogeny. Branch numbers refer to distinct branches on which each substitution in parallel (see Data Availability for tree files; Sessa et al., 2017 ; Jin et al., 2021 ; Koutroumpa et al., 2021 ; Landis et al., 2021 ). View this table: View inline View popup Download powerpoint Table S 4. Tests for correlation of amino acid substitution with biome shifts in Viburnum . For each site, models were fit to test independent evolution of amino acid and biome, amino acid transitions dependent on biome transitions, and biome transitions dependent on amino acid transitions. The model with the lower Akaike Information Criterion (AIC) out of the latter two models was chosen to be compared to the first model with a likelihood ratio test and shown in this table. ΔAIC, difference in AIC between dependent and independent models; p , p -value of likelihood ratio test. View this table: View inline View popup Download powerpoint Table S 5. Summary of structure, stability, and phylogenetic analyses for all sites under positive selection or convergent evolution. *, positive selection or convergent evolution with posterior probability > 0.95; +, increase in H-bond ability or increase in stability; –, decrease in H-bond ability or decrease in stability. ND, not determined because the amino acid site is not represented in the crystal structure. Download figure Open in new tab Figure S 1. Correlations between growing season temperature and change in free energy of folding (ΔΔ G fold ) for each genus ( n = 93, 110, 87, and 100, respectively, for Dryopteris , Limonium , Pinus , and Viburnum ). The solid line is a linear regression, and the grey shaded area represents the 95% confidence interval. Acknowledgements The authors thank Drs. Emily Sessa and Konstantina Koutroumpa for providing the Dryopteris and Limonium tree files and giving permission to provide these files in the supplemental data. The research was funded by a Queen Elizabeth II/Charles E. Eckenwalder Scholarship in Science and Technology to A.L. and NSERC Discovery Grants to B.S.W.C. and R.F.S. Funder Information Declared Natural Sciences and Engineering Research Council, https://ror.org/01h531d29 References ↵ Aigner , H. , R. H. Wilson , A. Bracher , L. Calisse , J. Y. Bhat , F. U. Hartl , and M. Hayer-Hartl . 2017 . Plant RuBisCo assembly in E. coli with five chloroplast chaperones including BSD2 . Science 358 : 1272 – 1278 . OpenUrl Abstract / FREE Full Text ↵ Amritkar , K. , B. Cuevas-Zuviría , and B. Kaçar . 2025 . Evolutionary dynamics of RuBisCO: emergence of the small subunit and its impact through time J. Echave [ed.] ,. Molecular Biology and Evolution 42 : msae268 . OpenUrl PubMed ↵ Andersson , I . 1996 . Large structures at high resolution: the 1.6 Å crystal structure of spinach ribulose-1,5-bisphosphate carboxylase/oxygenase complexed with 2-carboxyarabinitol bisphosphate . Journal of Molecular Biology 259 : 160 – 174 . OpenUrl CrossRef PubMed Web of Science ↵ Bielawski , J. P. , and Z. Yang . 2004 . A maximum likelihood method for detecting functional divergence at individual codon sites, with application to gene family evolution . Journal of Molecular Evolution 59 . ↵ Bouvier , J. W. , D. M. Emms , and S. Kelly . 2024 . Rubisco is evolving for improved catalytic efficiency and CO 2 assimilation in plants . Proceedings of the National Academy of Sciences 121 : e2321050121 . OpenUrl CrossRef PubMed ↵ Bouvier , J. W. , D. M. Emms , T. Rhodes , J. S. Bolton , A. Brasnett , A. Eddershaw , J. R. Nielsen , et al. 2021 . Rubisco adaptation is more limited by phylogenetic constraint than by catalytic trade-off J. Echave [ed.] ,. Molecular Biology and Evolution 38 : 2880 – 2896 . OpenUrl CrossRef PubMed ↵ Bouvier , J. W. , and S. Kelly . 2023 . Response to Tcherkez and Farquhar: Rubisco adaptation is more limited by phylogenetic constraint than by catalytic trade-off . Journal of Plant Physiology 287 : 154021 . OpenUrl CrossRef PubMed ↵ Bowes , G. , and W. L. Ogren . 1972 . Oxygen inhibition and other properties of soybean ribulose 1,5-diphosphate carboxylase . Journal of Biological Chemistry 247 : 2171 – 2176 . OpenUrl Abstract / FREE Full Text ↵ Bullock , T. H . 1955 . Compensation for temperature in the metabolism and activity of poikilotherms . Biological Reviews 30 : 311 – 342 . OpenUrl CrossRef Web of Science ↵ Buß , O. , J. Rudat , and K. Ochsenreither . 2018 . FoldX as protein engineering tool: better than random based approaches? Computational and Structural Biotechnology Journal 16 : 25 – 33 . OpenUrl ↵ Capó-Bauçà , S. , S. Whitney , C. Iñiguez , O. Serrano , T. Rhodes , and J. Galmés . 2023 . The trajectory in catalytic evolution of Rubisco in Posidonia seagrass species differs from terrestrial plants . Plant Physiology 191 : 946 – 956 . OpenUrl PubMed ↵ Cavanagh , A. P. , R. Slattery , and D. S. Kubien . 2022 . Temperature-induced changes in Arabidopsis Rubisco activity and isoform expression R. Sharwood [ed.] ,. Journal of Experimental Botany: erac 379 . ↵ Cen , Y.-P. , and R. F. Sage . 2005 . The regulation of rubisco activity in response to variation in temperature and atmospheric CO 2 partial pressure in sweet potato . Plant Physiology 139 : 979 – 990 . OpenUrl Abstract / FREE Full Text ↵ Chang , B. S. W. , J. Du , C. J. Weadick , J. Müller , C. Bickelmann , D. D. Yu , and J. M. Morrow . 2012 . The future of codon models in studies of molecular function: ancestral reconstruction and clade models of functional divergence . In G. M. Cannarozzi, and A. Schneider [eds.], Codon Evolution , 145 – 163 . Oxford University Press . ↵ Christin , P.-A. , N. Salamin , A. M. Muasya , E. H. Roalson , F. Russier , and G. Besnard . 2008 . Evolutionary switch and genetic convergence on rbcL following the evolution of C 4 photosynthesis . Molecular Biology and Evolution 25 : 2361 – 2368 . OpenUrl CrossRef PubMed Web of Science ↵ Delgado , J. , L. G. Radusky , D. Cianferoni , and L. Serrano . 2019 . FoldX 5.0: working with RNA, small molecules and a new graphical interface . Bioinformatics 35 : 4168 – 4169 . OpenUrl CrossRef PubMed ↵ Dong , Y. , M. Liao , X. Meng , and G. N. Somero . 2018 . Structural flexibility and protein adaptation to temperature: Molecular dynamics analysis of malate dehydrogenases of marine molluscs . Proceedings of the National Academy of Sciences 115 : 1274 – 1279 . OpenUrl Abstract / FREE Full Text ↵ Dong , Y. , and G. N. Somero . 2009 . Temperature adaptation of cytosolic malate dehydrogenases of limpets(genus Lottia ): differences in stability and function due to minor changes in sequence correlate with biogeographic and vertical distributions . Journal of Experimental Biology 212 : 169 – 177 . OpenUrl Abstract / FREE Full Text ↵ Donoghue , M. J. , D. A. R. Eaton , C. A. Maya-Lastra , M. J. Landis , P. W. Sweeney , M. E. Olson , N. I. Cacho , et al. 2022 . Replicated radiation of a plant clade along a cloud forest archipelago . Nature Ecology & Evolution 6 : 1318 – 1329 . OpenUrl PubMed ↵ Donoghue , M. J. , and E. J. Edwards . 2019 . Model clades are vital for comparative biology, and ascertainment bias is not a problem in practice: a response to Beaulieu and O’Meara (2018) . American Journal of Botany 106 : 327 – 330 . OpenUrl CrossRef PubMed ↵ Dungan , S. Z. , A. Kosyakov , and B. S. W. Chang . 2016 . Spectral tuning of killer whale (Orcinus orca) rhodopsin: evidence for positive selection and functional adaptation in a cetacean visual pigment . Molecular Biology and Evolution 33 : 323 – 336 . OpenUrl CrossRef PubMed ↵ Edgar , R. C . 2022 . Muscle5: High-accuracy alignment ensembles enable unbiased assessments of sequence homology and phylogeny . Nature Communications 13 : 6968 . OpenUrl PubMed ↵ Edwards , E. J. , D. S. Chatelet , B.-C. Chen , J. Y. Ong , S. Tagane , H. Kanemitsu , K. Tagawa , et al. 2017 . Convergence, consilience, and the evolution of temperate deciduous forests . The American Naturalist 190 : S87 – S104 . OpenUrl CrossRef PubMed ↵ Farquhar , G. D. , S. Von Caemmerer , and J. A. Berry . 1980 . A biochemical model of photosynthetic CO 2 assimilation in leaves of C 3 species . Planta 149 : 78 – 90 . OpenUrl CrossRef PubMed Web of Science ↵ Felsenstein , J . 1985 . Phylogenies and the comparative method . The American Naturalist 125 : 1 – 15 . OpenUrl CrossRef Web of Science ↵ Fersht , A . 1995 . Enzyme structure and mechanism. 2. ed., 13. print. Freeman, New York . ↵ Fick , S. E. , and R. J. Hijmans . 2017 . WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas . International Journal of Climatology 37 : 4302 – 4315 . OpenUrl CrossRef ↵ Fields , P. A . 2001 . Protein function at thermal extremes: balancing stability and flexibility . Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology 129 : 417 – 431 . OpenUrl CrossRef ↵ Fields , P. A. , Y. Dong , X. Meng , and G. N. Somero . 2015 . Adaptations of protein structure and function to temperature: there is more than one way to ‘skin a cat’ J. E. Podrabsky, J. H. Stillman, and L. Tomanek [eds.] ,. Journal of Experimental Biology 218 : 1801 – 1811 . OpenUrl Abstract / FREE Full Text ↵ Fields , P. A. , and D. E. Houseman . 2004 . Decreases in activation energy and substrate affinity in cold-adapted A 4 -lactate dehydrogenase: evidence from the Antarctic notothenioid fish Chaenocephalus aceratus . Molecular Biology and Evolution 21 : 2246 – 2255 . OpenUrl CrossRef PubMed Web of Science ↵ Fields , P. A. , and G. N. Somero . 1998 . Hot spots in cold adaptation: Localized increases in conformational flexibility in lactate dehydrogenase A4 orthologs of Antarctic notothenioid fishes . Proceedings of the National Academy of Sciences 95 : 11476 – 11481 . OpenUrl Abstract / FREE Full Text ↵ Flamholz , A. I. , N. Prywes , U. Moran , D. Davidi , Y. M. Bar-On , L. M. Oltrogge , R. Alves , et al. 2019 . Revisiting trade-offs between rubisco kinetic parameters . Biochemistry 58 : 3365 – 3376 . OpenUrl CrossRef PubMed ↵ Fukushima , K. , and D. D. Pollock . 2023 . Detecting macroevolutionary genotype–phenotype associations using error-corrected rates of protein convergence . Nature Ecology & Evolution 7 : 155 – 170 . OpenUrl PubMed ↵ Galmés , J. , P. J. Andralojc , M. V. Kapralov , J. Flexas , A. J. Keys , A. Molins , M. A. J. Parry , and M. À. Conesa . 2014a . Environmentally driven evolution of Rubisco and improved photosynthesis and growth within the C 3 genus Limonium (Plumbaginaceae) . New Phytologist 203 : 989 – 999 . OpenUrl CrossRef PubMed ↵ Galmés , J. , J. Flexas , A. J. Keys , J. Cifre , R. A. C. Mitchell , P. J. Madgwick , R. P. Haslam , et al. 2005 . Rubisco specificity factor tends to be larger in plant species from drier habitats and in species with persistent leaves . Plant, Cell and Environment 28 : 571 – 579 . OpenUrl CrossRef Web of Science ↵ Galmés , J. , C. Hermida-Carrera , L. Laanisto , and Ü. Niinemets . 2016 . A compendium of temperature responses of Rubisco kinetic traits: variability among and within photosynthetic groups and impacts on photosynthesis modeling . Journal of Experimental Botany 67 : 5067 – 5091 . OpenUrl CrossRef PubMed ↵ Galmés , J. , M. V. Kapralov , P. J. Andralojc , M. À. Conesa , A. J. Keys , M. A. J. Parry , and J. Flexas . 2014b . Expanding knowledge of the Rubisco kinetics variability in plant species: environmental and evolutionary trends . Plant, Cell and Environment 37 : 1989 – 2001 . OpenUrl CrossRef PubMed ↵ GBIF . 2024a . Occurrence Download - Dryopteris. 138877695. GBIF. 2024b. Occurrence Download - Limonium . 11409777. GBIF. 2024c. Occurrence Download - Pinus . 11409777 . ↵ GBIF . 2023 . Occurrence Download - Viburnum . 80522031 . ↵ Gunn , L. H. , E. Martin Avila , R. Birch , and S. M. Whitney . 2020 . The dependency of red Rubisco on its cognate activase for enhancing plant photosynthesis and growth . Proceedings of the National Academy of Sciences 117 : 25890 – 25896 . OpenUrl Abstract / FREE Full Text ↵ Hammond , E. T. , T. J. Andrews , K. A. Mott , and I. E. Woodrow . 1998 . Regulation of Rubisco activation in antisense plants of tobacco containing reduced levels of Rubisco activase . The Plant Journal 14 : 101 – 110 . OpenUrl CrossRef PubMed Web of Science ↵ He , P. , S. M. Gleason , I. J. Wright , E. Weng , H. Liu , S. Zhu , M. Lu , et al. 2020 . Growing-season temperature and precipitation are independent drivers of global variation in xylem hydraulic conductivity . Global Change Biology 26 : 1833 – 1841 . OpenUrl CrossRef PubMed ↵ Hermida-Carrera , C. , M. A. Fares , Á. Fernández , E. Gil-Pelegrín , M. V. Kapralov , A. Mir , A. Molins , et al. 2017 . Positively selected amino acid replacements within the RuBisCO enzyme of oak trees are associated with ecological adaptations T.-Y. Chiang [ed.] ,. PLOS ONE 12 : e0183970 . OpenUrl CrossRef PubMed ↵ Hermida-Carrera , C. , M. V. Kapralov , and J. Galmés . 2016 . Rubisco catalytic properties and temperature response in crops . Plant Physiology 171 : 2549 – 2561 . OpenUrl Abstract / FREE Full Text ↵ Heyduk , K. , J. J. Moreno-Villena , I. S. Gilman , P.-A. Christin , and E. J. Edwards . 2019 . The genetics of convergent evolution: insights from plant photosynthesis . Nature Reviews Genetics 20 : 485 – 493 . OpenUrl CrossRef PubMed ↵ Hochachka , P. W. , and G. N. Somero . 2002 . Biochemical adaptation: mechanism and process in physiological evolution . Oxford University Press , New York . ↵ Iñiguez , C. , S. Capó-Bauçà , Ü. Niinemets , H. Stoll , P. Aguiló-Nicolau , and J. Galmés . 2020 . Evolutionary trends in RuBisCO kinetics and their co-evolution with CO 2 concentrating mechanisms . The Plant Journal 101 : 897 – 918 . OpenUrl CrossRef PubMed ↵ Ishikawa , C. , T. Hatanaka , S. Misoo , C. Miyake , and H. Fukayama . 2011 . Functional Incorporation of sorghum small subunit increases the catalytic turnover rate of rubisco in transgenic rice . Plant Physiology 156 : 1603 – 1611 . OpenUrl Abstract / FREE Full Text ↵ Jeffrey , G. A. 1997 . An introduction to hydrogen bonding . Oxford University Press , New York Oxford . ↵ Jiang , C. , W. P. Quick , R. Alred , D. Kliebenstein , and S. R. Rodermel . 1994 . Antisense RNA inhibition of Rubisco activase expression . The Plant Journal 5 : 787 – 798 . OpenUrl CrossRef Web of Science ↵ Jin , W.-T. , D. S. Gernandt , C. Wehenkel , X.-M. Xia , X.-X. Wei , and X.-Q. Wang . 2021 . Phylogenomic and ecological analyses reveal the spatiotemporal evolution of global pines . Proceedings of the National Academy of Sciences 118 : e2022302118 . OpenUrl Abstract / FREE Full Text ↵ Jordan , D. B. , and W. L. Ogren . 1981 . Species variation in the specificity of ribulose biphosphate carboxylase/oxygenase . Nature 291 : 513 – 515 . OpenUrl CrossRef Web of Science ↵ Jordan , D. B. , and W. L. Ogren . 1984 . The CO 2 /O 2 specificity of ribulose 1,5-bisphosphate carboxylase/oxygenase: dependence on ribulose bisphosphate concentration, pH and temperature . Planta 161 : 308 – 313 . OpenUrl CrossRef PubMed Web of Science ↵ Kapralov , M. V. , and D. A. Filatov . 2007 . Widespread positive selection in the photosynthetic Rubisco enzyme . BMC Evolutionary Biology 7 : 73 . OpenUrl PubMed ↵ Khan , S. , and M. Vihinen . 2010 . Performance of protein stability predictors . Human Mutation 31 : 675 – 684 . OpenUrl CrossRef PubMed Web of Science ↵ Koutroumpa , K. , B. H. Warren , S. Theodoridis , M. Coiro , M. M. Romeiras , A. Jiménez , and E. Conti . 2021 . Geo-climatic changes and apomixis as major drivers of diversification in the Mediterranean sea lavenders ( Limonium Mill .). Frontiers in Plant Science 11 : 612258 . OpenUrl PubMed ↵ Kubien , D. S. , S. M. Whitney , P. V. Moore , and L. K. Jesson . 2008 . The biochemistry of Rubisco in Flaveria . Journal of Experimental Botany 59 : 1767 – 1777 . OpenUrl CrossRef PubMed Web of Science ↵ Landis , M. J. , D. A. R. Eaton , W. L. Clement , B. Park , E. L. Spriggs , P. W. Sweeney , E. J. Edwards , and M. J. Donoghue . 2021 . Joint phylogenetic estimation of geographic movements and biome shifts during the global diversification of Viburnum A. Antonelli [ed.] ,. Systematic Biology 70 : 67 – 85 . OpenUrl CrossRef PubMed ↵ Liao , M. , S. Zhang , G. Zhang , Y. Chu , G. N. Somero , and Y. Dong . 2017 . Heat-resistant cytosolic malate dehydrogenases (cMDHs) of thermophilic intertidal snails (genus Echinolittorina): protein underpinnings of tolerance to body temperatures reaching 55°C . Journal of Experimental Biology 220 : 2066 – 2075 . OpenUrl Abstract / FREE Full Text ↵ Lin , M. T. , W. D. Stone , V. Chaudhari , and M. R. Hanson . 2020 . Small subunits can determine enzyme kinetics of tobacco Rubisco expressed in Escherichia coli . Nature Plants 6 : 1289 – 1299 . OpenUrl PubMed ↵ Lorimer , G. H. , and T. J. Andrews . 1973 . Plant photorespiration—an inevitable consequence of the existence of atmospheric oxygen . Nature 243 : 359 – 360 . OpenUrl CrossRef Web of Science ↵ Mabry , M. E. , R. S. Abrahams , I. A. Al-Shehbaz , W. J. Baker , S. Barak , M. S. Barker , R. L. Barrett , et al. 2024 . Complementing model species with model clades . The Plant Cell 36 : 1205 – 1226 . OpenUrl PubMed ↵ Maddison , W. P. , and R. G. FitzJohn . 2015 . The unsolved challenge to phylogenetic correlation tests for categorical characters . Systematic Biology 64 : 127 – 136 . OpenUrl CrossRef PubMed ↵ Milward , S. E. 2018 . Interrogating plant Rubisco-Rubisco activase interactions. 1 vol . ↵ Minh , B. Q. , H. A. Schmidt , O. Chernomor , D. Schrempf , M. D. Woodhams , A. von Haeseler , and R. Lanfear . 2020 . IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era E. Teeling [ed.] ,. Molecular Biology and Evolution 37 : 1530 – 1534 . OpenUrl CrossRef PubMed ↵ Oh , Z. G. , T. A. Robison , D. H. Loh , W. S. L. Ang , J. Z. Y. Ng , F.-W. Li , and L. H. Gunn . 2024 . Unique biogenesis and kinetics of hornwort Rubiscos revealed by synthetic biology systems . Molecular Plant 17 : 1833 – 1849 . OpenUrl PubMed ↵ Orr , D. , A. Alcântara , M. V. Kapralov , J. Andralojc , E. Carmo-Silva , and M. A. J. Parry . 2016 . Surveying Rubisco diversity and temperature response to improve crop photosynthetic efficiency . Plant Physiology : pp. 00750 . 2016 . OpenUrl ↵ Pagel , M . 1994 . Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters . Proceedings of the Royal Society of London. Series B: Biological Sciences 255 : 37 – 45 . OpenUrl CrossRef ↵ Paradis , E. , and K. Schliep . 2019 . ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R R. Schwartz [ed.] ,. Bioinformatics 35 : 526 – 528 . OpenUrl CrossRef PubMed ↵ Park , S. , and J. G. Saven . 2005 . Statistical and molecular dynamics studies of buried waters in globular proteins. Proteins: Structure , Function, and Bioinformatics 60 : 450 – 463 . OpenUrl ↵ Rees , D. C. , and A. D. Robertson . 2001 . Some thermodynamic implications for the thermostability of proteins . Protein Science 10 : 1187 – 1194 . OpenUrl CrossRef PubMed Web of Science ↵ Revell , L. J . 2012 . phytools: an R package for phylogenetic comparative biology (and other things) . Methods in Ecology and Evolution 3 : 217 – 223 . OpenUrl ↵ Sage , R. F . 2002 . Variation in the k cat of Rubisco in C 3 and C 4 plants and some implications for photosynthetic performance at high and low temperature . Journal of Experimental Botany 53 : 609 – 620 . OpenUrl CrossRef PubMed Web of Science ↵ Salvucci , M. E. , J. M. Werneke , W. L. Ogren , and A. R. Portis . 1987 . Purification and species distribution of rubisco activase . Plant Physiology 84 : 930 – 936 . OpenUrl Abstract / FREE Full Text ↵ Sargent , D. , S. M. Whitney , K. Broughton , S. Y. Jaconis , D. T. Tissue , M. Bange , W. C. Conaty , and R. E. Sharwood . 2025 . Diversity in the Rubisco temperature response kinetics and photosynthesis among four linted cotton species E. Carmo-Silva [ed.] ,. Journal of Experimental Botany . ↵ Savir , Y. , E. Noor , R. Milo , and T. Tlusty . 2010 . Cross-species analysis traces adaptation of Rubisco toward optimality in a low-dimensional landscape . Proceedings of the National Academy of Sciences 107 : 3475 – 3480 . OpenUrl Abstract / FREE Full Text ↵ Scholander , P. F. , W. Flagg , V. Walters , and L. Irving . 1953 . Climatic adaptation in arctic and tropical poikilotherms . Physiological Zoology 26 : 67 – 92 . OpenUrl CrossRef ↵ Schott , R. K. , S. P. Refvik , F. E. Hauser , H. López-Fernández , and B. S. W. Chang . 2014 . Divergent positive selection in rhodopsin from lake and riverine cichlid fishes . Molecular Biology and Evolution 31 : 1149 – 1165 . OpenUrl CrossRef PubMed ↵ Schrodinger . 2015 . The PyMOL Molecular Graphics System, Version 2.5 . ↵ Searle , S. R. , F. M. Speed , and G. A. Milliken . 1980 . Population marginal means in the linear model: an alternative to least squares means . The American Statistician 34 : 216 – 221 . OpenUrl CrossRef Web of Science ↵ Sessa , E. B. , A. Juslén , H. Väre , and S. M. Chambers . 2017 . Into Africa: molecular phylogenetics and historical biogeography of sub-Saharan African woodferns ( Dryopteris ) . American Journal of Botany 104 : 477 – 486 . OpenUrl Abstract / FREE Full Text ↵ Sharwood , R. E. , O. Ghannoum , M. V. Kapralov , L. H. Gunn , and S. M. Whitney . 2016 . Temperature responses of Rubisco from Paniceae grasses provide opportunities for improving C 3 photosynthesis . Nature Plants 2 : 16186 . OpenUrl PubMed ↵ Spreitzer , R. J . 2003 . Role of the small subunit in ribulose-1,5-bisphosphate carboxylase/oxygenase . Archives of Biochemistry and Biophysics 414 : 141 – 149 . OpenUrl CrossRef PubMed Web of Science ↵ Spreitzer , R. J. , and M. E. Salvucci . 2002 . Rubisco: structure, regulatory interactions, and possibilities for a better enzyme . Annual Review of Plant Biology 53 : 449 – 475 . OpenUrl CrossRef PubMed Web of Science ↵ Studer , R. A. , P.-A. Christin , M. A. Williams , and C. A. Orengo . 2014 . Stability-activity tradeoffs constrain the adaptive evolution of RubisCO . Proceedings of the National Academy of Sciences 111 : 2223 – 2228 . OpenUrl Abstract / FREE Full Text ↵ Tamura , K. , G. Stecher , and S. Kumar . 2021 . MEGA11: Molecular Evolutionary Genetics Analysis version 11 F. U. Battistuzzi [ed.] ,. Molecular Biology and Evolution 38 : 3022 – 3027 . OpenUrl CrossRef PubMed ↵ Tange , O. 2023 . GNU Parallel 20231122 (’Grindavík’) . ↵ Taylor , S. H. , E. Gonzalez-Escobar , R. Page , M. A. J. Parry , S. P. Long , and E. Carmo-Silva . 2022 . Faster than expected Rubisco deactivation in shade reduces cowpea photosynthetic potential in variable light conditions . Nature Plants 8 : 118 – 124 . OpenUrl PubMed ↵ Tcherkez , G. , and G. D. Farquhar . 2021 . Rubisco catalytic adaptation is mostly driven by photosynthetic conditions – not by phylogenetic constraints . Journal of Plant Physiology 267 : 153554 . OpenUrl CrossRef PubMed ↵ Tcherkez , G. G. B. , G. D. Farquhar , and T. J. Andrews . 2006 . Despite slow catalysis and confused substrate specificity, all ribulose bisphosphate carboxylases may be nearly perfectly optimized . Proceedings of the National Academy of Sciences 103 : 7246 – 7251 . OpenUrl Abstract / FREE Full Text ↵ Wagner , G. P . 1996 . Homologues, natural kinds and the evolution of modularity . American Zoologist 36 : 36 – 43 . OpenUrl CrossRef ↵ Wagner , G. P. , M. Pavlicev , and J. M. Cheverud . 2007 . The road to modularity . Nature Reviews Genetics 8 : 921 – 931 . OpenUrl CrossRef PubMed Web of Science ↵ Weadick , C. J. , and B. S. W. Chang . 2012 . An improved likelihood ratio test for detecting site-specific functional divergence among clades of protein-coding genes . Molecular Biology and Evolution 29 : 1297 – 1300 . OpenUrl CrossRef PubMed Web of Science ↵ Whitney , S. M. , R. Birch , C. Kelso , J. L. Beck , and M. V. Kapralov . 2015 . Improving recombinant Rubisco biogenesis, plant photosynthesis and growth by coexpressing its ancillary RAF1 chaperone . Proceedings of the National Academy of Sciences 112 : 3564 – 3569 . OpenUrl Abstract / FREE Full Text ↵ Wostrikoff , K. , and D. Stern . 2007 . Rubisco large-subunit translation is autoregulated in response to its assembly state in tobacco chloroplasts . Proceedings of the National Academy of Sciences 104 : 6466 – 6471 . OpenUrl Abstract / FREE Full Text ↵ Wright , I. J. , N. Dong , V. Maire , I. C. Prentice , M. Westoby , S. Díaz , R. V. Gallagher , et al. 2017 . Global climatic drivers of leaf size . Science 357 : 917 – 921 . OpenUrl Abstract / FREE Full Text ↵ Xia , X . 2018 . DAMBE7: new and improved tools for data analysis in molecular biology and evolution S. Kumar [ed.] ,. Molecular Biology and Evolution 35 : 1550 – 1552 . OpenUrl CrossRef PubMed ↵ Yang , Z . 2007 . PAML 4: Phylogenetic Analysis by Maximum Likelihood . Molecular Biology and Evolution 24 : 1586 – 1591 . OpenUrl CrossRef PubMed Web of Science ↵ Yao , X. , Y. Tan , J. Yang , Y. Wang , R. T. Corlett , and J.-F. Manen . 2019 . Exceptionally high rates of positive selection on the rbcL gene in the genus Ilex (Aquifoliaceae) . BMC Evolutionary Biology 19 : 192 . OpenUrl PubMed ↵ Yeoh , H.-H. , M. R. Badger , and L. Watson . 1981 . Variations in kinetic properties of ribulose-1,5-bisphosphate carboxylases among plants . Plant Physiology 67 : 1151 – 1155 . OpenUrl Abstract / FREE Full Text ↵ Yeoh , H.-H. , M. R. Badger , and L. Watson . 1980 . Variations in K m (CO 2 ) of ribulose-1,5-bisphosphate carboxylase among grasses . Plant Physiology 66 : 1110 – 1112 . OpenUrl Abstract / FREE Full Text ↵ Yoon , J. , E. E. Nekongo , J. E. Patrick , T. Hui , A. M. Phillips , A. I. Ponomarenko , S. J. Hendel , et al. 2022 . The endoplasmic reticulum proteostasis network profoundly shapes the protein sequence space accessible to HIV envelope H. S. Malik [ed.] ,. PLOS Biology 20 : e3001569 . OpenUrl CrossRef PubMed ↵ Zomer , R. J. , J. Xu , and A. Trabucco . 2022 . Version 3 of the Global Aridity Index and Potential Evapotranspiration Database . Scientific Data 9 : 409 . OpenUrl PubMed View the discussion thread. Back to top Previous Next Posted October 08, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Convergence, stability, and thermal adaptation in the rubisco enzyme in plants Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Convergence, stability, and thermal adaptation in the rubisco enzyme in plants Arthur Leung , Belinda S. W. Chang , Rowan F. Sage bioRxiv 2025.10.08.681247; doi: https://doi.org/10.1101/2025.10.08.681247 Share This Article: Copy Citation Tools Convergence, stability, and thermal adaptation in the rubisco enzyme in plants Arthur Leung , Belinda S. W. Chang , Rowan F. Sage bioRxiv 2025.10.08.681247; doi: https://doi.org/10.1101/2025.10.08.681247 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Plant Biology Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41909) Biophysics (21435) Cancer Biology (18576) Cell Biology (25479) Clinical Trials (138) Developmental Biology (13366) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24301) Genetics (15598) Genomics (22482) Immunology (17726) Microbiology (40359) Molecular Biology (17162) Neuroscience (88529) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4820) Physiology (7636) Plant Biology (15125) Scientific Communication and Education (2044) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)