Patterns and Rates of Plastid rps12 Gene Evolution Inferred in a Phylogenetic Context using Plastomic Data of Ferns

The trans-splicing rps12 gene of fern plastomes (plastid genomes) exhibits a unique structure owing to its variations in intragenic exon location and intron content, and thus, it provides an excellent model system for examining the effect of plastid gene structure on rates and patterns of molecular evolution. In this study, 16 complete fern plastome sequences were newly generated via the Illumina HiSeq sequencing platform. We reconstructed the phylogeny of ferns and inferred the patterns and rates of plastid rps12 gene evolution in a phylogenetic context by combining these plastome data with those of previously published fern species. We uncovered the diversity of fern plastome evolution by characterizing the structures of these genomes and obtained a highly supported phylogenetic framework for ferns. Furthermore, our results revealed molecular evolutionary patterns that were completely different from the patterns revealed in previous studies. There were significant differences in the patterns and rates of nucleotide substitutions in both intron-containing and intron-less rps12 alleles. Rate heterogeneity between single-copy (SC) and inverted repeat (IR) exons was evident. Unexpectedly, however, IR exons exhibited significantly higher synonymous substitution rates (dS) than SC exons, a pattern that contrasts the regional effect responsible for decreased rates of nucleotide substitutions in IRs. Our results reveal that structural changes in plastid genes have important effects on evolutionary rates, and we propose possible mechanisms to explain the variations in the nucleotide substitution rates of this unusual gene.

Plastid ribosomes are ubiquitous organelles in plant cells and play a vital role in the biosynthesis of proteins. In higher plants, plastid ribosomes contain approximately 60 ribosomal proteins that are encoded in both the plastid and the nuclear genetic compartments 1 . Among the plastid-encoded ribosomal protein gene structures, rps12 is the most notable. The plastid ribosomal protein S12 encoded by the rps12 gene is a highly conserved protein located in the functional center of the 30S subunit of the ribosome 2 . In fern plastomes (plastid genomes), where rps12 is a trans-splicing gene, this gene is split into three exons by two introns in most, but not all, ferns, and one intron (intron I) is discontinuous. The first exon of the rps12 gene is generally located in the large single-copy (LSC) region, whereas the second and third exons reside in the inverted repeats (IRs); the two IR copies have identical sequences but opposite transcriptional directions. More importantly, the second intron is lacking from the rps12 gene of all species belonging to three basal fern lineages: Psilotales, Ophioglossales, and Equisetales 3,4 . Thus, two distinct rps12 gene types were identified in ferns based on the presence or absence of intron II, and type I and II genes corresponded to intron-containing and intron-less genes, respectively. Plastid genomic structural alterations, such as inversions, duplications, and gene or intron loss, are often accompanied by an increase in the rate of plastome sequence evolution 5,6 . Therefore, different copy numbers of exons and the presence or absence of introns have garnered substantial interest as an avenue to explore the evolutionary patterns of this unique gene.
Introns are highly stable components of land plant plastomes, and it is widely believed that a basic set of introns was established prior to the divergence between vascular and nonvascular plants, because this intron set is shared by different taxonomic groups (such as charophytes, bryophytes, and spermatophytes) 7 . The intron contents of fern plastomes are also highly conserved, with no gains and few losses during fern evolution 4,8 . The interspecies variation in the intron content of fern plastomes mainly reflects the complete loss of several intron-containing genes in specific lineages (e.g., intron-containing rps16 genes are absent in Psilotales, Ophioglossales, and Equisetales 4 ), or the presence of intron-less alleles in some lineages (e.g., intron-less rps12 genes are present in Psilotales, Ophioglossales, and Equisetales; Lygodium japonicum lacks rpoC1 introns that are commonly found in other ferns 8 ). Intron losses have been associated with elevated substitution rates in plastid genes, as previously described 5,9,10 , but there is an extreme bias in taxon sampling of these studies. The vast majority of these studies are based on angiosperms. Consequently, we cannot determine whether this pattern is a universal phenomenon in the plastome or an independent evolutionary event in certain lineages or genes. The plastid rps12 gene in ferns provides an excellent model system to study the effect of intron losses on evolutionary rate and test the generality of this evolutionary pattern.
Another distinctive feature of the rps12 gene in ferns is the variation in copy number among its exons. As mentioned previously, in most ferns, the first exon of rps12 is located in the LSC, whereas the second and third exons reside in the IRs. It is well known that two completely identical IR regions are prominent structural features of the plastomes of nearly all land plants. The regional effect of the evolutionary pattern of plastid genes, in which the evolutionary rate of genes located in the IRs is lower than that of genes in the single-copy (SC) regions, has been documented with extensive research [11][12][13][14][15] . The sequence identity of IRs can be maintained by copy-dependent DNA repair because when mutations are introduced into one IR copy, the other copy provides a template for error correction 13,16 , thereby suppressing the substitution rate in the IRs. Biased gene conversion, as an efficient mutation-correcting mechanism, can result in a genome with different copy regions that have different mutation rates 16,17 . There are up to 900 genomic copies in a single plastid 18 , whereas the duplicative property of the IRs provides an even greater number of copies; therefore, the frequency of gene conversion in the IRs should be higher than that in the SC regions, and this phenomenon might be responsible for a significantly lower evolutionary rate in the IRs than in the SC regions 14,16 .
However, with a growing quantity of plastome data available for comparative genomics, this hypothesis of the regional effect of the IRs has been refuted by several studies [19][20][21] . In addition to the impact of IRs on the evolutionary rate, the non-IR localization, locus-and lineage-specific also have significant effects on evolutionary rate heterogeneity in plastomes 19,[21][22][23] . Therefore, we suggest that each of these factors alone cannot sufficiently explain the rates and patterns of molecular evolution in plastid genes. For this reason, we turned our research focus to the split rps12 gene because of its location in both the SC and IR regions. The variation in the exon location and intron content of the rps12 gene in fern plastomes provides a unique opportunity to explore the effect of gene structure on sequence evolution. To thoroughly understand the pattern of rps12 gene evolution in ferns, in the current study, we employed greatly expanded taxon sampling, including 91 fern species, and undertook a broad survey to investigate the impact of structural variation on plastid gene evolution.

