Adaptive evolution and co-evolution of chloroplast genomes in Pteridaceae species occupying different habitats: overlapping residues are always highly mutated

Background The evolution of protein residues depends on the mutation rates of their encoding nucleotides, but it may also be affected by co-evolution with other residues. Chloroplasts function as environmental sensors, transforming fluctuating environmental signals into different physiological responses. We reasoned that habitat diversity may affect their rate and mode of evolution, which might be evidenced in the chloroplast genome. The Pteridaceae family of ferns occupy an unusually broad range of ecological niches, which provides an ideal system for analysis. Results We conducted adaptive evolution and intra-molecular co-evolution analyses of Pteridaceae chloroplast DNAs (cpDNAs). The results indicate that the residues undergoing adaptive evolution and co-evolution were mostly independent, with only a few residues being simultaneously involved in both processes, and these overlapping residues tend to exhibit high mutations. Additionally, our data showed that Pteridaceae chloroplast genes are under purifying selection. Regardless of whether we grouped species by lineage (which corresponded with ecological niches), we determined that positively selected residues mainly target photosynthetic genes. Conclusions Our work provides evidence for the adaptive evolution of Pteridaceae cpDNAs, especially photosynthetic genes, to different habitats and sheds light on the adaptive evolution and co-evolution of proteins. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-023-04523-1.


