Sampling of plant materials
Plant materials of 20 fern species for the present investigation were collected from Wuhan Botanical Garden, Fairy Lake Botanical Garden and South China Botanical Garden, respectively (Table 1).
Isolation of total genomic DNA
For each species, three pieces of fresh leaves were collected to isolate genomic DNA with the modified CTAB protocols . Each sample was dissolved in 50 μl TE buffer. Roughly, the quality was determined by 1% agarose/TAE gel electrophoresis and the quantity was estimated via DNA Ladder (Takara). The bright sample under UV-light with right size was used as the template in PCR reactions. The absorption at 260 and 280 nm of qualified template DNA was measured using a 752 spectrophotometer. The purity and concentration was resolved and calculated by the A260/A280 ratio and A260 absorption value.
PCR amplification and DNA sequencing
PCR amplification was carried out in 100 μl volumes containing 50 mM KCl, 10 mM Tris–HCl (pH 8.0), 0.1% Triton X-100, 1.5 mM MgCl2, 0.2 mM each deoxynucleoside triphosphate, 2 U Taq DNA polymerase, 0.3 μM primer, 30 ng genomic DNA and DNA-free water. The 3-step and 2-step PCR protocols employed species-specific and universal primers, respectively. Individually, the annealing temperature in 3-step PCR reaction (Additional file 1 Table S1) of species-specific primers was calculated via Primer 3 (http://frodo.wi.mit.edu/primer3/). The thermo-cycling program was set as: 5 min at 95°C (1 cycle); 45 s at 94°C, 60s at Ta°C, 90s at 72°C (34 cycles); 10 min at 72°C (1 cycle). However, the annealing temperature was ignored in 2-step PCR reaction (Additional file 1 Table S2). The thermo-cycling program was set as: 5 min at 95°C (1 cycle); 60s at 94°C, 150 s at 60°C (35 cycles); 20 min at 72°C (1 cycle). Except as normal reactions, the genomic DNA was excluded from the reaction mix for negative control. Then the molecular weight of PCR products was verified in 1% agarose/TAE gels.
Each qualified DNA fragment amplified by the above steps was recovered and purified with a quick PCR Purification Kit (Promega), and then cloned into PMD19-T vectors (Takara). The plasmids, composed of the vectors and the DNA fragments, were transformed to Escherichia coli strain DH5α. Plasmids within positive clones were extracted and sequenced with an ABI PRISM 3730 DNA analyzer. Three clones were sequenced for each amplicon to control Taq polymerase errors. The overlapping sequences from various amplification steps were assembled as a single contig. To ascertain the contigs’ locations among cp genomes, the sequencing results were submitted to DOGMA website .
Multiple sequence alignments and best-fit nucleotide substitution model
Six multiple sequence alignments (hereafter MSAs) were established for the present investigation (Additional file 1 Table S3). Nucleotide sequences obtained experimentally (Table 1) plus those retrieved from the public databases (Additional file 1 Table S4) were aligned using the MUSCLE software . Partially fern psbA sequences from GenBank were excluded from the present study.
The best-fit nucleotide substitution model for each MSA was selected via jModeltest package . And it was also estimated via the automatic model selection tool at the Datamonkey website (http://www.datamonkey.org) for the coding regions [43, 44].
Reconstruction of phylogeny along with time-scale
Isoetes flaccida (GenBank Accession GU191333) was selected as outgroup in the phylogenetic analyses. Two nodes were chosen to constrain for a rate consistent with the known fossils: 1) Since Osmunda fossils have been described from the Upper Triassic , the Osmundaceae clade was constrained to 199.6 million years ago (Mya, Node 4); 2) According to the fossil Gleichenia, node 14 were dated to have been originated 99.6 Mya . Moreover, following previous estimates, another node 21 and 22 were separately dated to 110.5 Mya and 42.5 Mya at the beginning of the running: 3) Pteridaceae clade  and 4) Polypodiaceae clade .
According to authors’ suggestion , to avoid the misspecification of dating and taxon sampling, the empty alignment was run before the real MSAs. Then BEAST 1.7 was allowed to infer topology, branch lengths, and dates for combined datasets . A uniform distribution is applied over the estimating of the absolute ages via the MCMC process. For each MSA, BEAST runs 6 × 107 generations, saving data every 1,000 generations, producing 60,000 estimates of dates under a Yule speciation prior under the uncorrelated lognormal distributed relaxed clock model. Convergence statistics was analyzed in Tracer v 1.5, resulting in 54,000 post-burn-in trees. Before the consensus tree was graphically illustrated by Figtree v.1.3.1, TreeAnnotator v.1.6.1 was utilized to produce maximum clade credibility trees from the post-burn-in trees and to determine the 95% probability density of ages for all nodes in the phylogenetic tree.
Detection of positively and negatively selected codons
Since identification of positive/negative (non-neutral) evolution is fundamental to the understanding of the process of diversifying/purifying selection, this subject has been the focus of several decades of mathematical and computational efforts. Different scientists have developed numerous analysis models and methods, and each has its own advantages [42–44, 65]. However, the general consensus of them is that non-synonymous nucleotide substitutions (d
N), whose alternatives leading to a change in the codon and its corresponding amino acid, can be time-scaled by the number of synonymous replacements (d
S), which are nucleotide changes that only change the codon but not the amino acid and are consequently neutrally fixed and proportional to the divergence time between the sequences.
The random-site models (M0 vs. M3, M1a vs. M2a, M5 vs. null test, M7 vs. M8, and MEC vs. null test), contained in PAML package version 4.1 and Selecton version 2.2, allow the ω ratio (ω = d
S) to shift among codons within the MSA and this parameter is estimated by maximum-likelihood value via Bayesian inference approach [41, 42]. Additionally, the results from Selecton version 2.2 are visualized with seven-color scale for representing the different types of selection. To identify the statistical significant levels of the results, the LRT was conducted to compare the nested models .
Besides, another three models for detecting codons under selection are implemented on Datamonkey website : SLAC, FEL and REL. SLAC originated from the Suzuki–Gojobori counting approach  and is quite efficient on detecting non-neutral evolution in large MSAs. Less conservative than SLAC, FEL fits a site-specific d
N and d
S in the context of codon substitution models and tests whether d
N = d
S, outperforming other two models in MSAs of intermediate size. As the most powerful model of the three, REL is improved based on the Nielsen–Yang approach . Before the analysis of Natural selection, the best-fit nucleotide substitution models in the MSAs were calculated via the model-selection molecule online. Different from the posterior LRTs in the nested models, the parameter for statistical significant level (p value or Bayesian factor) was pre-set prior to the estimating processes .
Analysis of inter-dependent evolutionary sites
To understand the broad implications of the amino acid replacements in D1 proteins from fern species, the analysis of the evolutionary dependencies among sites to identify functional/structural dependencies among residues were conducted via five.
Based on a tree-ignorant strategy, CAPS outperforms the tree-aware strategy methods as reported by previous published work , which compares the correlated variance of the evolutionary rates at 2 sites corrected by the time since the divergence of the 2 sequences . The significance of the results was evaluated by randomization of pairs in the alignment, calculation of their correlation values, and comparison of the real values with the distribution of 10,000 randomly sampled values. The step-down permutation procedure was employed to correct for multiple tests and non-independence of data . An alpha value of 0.001 was applied to minimize Type I errors. The correlated variability between amino acid sites was weighted by the level of substitutions per synonymous site in order to normalize parameters by the time of sequence divergence .
The mathematical models for detecting the co-evolving residues in protein from InterMap3D and Datamonkey websites are based on a tree-aware strategy [47, 73]. Currently, three different models, namely Row and Column Weighing of Mutual Information (RCW-MI), Dependency and Mutual Information/Entropy (MI/E), were implemented at the former website , and one (Spidermonkey) at the later [43, 73]. RCW-MI and Dependency extract both the entropy dependency and the phylogenetic signal [74, 75], while the MI/E extracts the entropy dependency from the signal by dividing mutual information by the joint site’s entropy [76, 77]. Spidermonkey molecule reconstructs the substitution history of the MSAs by ML-based methods, and analyses the joint distribution of substitution events by Bayesian graphical models .