Results and Discussion
Organization and dynamic structural evolution of plastid genomes in ferns. Whole-genome sequencing using an Illumina HiSeq platform generated 6,794,240-23,309,670 raw reads for 16 samples. We obtained 6,026,844-21,215,174 clean reads by removing adaptors and low-quality read pairs (Table 1). Following de novo and reference-guided assembly strategies, 16 newly sequenced fern plastomes were each assembled into a single circular molecule. Both genome size and GC content were relatively conserved among all species ( www.nature.com/scientificreports www.nature.com/scientificreports/ repeats (IRs, 22,528-32,017 bp) ( Fig. 1 and Table 2). Across all sequenced ferns, there were 84-86 protein coding genes, 27-29 tRNA genes and 4 rRNA genes. The polypods and tree ferns had similar coding gene contents, with a few notable distinctions. Compared with the 84 protein coding genes inferred to be present in the ancestral plastomes of polypods, tree ferns had an additional psaM gene that was shared with most other non-polypods. The majority of the sequenced samples displayed typical fern intron contents. The rps12 gene in all the sequenced species was classified as type I, and the first exon was located in the LSC, far away from the second and third exons of rps12, which were present in two copies in the IRs.
Although the plastome structure remained relatively stable, some species showed exceptional variability in tRNA gene content ( Table 2). Several tRNA gene losses from multiple independent lineages were detected, including the loss of trnR-UCG in Platycerium bifurcatum, Tectaria decurrens, and Nephrolepis cordifolia; the loss of trnV-UAC in P. bifurcatum, D. sieboldii, and Brainea insignis; the loss of trnN-GUU in T. decurrens; and the pseudogenization of trnT-UGU in Sphaeropteris lepifera and Plagiogyria subadnata ( Table 2). All of these genes have also been parallelly lost in other polypods plastomes with the exception of the trnN. The trnN is one of the core set genes contained in IRs, generally adjacent to either ndhF or chlL at the IR/SSC borders. In ferns, IRs are generally considered to be the most stable part of the plastome because genomic rearrangement rarely occurs in these regions. However, in contrast, a recent study indicated that the IR sequences and gene contents were highly variable in polypods 24 . Our results are in accordance with the latter findings and provide additional evidence for the dynamic evolution of IR regions among closely related polypod plastomes.  Table 2. Plastome features of the sequenced species. Plus and minus signs denote genes that are present and absent, respectively, in the corresponding species. ψ represents pseudogenes. Polypods are the lineage of most derived ferns that diversified in the Cretaceous period, displaying an ecologically opportunistic response to the diversification of angiosperms 25 . The plastomes of polypods have undergone multiple complex genomic reconfigurations during fern evolution, and thus, their plastomes differ substantially from the plastomes of basal ferns (Psilotales, Ophioglossales, Marattiales, and Equisetales). Plastome evolution among polypods is considered relatively static compared with that in lineages other than polypods 26 . Surprisingly, distinct genome organizations were identified in S. yakushimensis based on inversions and IR boundary variation. A major variation in the S. yakushimensis plastome relative to the genomes of other core leptosporangiates is the altered location from the ndhF-ccsA (Fig. 2). In addition, IR expansion into the SSC resulted in a duplication of the ycf1, chlL, and chlN genes in the S. yakushimensis plastome. The extent of IR expansion in the S. yakushimensis plastomes is unprecedented among ferns. The IRs in S. yakushimensis are up to 32 kb in length ( Fig. 1 and Table 2), whereas the longest IR found previously was 29 kb in Cibotium barometz 27 . These rare genome structure variations in polypods were detected for the first time.
Phylogenetic analyses. Both the 93-and 84-taxon datasets based on 50 protein-coding genes showed consistent phylogenetic framework, differing only in the support values for some nodes ( Fig. 3 and Supplementary  Fig. S1). Along the backbone, Equisetales was uncovered as a sister clade to the remaining ferns with strong support, followed by a highly supported joint Ophioglossales + Psilotales clade, itself sister to a clade with Marattiales (Fig. 3). The phylogenetic relationships among these four basal fern orders are the most debated topics in fern phylogeny. Most previous studies using nuclear genes [28][29][30] , a combination of mitochondrial and plastid sequences 31 , several plastid genes 32,33 , and whole plastome sequences also obtained this topology 34,35 . In contrast, other phyloplastomic studies tend to support grouping Equisetales and Ophioglossales + Psilotales together as a monophyletic group and sister to the remaining ferns 36,37 . In addition, the recent phylotranscriptomic analyses have also revealed a distinct topology of relationships among the basal fern orders, showing that Equisetales is the sister group to all other ferns, whereas Marattiales and Psilotales + Ophioglossales form a monophyletic group 38 . Given that plastid genes generally evolve more slowly than the nuclear genes, these topological differences may be due to different numbers of phylogenetically informative sites contained within the diverse molecular data 39 .
The diversification of leptosporangiates occurred after that of eusporangiate ferns, within which Osmundales was the earliest-diverging lineage, and then, Gleicheniales and Hymenophyllales formed a monophyletic group sister to the remaining non-Osmundales ferns (Fig. 3). The phylogenetic position of Hymenophyllales remains debated. Most phylogenetic analyses have depicted two alternative relationships: Hymenophyllales as a sister clade to Gleicheniales + the rest of non-Osmundales leptosporangiates 25,32,33,[40][41][42][43] , or Hymenophyllales and Gleicheniales together form a clade that is sister to the remaining non-Osmundales leptosporangiates 29,38,44 . Interestingly, both topologies have been found in more recent plastid phylogeny reconstructions using different data types and partition schemes 37 . Kuo et al. 37 showed that the use of plastome organization features also fails to provide additional support for either of these two topologies. This phylogenetic uncertainty is probably due to sparse taxon sampling, data types, model selection, and tree inference methods. Although the relative positions www.nature.com/scientificreports www.nature.com/scientificreports/ of the order Hymenophyllales remains inconclusive, our results should not be ignored based on the feasibility and effectiveness of plastomes for inferring phylogenies.
As in previous studies, Schizaeales was clearly identified as the sister clade to the core leptosporangia tes 29,33,36,38,44 . In the core leptosporangiates, Salviniales and Cyatheales are successive sisters to Polypodiales, and each of these orders is clearly monophyletic, with moderate to high support ( Fig. 3 and Supplementary Fig. S1). Our results showed that the earliest-diverging clade of Polypodiales is Lindsaeaceae, followed by Pteridaceae. These clades are well supported as monophyletic. Although Dennstaedtiaceae was weakly supported as the sister lineage to the eupolypods in our study, the relationship found here is congruent with the findings of most recent studies 28,36,38,45 (Fig. 3 and Supplementary Fig. S1). Eupolypods account for well over half of the extant fern diversity, and determining their sister group has been difficult. Earlier studies reported that Pteridaceae, instead of Dennstaedtiaceae, was sister to the eupolypods; however, the relationships among these three lineages were not well-resolved because of low support 41,42 . These relationships need to be further studied to ascertain which family is the sister lineage to the eupolypods. www.nature.com/scientificreports www.nature.com/scientificreports/ impact of intron loss on the rps12 evolutionary rate. The occurrence of rps12 intron loss in three basal fern lineages is considered to be an important evolutionary event in ferns. If the phylogenetic relationships among ferns are indeed consistent with our analysis, which recovered Equisetales as a sister clade to the remaining ferns, then the intron of the rps12 gene would have been independently lost at least twice during fern evolution 4 . Research suggests that the rate of plastome sequence evolution is generally affected by structural changes 5,6 ; thus, we sought to investigate the impact of intron loss on the evolutionary rate of the rps12 genes.
Branch length in a phylogenetic tree represents an estimate of the amount of sequence divergence in the corresponding lineage, which is equal to the product of the absolute substitution rate and time 46 . For this reason, we cannot ignore the timescale over which molecular rates change in comparisons of substitution rates among taxa. Therefore, we tested for an intron effect on substitution rate changes by comparing the absolute substitution rates of the type I and type II genes. Our results showed that the trans-splicing rps12 genes in all examined species were highly conserved with a size of 372 bp, encoding a total of 123 amino acids. The differences in the rates between the two types of genes were mainly reflected in the rates of synonymous substitution (R S ). Wilcoxon rank sum tests showed that the values of R S for type II genes were significantly higher than those for type I genes (P < 0.01), whereas the values of the rates of nonsynonymous substitution (R N ) were not significantly different (P = 0.8019) ( Fig. 4; Supplementary Fig. S2 and Table S1). Moreover, in order to detect whether the selection pressures acting on the two types of genes were significantly different, we compared dN/dS ratios in a phylogenic context via a model-based approach. A likelihood ratio test (LRT) was used to compare the fits of two models: the null model, where values of dN/dS were not significantly different between type I and type II, and the alternative model, where type I had different dN/dS ratios relative to type II. Overall, the alternative model was significantly different from the null model (P < 0.05).
Currently, the most widely accepted mechanism of intron loss is a reverse transcriptase (RT)-mediated model (namely, retroprocessing), which is a process of integrating intron-less cDNA generated by reverse transcription of the corresponding mRNA into the genome by homologous recombination [47][48][49] . Another two possible mechanisms of intron loss include genomic deletion 50 and exonization of intronic sequences 51 . In these cases, the intron is often removed imprecisely, changing the intron/exon borders. Consequently, the rps12 intron loss that occurred in the plastomes of ferns is indicative of the first mechanism because the exon boundaries of the intron-less genes have been shown to be perfectly matched to that of the intron-containing genes. Previous studies have reported that the accelerated evolution of clpP1 in many seed plant lineages is associated with intron loss 9,10,52 . A possible explanation for this acceleration in clpP1 evolution is that reverse transcriptases and/or RNA polymerases have higher error rates during retroprocessing 53 . Although the rapid evolution of the clpP1 gene could also relate to a hybrid effect (e.g., duplications, indels, and pseudogenization), we found a significantly higher R S for type II genes than for type I genes in the fairly conserved gene rps12, indicating that accelerated gene evolutionary rates are correlated with loss of introns. Therefore, this phenomenon may be a prevalent genome-wide pattern. www.nature.com/scientificreports www.nature.com/scientificreports/ Furthermore, we analyzed the two well-known types of genetic mutations, transitions (TIs) and transversions (TVs), to test whether there were differences in mutation pattern between type I and type II genes. Our results showed that TIs occurred more frequently than TVs in both type I and type II genes ( Fig. 4; Supplementary  Fig. S2 and Table S2). From the perspective of natural selection, this "transition bias" phenomenon is considered to indicate that selection disfavors transversions, as transversions are more likely to alter the amino acid sequence of proteins than transitions 54 . Comparison between the TI and TV values showed that both values were significantly higher in type II than in type I genes (P < 0.001) ( Fig. 4; Supplementary Fig. S2 and Table S2). An LRT was also implemented to identify rate shifts between the type I and type II branches, and the results indicated that TI/ TV was significantly different in the type I branches than in the type II branches (P < 0.001). That is, transition and transversion events that occurred in type II genes resulted in more synonymous changes than those in type I genes. This finding may explain the significant change in R S between the type I and type II genes, whereas R N did not exhibit significant changes.
Complicated rate variation in rps12 exons. Due to a series of genomic inversions across the fern phylogeny, 84 of 93 sampled species showed the same pattern of exon distribution, which was divided into IR and SC exons. The consensus sequences of the SC and IR exons of the fern rps12 genes were 114 bp (encoding the 1st to 38th amino acids) and 232 bp in length (encoding the 39th to 123rd amino acids), respectively. Since nonsynonymous mutations are more strongly affected by selective constraints, while synonymous mutations are largely invisible to natural selection, a comparison of the synonymous substitution rates provides a better understanding of DNA sequence evolution 55 . In order to compare with previous research more intuitively, we used the same pairwise comparison method as the previous studies when investigating the difference between the rates of IR and SC exons 12,21 . Our results revealed a more complicated rate variation in the rps12 exons than expected. The dS values for the IR exon were significantly higher than those for the SC exon (P < 0.001), a pattern that contrasts the regional effect responsible for decreased rates of nucleotide substitutions in the IRs (Fig. 5; Supplementary  Table S3). The results of an LRT between the one-rate and two-rate models showed significant rate changes relative to the position of an exon in either the IR or SC region (P < 0.01).
A previous study indicated that the total substitution rates decreased after gene translocation into the IRs by using a model-based approach across the fern phylogeny 11 . Moreover, Zhu et al. 12 compared the dS values of IR and SC genes across land plants using the pairwise estimation method, including some representative fern species. Based on very limited taxon sampling, they observed that former IR genes that moved into the SC region due to genome rearrangements underwent corresponding dS accelerations. However, we discovered that the IR exons of a given gene in fern species had, in fact, significantly higher dS rates on average than the corresponding SC exons. The most substantial distinction between our study and previous studies is that we compared the IR exons of rps12 with the SC exons of the same gene rather than other SC genes. Furthermore, we employed vastly expanded taxon sampling that included not only the major fern lineages but also additional species of derived lineages. Comparing the rate changes in exons located in different regions of a given gene can effectively avoid rate heterogeneity caused by lineage effects 8 . The intragenic location changes in the rps12 exons present a unique opportunity to test the effect of IRs on evolutionary rate variation. Thus, the pattern of substitution rate inhibition in the IRs may not be applicable to the fern rps12 gene. Indeed, several studies have shown that a pattern of decreased IR substitution rates is not universally suitable among vascular plants, as increased IR substitution rates were previously observed in some genes of Pelargonium 19,21 , Plantago, and Silene 12 . Several possible mechanisms have previously been proposed to explain substitution rate increases in plastid genes. These abnormally evolving genes could result from local hypermutation potentially induced by a high level of error-prone double-strand break repair 12 . A subsequent study by Weng et al. 21 expanded taxon sampling to characterize rates and patterns of evolution in Pelargonium, and the results showed that the anomalous rate acceleration observed in Pelargonium plastomes could be explained by a mixture of locus-specific, lineage-specific, and IR-dependent effects.
Thus, the mechanisms responsible for generating substitution rate variation may be different in each case. Considering that homologous recombination could enhance copy-correction activity in the IRs via gene conversion, the most direct explanation is the reduction in recombinant activity in IRs. Furthermore, a model of recombination repair 16,56 could serve as another possible explanation for their observations of increased dS. Recombination repair is also known as post-replication repair 56 , which is a process of repairing impaired www.nature.com/scientificreports www.nature.com/scientificreports/ molecules by using undamaged molecules as donors, and if this repair is error-prone, single base pair substitutions will be generated because the repair process involves gap-filling DNA synthesis 16 . Additional copies of IR exons in a cell or organelle can act as donors to ensure the efficiency of recombination repair. Consequently, we speculate that recombinant activity and recombination repair may have larger effects than the IR on the substitution rate variation in the rps12 gene.
Furthermore, values of TI and TV for SC and IR exons were also estimated. The results showed that both the TI and TV values in the SC exons were significantly higher than those in the IR exons ( Fig. 5; Supplementary  Table S3), and an LRT rejected the one-rate model in favor of the two-rate model (P < 0.001). Likewise, comparisons of the DNA sequences of the rps12 gene showed that the base composition in the two exons was not uniform and an excess of transitional over transversional substitutions was present. One reason for the transition bias observed in these exons could be a mutational bias due to the intrinsic properties of DNA, as purines and pyrimidines have different conformational sizes.

Materials and Methods
Taxon sampling and DNA extraction. Ferns are a species-rich lineage of vascular plants, occupying a high diversity of ecological niches, and are a major component of the earth's land flora 57 . However, only 132 complete plastomes are available from GenBank. All 132 complete fern plastome sequences in the NCBI RefSeq collection as of 10 Feb 2019 were downloaded from GenBank, and the annotated coding sequences and the number of exons in rps12 were extracted from these sequences. Because the rps12 gene is highly conserved in fern plastomes and no substantial variation in the rps12 sequence was observed among congeners, only one sample was chosen from each genus (except for two samples that were selected from Equisetum) among the previously published fern plastomes to reduce redundancy in the dataset. Then, additional plastome sequences from 16 fern species were sequenced to ensure the coverage of additional fern clades. The increased sampling was undertaken to better understand the evolutionary patterns of rps12 gene in ferns. Sequenced species mainly included tree ferns (Cyatheales) and polypod ferns (Polypodiales), which contain most of the extant fern diversity. This sampling strategy resulted in 93 samples representing all 11 extant fern orders and 32 families and two plastomes of lycophytes (outgroup) (Supplementary Table S4).
In this study, fresh leaves of 16 newly sequenced fern species were sampled from Wuhan Botanical Garden, Chinese Academy of Sciences (CSA), and Fairy Lake Botanical Garden (CSA), respectively (Table 3). Samples were taken from young leaves of each plant and were either flash-frozen in liquid nitrogen or placed in paper envelopes and dried with silica gel. Genomic DNA was extracted from the silica-dried or lyophilized leaf tissue of each sample using the Plant Genomic DNA Kit (Tiangen Biotech., Beijing, China). The DNA concentration and purity assessments were performed using agarose gel electrophoresis and a NanoDrop spectrophotometer (Thermo Scientific, Carlsbad, CA, USA). Isolations with concentrations ≥ 150 ng/μl were chosen for Illumina sequencing.
Genome sequencing, assembly, and annotation. Genomic DNA was sheared using a Covaris M220 focused-ultrasonication device (Covaris Inc., MS, USA) to a mean fragment size of 300 bp. Paired-end libraries were prepared using the NEBNext Ultra DNA Library Prep Kit (New England Biolabs, Ipswich, MA) for sequencing on an Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA). Subsequently, 150 bp paired-end reads were produced with an insert size of ~300 bp. Following enrichment, we obtained the raw data for 16 species, ranging from 2.04-6.99 Gb (Table 1). Raw data were then filtered, and adaptor sequences were removed with Trimmomatic v0.33 58 . Trimming was performed from both ends of each read, removing all bases with a quality lower than Q20, and keeping only reads of 50 nt or longer. The quality-filtered reads were subjected to de novo www.nature.com/scientificreports www.nature.com/scientificreports/ assembly using Velvet v1.2.08 59 . Each assembled contig was BLAST searched against previously reported fern plastomes to identify plastid contigs. Gaps that remained in the assembled draft sequence were filled by polymerase chain reaction (PCR) using specific primers that were designed from the regions flanking the gaps.
We annotated plastomes in DOGMA 60 . Plastid genes were corrected in each draft genome by performing homology searches with BLASTX or BLASTN against the previously published fern plastomes. Genes were annotated as pseudogenes if they showed disruptions in reading frames or frameshifts caused by nontriplet insertions or deletions. Transfer RNA (tRNA) genes were verified using two online programs, ARAGORN 61 and tRNAscan-SE 62 . Plastome maps for all species were drawn using OGDRAW 63 . All fully annotated plastome sequences were deposited in NCBI GenBank under accession numbers MN623352-MN623367 (Table 1).