Background
The chloroplast genome of eukaryotes has been substantially reduced due to gene transfer to the host genome and simple gene loss [1].Therefore, fewer than 5% of genes from the ancestral cyanobacteria remain in chloroplast DNA (cpDNA), mainly those encoding housekeeping and photosynthesis-associated proteins [2].Nonetheless, chloroplasts are the sites of many important cellular processes, including photosynthesis and lipid, amino acid, chlorophyll, and carotenoid biosynthesis [3,4].Moreover, chloroplasts play important roles in plant adaptation to environmental stress, such as drought [5], salinity [6], extreme temperature [7], high light [8], and heavy metal stress [9].Because chloroplasts act as the energy hubs of plant cells, their homeostasis is readily affected by environmental stress [3,10].Chloroplasts function as environmental sensors, transforming fluctuating environmental signals into physiological responses [11][12][13].For instance, excessive reactive oxygen species (ROS) in plants can trigger retrograde signaling to subcellular compartments to readjust whole-cell metabolism in order to repair or turn over damaged macromolecules [7,14].In addition, the timely expression of chloroplast genes can promote adaptation to environmental fluctuations [15].
Genes encoding proteins with important functions are usually subject to strong evolutionary selection pressure.To measure selection pressure at the molecular level, the nucleotide substitution rate can be used to reflect the changes in protein-coding sequences: non-synonymous substitutions (dN) cause amino acid changes, while synonymous substitutions (dS) do not [16,17].dN/dS (ω) is a measure of natural selection that is widely used to detect genes related to environmental adaptation; a ω value greater than 1 is commonly believed to represent positive selection acting on a gene or gene lineage; a ω equal to 1 represents neutral selection, and a ω less than 1 represents purifying selection [18][19][20].Many amino acids in functional proteins are conserved due to strong structural and functional constraints; thus, the average dN rate is rarely higher than the average dS rate [21,22].However, some genes may be positively selected under certain environmental conditions, resulting in higher ω values [23].
There may be an evolutionary association between amino acid residues of a given protein: co-evolution may occur between sequences within a protein that physically interact or are functionally relevant, so a change in residues at one site in the molecule may lead to a change in selection pressure at another site.However, amino acid residues located in crucial functional/structural regions will be subject to stronger selection restrictions because these sites may have a huge influence on the protein's function [24].Harmful mutations at such sites will be immediately eliminated from the population, but for a mutation that is actively selected for, compensation and substitution may occur to restore any fitness loss due to the mutation [25].In other words, amino acid co-evolution or correlated mutation is a phenomenon whereby a deleterious substitution at one position is compensated for by another substitution elsewhere in the protein, so that the structure and function of the protein remain stable [26][27][28][29].Notably, if the adaptive positively selected sites of a protein map to positions important for protein structure or function [30][31][32][33], and these sites show coevolutionary relationships, this may indicate their functional/structural dependence [34].Therefore, identifying the strategies employed by plants for adapting to different environments at the molecular level, along with their co-evolutionary dynamics, can shed light on the evolutionary pattern of species as well as the complex co-adaptation between residues in proteins.
Photosynthesis is particularly sensitive to environmental conditions and plays an important role in the unique niche occupied by land plants [35,36].Despite the wide variations in plant morphology and habitat throughout the plant kingdom, the structure and function of most cpDNA-encoded proteins are relatively well conserved.This conservation is largely believed to be driven by strong selection pressure related to the functional requirements of the proteins involved in photosynthesis [37].Some functional genes with key roles in photosynthesis have undergone adaptive evolution with the "radiation"-type differentiation of species.For instance, genes encoding Rubisco components evolved under positive selection in most terrestrial plant lineages [38], and adaptive evolution of rbcL has also been observed in some aquatic plant lineages [39,40], suggesting that RbcL subunits may have undergone continuous "fine-tuning" in different ecosystems.Similarly, atpB, psaB, and rbcL in Chlamydomonas sp.ICE-L show conserved adaptive evolution in extreme environments, which supports the notion that the adaptation of algae to extreme environments is related to the mode of selection of plastid proteins [41].
The Pteridaceae family of ferns contains five subfamilies, 53 genera, and an estimated 1,211 species, contributing to ~ 10% of extant leptosporangiate fern diversity [42,43].Pteridaceae have a diverse habitat distribution, ranging from wet tropical to arid regions, and they occupy an unusually broad range of ecological niches [44][45][46].Most Pteridaceae species are terrestrial and inhabit open, often rocky environments, but representatives of some Adiantum and Pteris species are frequent in forests [47].The Parkerioideae subfamily, represented by Ceratopteris, thrives in aquatic habitats, and Acrostichum is often associated with mangroves in the intertidal zone [48].Species in the Cheilanthoideae subfamily are significant components of arid terrestrial habitats [49].Vitarioid ferns, as sisters to the genus Adiantum L., are highly simplified and predominantly epiphytic species [50].This provides an ideal system for studying the adaptive evolution and intra-molecular co-evolution of cpDNA in a wide range of niches.Previous studies on Pteridaceae have mainly focused on exploring their phylogenetic relationships based on their DNA sequences; however, few studies have focused on the molecular evolution of this family, especially between species inhabiting different environments.In this study, we divided 41 species of Pteridaceae into the terrestrial clade, the epiphytic clade, and the aquatic clade and analyzed the evolutionary dynamics of cpDNA of different habitat groups using common protein-coding sequence data sets to analyze (1) the molecular mechanism of adaptive evolution of Pteridaceae cpDNAs in different groups and (2) the intramolecular co-evolution patterns of genes involved in this adaptive evolution.