Phylogenetic analyses.
Phylogenetic analyses were performed on two separate protein-coding gene matrices. The first included 50 shared genes from all the sampled taxa. The second dataset included 50 gene sequences from 84 taxa. This matrix included the same genes as the first dataset except that those taxa with a single copy of rps12 were removed (Supplementary Table S4). All sequence alignments for each gene were conducted using MAFFT v7.310 64 , and trimmed to exclude poorly aligned positions using Gblocks v0.91b 65 with default settings. The trimmed alignments were concatenated in SequenceMatrix v1.0 66 . Phylogenetic trees from the concatenated dataset were estimated using the maximum likelihood (ML) method in RAxML v8.2.4 67 with the standard GTRGAMMA model by 1,000 rapid bootstrap replicates.
Estimation of evolutionary rates. Nucleotide sequences for each part exon of the rps12 gene from 91 fern species and two lycophyte outgroups, Huperzia lucidula and Isoetes flaccida, were extracted from plastomes and spliced into a complete coding region. For this step, the exon/intron boundaries of each sample must be carefully examined, and incorrectly annotated genes were manually adjusted if necessary because misalignment would result in the erroneous inference of substitutions rates. The coding sequences of the rps12 gene were aligned at the protein level using the Align Codon option of ClustalW as implemented in MEGA7 68 .
Absolute substitution rates for rps12 were calculated according to the methods that have been described previously 69,70 . In brief, divergence times for all nodes within ferns were estimated using the penalized likelihood approach in the program r8s 71 , with a time constraint of 354 million years for the crown group age of ferns 25 . The ML tree constructed from the dataset of all sampled taxa was used as the starting tree. We performed 100 bootstrap replicates to calculate standard errors for the divergence time of each node. The branch lengths of dS (number of substitutions per synonymous site) and dN (number of substitutions per nonsynonymous site) trees were computed using the codeml module in PAML v4.8 72 . The standard errors of dN and dS branches were calculated from the standard errors of their corresponding total branch lengths as reported by PAML. Absolute rates of synonymous (R S ) and nonsynonymous (R N ) substitution for each branch were calculated by dividing the synonymous/nonsynonymous branch length by the length of time corresponding to that branch.
To compare the nucleotide mutation patterns of the two gene types, the TI and TV rates for rps12 were also estimated in the phylogenetic context of ferns using the HKY85 model (the local parameters option) as implemented in HyPhy 73 . The significance of the differences in R S , R N , TI, and TV between the type I and type II branches were evaluated using Wilcoxon rank sum tests. We used a model-based method to test whether the dN/ dS changes differed between the type I and type II branches in a phylogenetic context. The phylogenetic tree of all 93 species generated using ML was used as a constraint tree. The null model, where one dN/dS ratio was fixed across ferns, was compared with an alternative model, where the type I and type II clades were allowed to have different dN/dS ratios. LRTs were performed in PAML v4.8 72 to compare the goodness of fit of the two models. We also tested for possible rate heterogeneity in the TI/TV ratios across different branches in a phylogenetic context. The LRT scheme for detecting a significant difference in TI/TV between type I and type II genes was conducted using the HKY85 model in HyPhy 73 , in which the null and alternative models specified were the same as those used when detecting dN/dS. Furthermore, to investigate whether the evolutionary rate variations in the rps12 exon corresponded to its location, that is, the IR or SC regions, we performed a parallel analysis of 84 species to determine the effects of IR location on the substitution rate of a given gene. Nine samples were excluded from this analysis because their rps12 exon sequences were all single copies (Supplementary Table S4). Based on their locations, the rps12 exon sequences were extracted from the plastomes and aligned at the protein level, and the alignments were concatenated separately. The pairwise dS rates for the IR and SC exons between the two lycophytes and each fern species were calculated using the codeml module in PAML v4.8 72 (seqtype = 1, runmode = −2, and CodonFreq = 2 in the codeml.ctl files). The ML tree inferred from 50 plastid genes of 82 taxa was used as a constraint tree. Likewise, we calculated the pairwise TI and TV of the two exons by using the HKY85 model implemented in HyPhy to compare the nucleotide mutation patterns. Wilcoxon rank sum tests were used to test the rate differences between the SC and IR exons. To investigate whether the dN/dS changes in exons corresponded to their location, on the phylogeny, a model-based rate analysis was used. LRTs were performed using the MG94xHKY85 codon model in HyPhy 73 , in which the one-rate model specified a shared dN/dS across the whole gene, whereas the two-rate model allowed the IR and SC exons to have different dN/dS values. The same setting was applied to the LRTs for TI/TV, except that a HKY85 substitution model was applied.

Data availability
The complete plastome sequences generated in current study are available in GenBank, https://www.ncbi.nlm. nih.gov/genbank/ (accession numbers are described in the text).