Construction of a phylogenetic tree for adaptive evolution and co-evolution analyses
We constructed a phylogenetic tree using protein sequences inferred from cpDNAs of 41 Pteridaceae species (Table S1).The topologies and support values from maximum likelihood (ML) and Bayesian inference (BI) analyses were similar, and the 41 species were resolved into five highly supported "subfamily clades": Parkerioideae, Pteridoideae, Cryptogrammoideae, Cheilanthoideae, and Vittarioideae.Pteris and Adiantum were each resolved as monophyletic in this analysis, consistent with previous reports [43,[51][52][53].We categorized this system into three clades based on their respective habitats: the aquatic clade belonging to the Parkerioideae subfamily, the epiphytic clade consisting of vittarioid ferns, while the terrestrial clade composing of the others (Fig. 1).

Analysis of selection pressure
We used the basic model to detect the selection pressure range in 41 Pteridaceae species.The ω values of Pteridaceae chloroplast protein-coding genes ranged from 0.0046 to 0.4432 (Table S2), and all 76 common protein-coding genes had ω less than 0.5, indicating that the protein-coding genes of Pteridaceae chloroplasts are dominated by purifying selection.
To test whether the evolutionary pattern of the chloroplast protein-coding genes in the Pteridaceae family is related to adaptability to different habitats, we set the terrestrial, aquatic, or epiphytic ancestral clade as the foreground clade successively when using the branch model and the other two non-target clades as background clade (Specific methods are detailed in Materials and Methods).The branch model will detect positive selection in a lineage only when the average dN of all sites is higher than the average dS [21].We discarded five genes with extremely low dS values to avoid misestimating the ω value (999): accordingly, we eliminated infA from the terrestrial clade, psbJ from the epiphytic clade, and psaI, psaJ, and psbT from the aquatic clade.Using the remaining data sets, only when the terrestrial ancestral clade was used as the foreground clade did rbcL accept the M2 model, and the ω values became significant after p-value correction (ω = 0.0047, q-value = 0.0154).When the aquatic ancestral clade or the epiphytic ancestral clade was used as the foreground clade, we observed no significant positive selection after p-value correction.No matter which habitat was used as the foreground clade, the test results were mainly ω < 0.5 (Table S3), which also indicates that the Pteridaceae chloroplasts are dominated by purifying selection.

RELAX analysis
We identified a few genes with high but not statistically significant ω values in the above data.For instance, the ω values of psbL in the epiphytic ancestral clade and rps12 and rpl32 in the aquatic ancestral clade were greater than 1 and far higher than their respective background clades (Fig. 2).A higher ω value may be caused by the positive selection of specific genes or the relaxation of natural selection.Based on this result, we tested the selection intensity of Pteridaceae species from different habitats.We determined that the higher ω value of psbL in the epiphytic clade may be caused by the relaxation of selection pressure, while those of rps12 and rpl32 in the aquatic clade were caused by the intensification of selective pressure.In addition, rbcL with significant purifying selection in the terrestrial clade showed significant relaxation of selective pressure (Table S3).

Analysis of residues under positive selection
As most amino acids in a functional protein are under structural and functional constraints, and adaptive evolution likely affects only a few sites at a few time points [21], we used the site model, which ignores the variation in the ω value between branches, to determine whether any residues have undergone positive selection.Using the Bayes empirical Bayes (BEB) method in the site model with residue reference sites based on Pteris arisanensis, we detected positively selected residues within 20 chloroplast protein-coding genes that were statistically significant (Table 1).
The branch-site model provided evidence of positive selection on different habitat groups of Pteridaceae.Using the BEB method, we did not detect positively selected residues in the terrestrial ancestral clade.Instead, we identified 3 positively selected residues in the epiphytic ancestral clade and 2 in the aquatic ancestral clade (P > 95%, with residue reference sites based on P. arisanensis) (Table 1).These residues are encoded by ndhI (S50), psaB (T511), and rps19 (R17) in the epiphytic ancestral clade.In the aquatic ancestral clade, they are encoded by ndhG (W121) and rpl33 (T16).These residues showed strong positive selection in the branch-site model, suggesting they may be critical for adaptive evolution in the epiphytic and aquatic ancestral clades of Pteridaceae.

Relationship between adaptive evolution and co-evolution
Having more positively selected residues makes it easier to observe the relationship between positively selected residues and co-evolved residues.Therefore, we used the union of the M2a and M8 model results in the following analyses.We analyzed the relationship between positively selected residues and co-evolved residues by predicting the tertiary structures of protein sequences.The coevolved residues were scattered throughout the protein sequences (Fig. 3), and some sites were not in a direct co-evolutionary relationship even if they had a relatively close linear distance in the primary structure.For example, in the AtpA protein sequence, the co-evolved residues L82 and T83 were both predicted to interact with the distal sites in the tertiary structure (L82 & M220, distance: 35.4 Å; T83 & K255, distance: 43.8 Å) (Table 2).In addition, the positively selected residues and co-evolved residues within the protein sequence were mostly independent of each other, with only a few overlapping sites: AtpA (L308), AtpB (H123, S252, and S474), AtpH (L80), MatK (E118), NdhF (L629), PetB (S2 and R118), and RpoB (R575) (Fig. 3).After visualizing the co-evolved residues using Jalview [54], we determined that the coevolved residues that overlapped with positively selected residues tended to be highly mutated, while the nonoverlapping co-evolved residues had fewer mutations (Fig. 4).In this study, we defined a residue site as a high mutation site when the same residue was present in less than half (The results of multiple sequence alignment for specific co-evolved and positively selected residues can be found in Figure S1 and Figure S2, respectively).We employed a chi-square test to assess the conservativeness of residues in overlapping positions.The results indicated a statistical preference for overlapping residues to be located in positions with a low conservation ( χ 2 : 52.30, p-value: 0).

Discussion
Adaptive evolutionary responses to changing environmental conditions result in accelerated evolution and the functional evolution of specific stress-response proteins that favor improved fitness in the new environment [55][56][57].In this study, we observed that all protein-coding genes in the Pteridaceae chloroplast genome had ω < 0.5, indicating that Pteridaceae cpDNA has mainly undergone purifying selection.However, only rbcL showed a ω value (0.0047) in the terrestrial ancestral clade significantly different from that in the background clade, indicating that this gene has experienced significant purifying selection.Rubisco is crucial for photosynthesis, as it not only fixes CO 2 [58,59] but also improves plant growth performance under stress [60].Truncations and mutations of the conserved N terminus of RbcL can substantially affect Rubisco activity [61].RbcL has evolved adaptively  Table 2 Information about intra-molecular co-evolved residue pairs for protein sequences containing overlapping residues in co-evolution and adaptive evolution analyses in different environments and is thought to have undergone "continuous fine-tuning" in different ecosystems [38][39][40].However, in the terrestrial clade of the Pteridaceae, rbcL is under significant purifying selection and is in a state of relaxed selection.On the other hand, the branch-site model did not detect any positively selected residues in the terrestrial ancestral clade.Therefore, we speculate that the ancestors of Pteridaceae likely grew on land and had completed adaptive evolution, so the terrestrial clade of Pteridaceae is now under purifying selection.Under purifying selection, most non-synonymous substitutions will be purified, so rbcL exhibits significant purifying selection.Distance: distance between the two alpha carbon atoms.Bold font represents overlapping residues identified by co-evolution and adaptive evolution analyses.

Table 2 (continued)
Fig. 3 Distribution of positively selected residues and co-evolved residues in the AtpA, AtpB, AtpH, MatK, NdhF, PetB, and RpoB proteins of Pteridaceae.All protein tertiary structures were predicted by homology based on P. arisanensis.Purple represents co-evolution residues, magenta represents adaptive evolution residues, and firebrick red represents both and is labeled accordingly When we applied the likelihood ratio test (LRT) and BEB test without distinguishing between lineages, 20 genes showed residues under positive selection (significant at the 95% level) (Table 1).These residues were mainly encoded within photosynthetic genes, including genes encoding subunits of ATP synthase (atpA, atpB, atpE, and atpH), the predominant site of photosynthetic flux control [62]; NADH dehydrogenase (ndhB, ndhF, ndhG, and ndhK), which is important for adaptation of the photosynthetic mechanism to abiotic stress [63]; the cytochrome b6f complex (petB); and photosystem II (psbF, psbL, and psbZ).These genes encode indispensable components of photosynthesis, and the finding that they are under positive selection indicates that photosynthetic genes play key roles in the adaptive evolution of Pteridaceae.The results of the site model also revealed positively selected genes related to plastid-encoded RNA polymerase (rpoA, rpoB, and rpoC2), which are involved in the transcription of photosynthesis-related genes [64,65]; and small and large ribosome subunits (rpl16 and rps7).The other functional genes showing positive selection were cemA, matK, and ycf3.
To further elucidate the evolutionary pattern of Pteridaceae chloroplast-encoded proteins, we conducted a coevolution analysis of the protein sequences encoded by 20 genes including positively selected residues, finding that 17 showed intra-molecular co-evolution.Among these, AtpB, NdhB, RpoB, and RpoC2 encoded abundant coevolved residues (over 50 co-evolved pairs) and formed a complex co-evolution network (Figure S3); this complexity may be due to their large size, as larger proteins have more residue interactions and more possibilities for co-evolution [25].The co-evolution approach emphasizes the role of spatial structure information in the study of protein evolution and provides strong evidence to reveal the structural and functional correlation of residue sites [66].Co-evolved residue pairs often correspond to spatially proximal residues in protein structures [67].In this study, in addition to the co-evolved pairs with relatively close linear and structural distances, we also observed many co-evolutionary relationships between distal residues at the tertiary structure level, which far outnumbered those of proximal residues (Table 2), suggesting that these co-evolved pairs may be more prone to functional co-evolution [68].
When a species is faced with an extreme environment, some residues within proteins undergo adaptive changes [30][31][32][33].To maintain protein stability, this effect may be mitigated by compensatory changes in other amino acid residues [25].Therefore, the adaptive evolution of residues in a protein and the intra-molecular co-evolution will tend to be linked.However, our results were contrary to this expectation: residues undergoing adaptive evolution and co-evolution are scattered in the three-dimensional structure of proteins, and their overlap is limited (overlapping residues/adaptively evolving residues: 10/64; overlapping residues/co-evolved residues: 10/741), indicating that the positions of these two categories of sites in the protein sequence are nearly independent.After visualizing multiple amino acid sequences, whereas adaptively evolving residues generally exhibit lower conservation (Figure S2), co-evolved residues tend to display relatively higher conservation (Figure S1).Interestingly, the co-evolved residues that overlap with adaptively evolving residues generally show low conservation: that is, the sequences of the overlapping sites were always Fig. 4 The proportion of differential residues in multiple sequence alignment of co-evolved sites, positively selected sites, and overlapping sites between the two highly mutated.On the one hand, non-synonymous substitutions are more strongly affected by natural selection than are synonymous substitutions [69].If adaptive evolution leads to an increase in the non-synonymous substitution rates, it may indicate that the species is actively adapting to the new environment.Therefore, adaptive evolution is expected to increase the rate of non-synonymous substitutions.Notably, functionally/structurally important residues are often subject to strong selection constraints and thus tend to be conserved, even these are subject to mutation under extreme environmental stress.On the other hand, because the structure and function of a macromolecule depend on complex interactions between residues, such changes can alter crucial interactions with other residues.Thus, the fitness effect of a mutation at a given location may depend on the state of interacting residues, resulting in non-independent evolution (co-evolution) [70].This co-evolution should generally be quite slow, with these residues showing greater conservation [25].Moreover, the geometry of interactions in the tertiary structures of proteins may vary in time and sequence, with the residues involved potentially changing over time during the evolutionary history of a molecule [25].Thus, co-evolved residue pairs and adaptively evolved residue sites may vary in time and space.Therefore, one possible explanation for the presence of these overlapping residues is that adaptive evolution and intra-molecular co-evolution in species may be driven by distinct mechanisms, and they happen to overlap.Notably, these overlapping residues tend to exhibit lower conservation.
We also detected genes encoding positively selected residues in the epiphytic and the aquatic clades (95% level, BEB method): ndhG (W121) and rpl33 (T16) in the aquatic ancestral clade and ndhI (S50), psaB (T511), and rps19 (R17) in the epiphytic ancestral clade.The habitats of epiphytic species in a tree canopy are diverse, with some limited to smaller branches and others to larger branches or trunks [71].The aquatic A. speciosum usually grows in the shady understory of a mangrove forest and is frequently flooded by tides [72], while Ceratopteris species are often restricted to aquatic habitats, such as those near ponds and streams [73].Such shade-living species are especially vulnerable to damage to the photosynthetic machinery due to excessive light.The genes that were positively selected in these clades encode proteins crucial to this machinery: NdhG, NdhI, and PsaB are involved in photosystem I (PSI) cyclic electron transport [74][75][76], and the NDH complex appears to be particularly important for enabling the photosynthetic machinery to adapt to stress conditions [63] and/or alleviating stress [75].Rpl33 and Rps19 are ribosomal proteins, with Rpl33 being necessary to maintain sufficient cpDNA translation under cold conditions [77].However, we did not detect positively selected residues in the terrestrial ancestral clade, possibly because the Pteridaceae are descended from a terrestrial ancestor that was previously adapted to such environments.In summary, some residues encoded within Pteridaceae chloroplast genes showed positive selection, indicating that they play an important role in the adaptive evolution to different habitats.

Conclusion
Pteridaceae cpDNAs have mainly undergone purifying selection.Only rbcL in the terrestrial clade showed significant purification selection and is in a state of relaxed selection.Regardless of whether the lineage was distinguished, the positively selected residues were mostly encoded by photosynthetic genes, indicating that photosynthetic genes play an important role in the adaptive evolution of Pteridaceae species.We carried out a co-evolution analysis of 20 genes encoding adaptively evolved residues to explore the complex evolutionary pattern of proteins.These adaptively evolved sites and co-evolved sites were mostly independent, with only a few overlapping sites, and the amino acid sequences of these overlapping sites were always highly mutated.These overlapping sites may be due to the different mechanisms between adaptive evolution and co-evolution, and they may overlap in different spaces and times.Obtaining more structural/functional information for these overlapping sites will be crucial for a deeper understanding of the relationship between adaptive evolution and co-evolution.Here, we present evidence at the molecular level about the adaptive evolution of Pteridaceae cpDNAs to different habitats and provide insight into the adaptive evolution and co-evolution of proteins.

Sample collection
Fresh leaves of Pteris ensiformis, P. arisanensis, and Taenitis blechnoides were sampled from the campus of Shenzhen Fairy Lake Botanical Garden (location: E114°09′, N22°34′; altitude: 944 m), quickly frozen in liquid nitrogen, and stored at − 80 °C until use.The plant materials used in the study were identified by Ting Wang, and specimens were stored in the Herbarium of the College of Life Sciences, South China Agricultural University (Specimen numbers of P. ensiformis, P. arisanensis, and T. blechnoides are GXL20210903, GXL20210904 and GXL20210905, respectively).

Library preparation, sequencing and genome assembly
DNA was extracted from the samples using a Tiangen Plant Genome DNA Kit (Tiangen Biotech Co., Ltd., Beijing, China) according to the manufacturer's instructions.The Illumina NovaSeq6000 platform was used for sequencing.The complete chloroplast genomes of P. arisanensis and T. blechnoides were assembled using GetOrganelle [78].However, the complete chloroplast genome of P. ensiformis failed to assemble into a circle using GetOrganelle and was instead assembled using Novoplasty [79].The sequences were submitted to the National Center for Biotechnology Information (NCBI) under GenBank accession numbers OP441371 (P.arisanensis), OP743918 (P.ensiformis), and OP743919 (T.blechnoides).

DNA sequence alignment and phylogenetic analysis
The complete chloroplast genomes of 38 Pteridaceae species and Alsophila denticulata (outgroup) were downloaded from GenBank (Table S1).Combining our three newly sequenced species, a total of 41 Pteridaceae species were examined, covering all subfamilies.Seventy-six common but different protein-coding sequences of these species were retained (Table S1), MAFFT [80] was used to perform sequence alignment, and the gap area was deleted to exclude poorly aligned positions.PhyloSuite [81] was used to concatenate these sequences into a dataset for phylogenetic analysis.The maximum likelihood (ML) tree was inferred using RAxML [82], GTRGAM-MAI was selected as the nucleotide substitution model, and bootstrap values for each branch were obtained by performing 1,000 bootstrap replicates.The Bayesian inference (BI) tree was established by MrBayes [83] and was estimated by running 2,000,000 generations (Nst = 6, rates = invgamma).

Analysis of selection pressure and adaptive evolution
The 76 common non-repeating protein-coding sequences were used to build independent data sets (Table S1).The 41 Pteridaceae species were classified into three clades based on their habitat: the terrestrial, aquatic, and epiphyte clades.Before performing sequence alignment, the stop codon caused by RNA editing inside the sequences was modified, and the tail stop codon was removed.Alignment gaps and uncertainties were deleted to avoid false positives.
Selection pressure was analyzed using the CODEML [16] program.The ω values of each common proteincoding gene were calculated under the basic model (model = 0, Nsites = 0, which assumes no site-wise or branch-wise dN/dS variation).The branch model was applied by comparing the single-ratio model (M0) and two-ratio model (M2) (using the likelihood ratio test (LRT) with a χ 2 distribution); the M0 hypothesis was rejected if P < 0.05.The M0 model assumes that all Pteridaceae clades have the same ratio of non-synonymous to synonymous substitution rates; the M2 model assumes different ω between the foreground clade and background clade.The false discovery rate (FDR) correction was applied to the P values calculated above [84].
Positive selection models (M2a and M8) and null hypothesis models (M1a, M7, and M8a) provided by PAML [16] were used to perform site adaptive evolutionary analyses on a dataset of shared genes from 41 Pteridaceae species.Three sets of nested models (M2a vs. M1a, M8 vs. M7, and M8a vs. M8) were used to infer the genes experiencing positive selection.If the positive selection model significantly outperformed the null hypothesis model (p < 0.05), the genes were assumed to be under positive selection, where M8a vs. M8 resulted in lower false positive results.
The branch-site model allows different site-encoding genes to have different values of ω in different branches of the phylogenetic tree in order to test whether positive selection acts on certain sites in the foreground clade.Model A (model = 2, Nsites = 2, fixed omega = 0, omega = 2) assumes that only the foreground clade undergoes positive selection; Model A null (model = 2, Nsites = 2, fixed omega = 1, omega = 1) fixes the ω of the foreground clade in model A to 1.If Model A is significantly better than Model A null (p < 0.05), then the gene has undergone positive selection in the foreground clade.
In the above models, the LRT was compared to a χ 2 null distribution with the corresponding degrees of freedom.All sites under positive selection were retrieved using the Bayes empirical Bayes (BEB) method.

Natural selection pressure analyses
Parameter K was calculated for the data by running RELAX (implemented in HYPHY v2.5.42) [85] to test the relaxation of natural selection, the selection pressure in this context includes both purifying selection and positive selection.The K parameter is related to the value of ω as (ω background) K = (ω foreground) [85].K < 1 is indicative of relaxed natural selection and K > 1 suggests intensification in the test compared to the background.RELAX was used for LRT analysis by comparing the model with K = 1 to the model with K < 1 (or K > 1).If K < 1, P < 0.05 indicates that the test clade was under significantly relaxed selection, and if K > 1, P < 0.05 indicates that the test clade was under significantly intensified selection.

Co-evolution analysis
To further explore adaptive evolution at the molecular level, co-evolution analysis was conducted on the protein sequences encoded by genes that have undergone adaptive evolution in a specific lineage, and the structures of proteins with co-evolved residue pairs were predicted in order to understand their co-evolutionary mechanism.Co-evolution analysis was performed using the program CAPS [24], which reveals structural and functional correlations between sites by detecting whether amino acid site variants are associated.The alpha-value in the program was set to 0.01; the random sampling value was set to 1000.Protein 3D structure was predicted using Phyre2 [86] online in INTENSIVE mode and visualized using PyMOL (http://www.pymol.org).Multiple sequences of co-evolved residues were checked using Jalview [54].

Fig. 2
Fig. 2 The ω value of each retained protein-coding gene in different foreground clades under the branch model.Purple, blue, and green backgrounds represent genes related to the photosynthetic system, genetic system, and other functions respectively.* marked gene represents significance after pvalue correction.Back.: Background clade; Ter.: Terrestrial clade; Aqu.: Aquatic clade; Epi.: Epiphytic clade

Table 1
Parameter estimates for different model selection tests