Molecular evolution of psbA gene in ferns: unraveling selective pressure and co-evolutionary pattern

Background The photosynthetic oxygen-evolving photo system II (PS II) produces almost the entire oxygen in the atmosphere. This unique biochemical system comprises a functional core complex that is encoded by psbA and other genes. Unraveling the evolutionary dynamics of this gene is of particular interest owing to its direct role in oxygen production. psbA underwent gene duplication in leptosporangiates, in which both copies have been preserved since. Because gene duplication is often followed by the non-fictionalization of one of the copies and its subsequent erosion, preservation of both psbA copies pinpoint functional or regulatory specialization events. The aim of this study was to investigate the molecular evolution of psbA among fern lineages. Results We sequenced psbA , which encodes D1 protein in the core complex of PSII, in 20 species representing 8 orders of extant ferns; then we searched for selection and convolution signatures in psbA across the 11 fern orders. Collectively, our results indicate that: (1) selective constraints among D1 protein relaxed after the duplication in 4 leptosporangiate orders; (2) a handful positively selected codons were detected within species of single copy psbA, but none in duplicated ones; (3) a few sites among D1 protein were involved in co-evolution process which may intimate significant functional/structural communications between them. Conclusions The strong competition between ferns and angiosperms for light may have been the main cause for a continuous fixation of adaptive amino acid changes in psbA , in particular after its duplication. Alternatively, a single psbA copy may have undergone bursts of adaptive changes at the molecular level to overcome angiosperms competition. The strong signature of positive Darwinian selection in a major part of D1 protein is testament to this. At the same time, species own two psbA copies hardly have positive selection signals among the D1 protein coding sequences. In this study, eleven co-evolving sites have been detected via different molecules, which may be more important than others.


Background
The photosynthetic oxygen-evolving photo system II is a unique biochemical system that is capable of oxidizing water molecules [1] and is responsible for producing the almost totality of oxygen on earth [2,3]. psbA gene, along with three other chloroplast (cp) genes, namely psbB, psbC and psbD, encodes the core proteins complex in the chlo-roplasts of ferns [4][5][6][7]. Precisely, psbA gene encodes D1 protein (also known as PsbA protein) in ferns. Physical mapping and pastime sequencing unraveled a set of genome rearrangements around psbA gene in "higher" fern lineages [8][9][10]. Based on completely and partially sequenced fern pastimes [11][12][13], the psbA gene often locates either to inverted repeat regions (IRs) or to the large single copy region (LSC). Moreover, gene order around it is only found in one state in IRs (rps7-psbA-trnH), but in two possible arrangements in LSC (trnK-psbA-trnH or matK-psbA-trnH). Importantly, psbA gene duplicated when its location shifted from LSC to IRs [13]. This shift mainly affected four fern orders: Schizaeales, Salviniales, Cyatheales and Polypodiales, accounting for over 90% of extant fern diversity. Because both gene copies were preserved Angiopteris evecta (Forst.) Hoffm.

Ophioglossum vulgatum L.
Psilotum nudum (L.) Beauv.  Table S3, Additional file 2), the phylogenetic tree of 28 investigated species was reconstructed via BEAST packages. 27 nodes were numbered and tagged respectively. Length of each branch was in accordance with the estimated divergence time. Estimated parameters of each branch were showed in Table 2. Specific orders mentioned in the text were marked in phylogeny and species names within different orders were illustrated under phylogenies.

Isoetes flaccida Shuttlw. ex A. Braun
adaptations [14][15][16]. The opportunity for generating novel functions is, however, often balanced with the effects such duplications have on gene dosage [17][18][19][20]. The photosynthesis environment changed mainly as a result of the emergence of flowering plants [21], and as the open sunny ground transformed into the closed shadowy canopy. For this reason, the biodiversity of other vascular plants largely decreased, and some lineages underwent extinction [22][23][24]. Remarkably, rather than the result of being adept at holding on in the face of angiosperm domination, the leptosporangiate ferns may have the ability to capitalize upon it [25]. The underlying molecular functioning of this ability remains a mystery for evolutionary biologists, and the present study sheds some light on it. The colonization by these species of canopies angiosperm-dominated light environments sparks the idea that such extraordinary diversification may have been fuelled by functional innovation following psbA duplication. To disclose this obscurity, we performed a comprehensive analyze of the action of natural selection following the duplication in psbA.
Coevolving sites are non-independent amino acid sites during the course of protein evolution [26,27]. Among the different sites with physical or functional relationships in one protein, one mutation is likely to trigger the corresponding mutations at related sites [28]. Natural selection theory and nearly neutral theory have different explanations for the mechanism of these dynamic changes [29][30][31]. Despite the debates on the underlying mechanism, the compensatory mutations among the interdependent amino acid sites provide an important approach to understand protein structure and function [32]. The identification of coevolving amino acid sites will not only play a prominent role in the annotations of the function of D1 protein, but will also reflect the evolutionary pattern of the particular protein.
To achieve a better understanding on the evolutionary biology of psbA genes in ferns, we focused on three aims: 1) identifying gene order states around psbA genes in 11 fern orders; 2) detecting the selective constraints in consequence of psbA gene duplication; and 3) unraveling the co-evolution pattern of D1 protein among fern lineages.

Sampling and sequencing
Mainly based on Smith's system [33] and Lehtonen's recent study [34], we sampled more species in Order Polypodiales (5 species) and Cyatheales (4 species), which have significantly high current fern diversity on earth, and we selected at least one species from the 11 orders (Figures 1, 2, and Table 1). To sum up, these 27 psbA sequences represented 11 orders and 14 families (100% of fern orders and 40% of families) providing reasonable coverage of the most taxon-rich lineage. All the sequenced 20 coding regions have 1059nt, coding 353 amino acid residues. Along with their stop code, 20 currently determined psbA sequences were uploaded to GenBank (Accession: JQ684679 -JQ684698, Table 1).
Since the analysis of molecular adaptive evolution strictly rejects the termination codon (i.e. TAA, TGA, TAG) in a sequence [35], coding sequences in dataset 1, 2 and 3 showed in Additional file 1 Table S3 were applied in the investigation of selective pressures.

Phylogenetic analysis
The reconstruction of phylogenetic trees is non-sensitive to the stop codes [36][37][38]; moreover, intergenic sequence of psbA-trnH was widely accepted as good indicator for barcode of land plant [39,40]. Coding sequences along with intergenic regions of psbA-trnH in dataset 4, 5 and 6 (Additional file 2)were utilized for the reconstruction of the phylogenetic structures for better resolution ( Table 2).
In accordance with the phylogenetic structure in previous documents and present analyses ( Figure 2), we found that the combined dataset 6 might not be an accurate indicator for the reconstruction of fern phylogenetic structure under geological timescale. Even though the currently gained tree partially coincided with other relevant results from Smith's system [33] and Lehtonen's conclusion [34], we could found several disagreements.
The current phylogenetic tree along with timescale was not accurate according to other published literature [25]. i) Partial structure was not the same with the known phylogeny ( Figure 2): the Order Hymenophyllales was near to Marattiales in our study rather than Gleicheniales in well-known tree. ii) Several best estimated node ages were not precise, such as node 21 and 22. iii) Four nodes had low posterior possibilities (e.g. node 19, 20, 21 in Table 4).

Selective pressure in different psbA genes
As showed in Figures 1 and 2, we tested the action of natural selection in D1 protein using three datasets (Dataset 1, 2 and 3 in Additional file 1 Table S3) comprising both presently determined sequences (20 species in Table 1) and sequences retrieved from GenBank (7 species in Additional file 1 Table S4).
We can found subtle positive selection signals in psbA genes from fern species (Figure 3, Table 4 and Additional file 1 Table S6). However, when the species were divided into two groups as single and duplicated psbA genes, only those from single copy harbored alike signals (Table 4 and Additional file 1 Table S6). As showed in Table 4 and Additional file 1 Table S6, three codons (site 4, 91 and 155) were found as positively selected sites from dataset 1 and none from dataset 2 throughout different mathematics models.
For the detection of positive selection we used nested maximum likelihood models allowing for variation in 40 Table S3) via UCLD models, which is a summary structure of Figure 1. Right phylogeny is accepted fern cladogram based on several sources [33,34].
the ratio of non-synonymous to synonymous substitutions rates (d N /d S ) across codons, as implemented in PAML and Selecton [41,42], and three models in Datamonkey as well [43,44]. For each dataset, one Likelihood Ratio Test (LRT) was performed for d N /d S heterogeneity across codons (M0/M3 test). The Discrete model (M3) fitted the data significantly better than the one-ratio model (M0), which suggested the existence of significant variation in selective constraints among codons (Table 4,  Additional file 1 Table S5). Further, four LRTs were performed for the positively selected codons: M1a/M2a, M7/M8, M8a/M8 and M8a/Mechanistic empirical model (MEC). In order to take the differences between amino-acid replacement rates into account, MEC model employed a cpREV matrix to expand a 20 by 20 amino acid replacement rate matrix into a 61 by 61 sensecodon rate matrix [45]. By this means, a position with radical replacements will obtain a higher Ka value than a position with more moderate replacements ( Figure 3). Log-likelihood values and estimates of parameters under various models were given in Table 4, and the likelihood ratio tests in Additional file 1 Table S5. In dataset two, alternative models (M2a and M8), permitting ω > 1, failed to detect positively selected codons and showed non-significance against their null tests (M1a and M7). In M2a, though the estimated value (ω 2 = 34.2681) is  *The positions of positively selected codons were showed in brackets. p value stands for the statistical significant levels of the models, which is preset before the data manipulation. ** The contexts of the dataset were introduced in Additional file 1 Table S3.
greater than one, no codon (p 2 = 0) belongs to this kind of ω; and the estimated value for ω was not greater than 1 in M8. Meanwhile, M8 significantly outperformed M7 in both datasets; M2a outperformed M1a as well (Additional file 1 Table S5). Collectively, we can conclude three main points under nested maximum likelihood models (Table 4): 1) subtle positive selection signal was detected in D1 protein from fern species; 2) codons within D1 protein encoded by duplicated psbA genes were mainly under negative selection; 3) several codons encoded by single psbA gene might have undergone adaptive evolution.
The random-site models, SLAC, FEL, REL and MEC methods were employed to examine the adaptive evolution of the D1 protein in eleven fern orders [41][42][43][44]. For dataset one (Table 2 and Additional file 1 Table S6), one positively selected codon (91 L) has been detected via maximum likelihood methods (Table 4), none via SLAC model, one (4 T) via FEL model (p = 0.05), two (4 T and 155 T) via REL model (PP > 95%). Our results indicate that 96.99% codons are highly conserved, 2.69% evolve neutrally, and a few (0.31%) are positively selected with ω > 1 ( Table 4). In the evolution of singleton psbA gene, positive selection at a handful of codons has played an important role in the evolutionary dynamics. On the contrary, the majority of codons are under negative selection in dataset two, while a small part is under neutral selection. No positively selected codon was found in the four fern lineages (Order Schizaeales, Salviniales, Cyatheales and Polypodiales). Our results showed that most of the sites are under purifying selection, while a small part is under neutral selection. No sites in the duplicate psbA gene were found under positive selection, while 4 have statistical significance in the singleton.
The selective relaxation in duplicated psbA genes could be observed in the adaptive selection analyses based on single likelihood ancestor counting (SLAC) and fixed effects likelihood (FEL) methods. Results in Table 2 showed a decrease of selective pressure in duplicated psbA genes (Table 2) (Additional file 1  Table S6), while two amino acid sites (4 T and 155 T) were identified in dataset 1.

Co-evolutionary pattern among D1 protein
Multiple amino acid sites were involved in the coevolutionary network within D1 protein. CAPS indicated that one co-evolution pair was located between the N-terminal and α helix (Site 4 and 71) [46]. In Datamonkey [43,44], two kinds of methods for detecting coevolution (one parent and two parents) indicated that other pair has undergone co-evolution (Site 19 and 350). Moreover, the results from InterMap3D [47] indicated that ten pairs have undergone coevolution processes during their evolution dynamics (Table 5).

Phylogenetic marker
It is a great challenge to develop DNA barcodes for land plants [40,48]. Kress and Erickson (2007) recommended  Figure 1 and Figure 2, all the nodes in the phylogenetic tree were numbered and marked. ** Diverge time of each branch was estimated in BEAST software and listed accordingly. The geological time scale (era and period) were inferred from the estimated diverge time of each nodes. *** Posterior probabilities (PP) of every node were inferred via BEAST. Four of them were lower than 95%: 1) node 11; 2) node 19; 3) node 20; and 4) node 21.
rbcL gene and psbA-trnH spacer region as universal markers [49]. Our results indicated that the psbA and psbA-trnH sequences have high resolution at order and family level, but low resolution at several genera level (Figures 1 and 2, Table 5). Further investigation with extra locus will have a deeper reflection on this task.

Genome rearrangement and psbA duplication
The PCR amplification results coincide with the previous conclusion based on the completely and partially sequenced plastomes [10][11][12][13]: three kinds of gene order (trnK-psbA-trnH, matK-psbA-trnH and rps7-psbA-trnH) were found in the 20 species from 8 fern orders (Figure 1). The data of the extant complete fern cp genomes showed that the trnK-psbA-trnH and matK-psbA-trnH fragments are located in the LSC and the rps7-psbA-trnH fragment is in the IRs [10,13,50]. Since the large-scale inversion events involving the duplication of psbA gene are rare during the evolution history of fern species [8,9], one copy of psbA exists in those species of trnK-psbA-trnH and matK-psbA-trnH types and two in rps7-psbA-trnH. Their full-length psbA encoding sequences have not been determined in previous studies.

Evolutionary trajectory and survival strategy
The strong competition between ferns and angiosperms for light may have been the main cause for the evolutionary trajectory of psbA gene. Nevertheless, the occurrence of psbA duplication provided an alternative *The contexts of the dataset were introduced in Additional file 1 Table S3. Codons with PP > 95% were in bold. The phylogeny trees applied in the estimation were reconstructed throughout dataset introduced in Additional file 1 Table S3 accordingly. **Positions and posterior probability of the positively selected codons were showed in this column. The posterior probabilities were given via Bayes Empirical Bayes (BEB) method.
survival strategy versus the molecular adaptive evolution at the codons of D1 protein. Several positively selected sites were found in single psbA gene copy and none in duplicated psbA (Table 2). These results indicated that single psbA gene might have functional adaptation via the replacement at certain positions among the D1 protein while the duplicated one might have not. Conversely to the case of single gene copies, the leptosporangiate ferns with duplicate psbA gene may have acquired further fitness gain through the existence of a new transcription locus for the synthesis of D1 protein, which might directly increase the efficiency of photosynthesis by protein dosage effects [20,51,52]. The structure and function of the D1 protein are conserved among cyanobacteria, red algae and plants [53,54]. The purifying selection of psbA in Lejeuneaceae has been noticed in a recent study [55]. Regretfully, not like the elegant investigation in rbcL gene [56], the detail function of the codons in D1 protein remains obscure. Mutation experiments of D1 protein will unveil the functional importance of the positively selected positions. Moreover, the further research on the relationships between positively selected mutations among psbA gene from fern species and the biodiversity will have great impacts on the understanding of fern evolutionary biology.

M T A T L E R R E S A S L W G R F C N W T T S T E N R L Y I G W F G V L M I P T L L T A T S V F I I A F I A A P P V D I D G I P V S G S L L Y G N N I I S G A I I P T S A A I G L H F Y P I W E A
Positive Selection Negative Selection  *Three packages were applied in the current researches [46,47,73]. **Site1, Site 2: Coordinates of each of the interacting residues. ***Two kinds of molecules (two parents and one parent) were applied and gave the same pair.

Complicated intra-molecular evolution under selective pressure
Eleven co-evolving sites (Table 5) have been detected via  different molecules, and they (site 4, 19, 53, 71, 72, 92,  281, 350, 351, 352 and 353) may be more important than others in D1 protein. Point mutation experiments may have distinctive outcomes, some could cause severe functional consequences and others could result in completely undetectable change. This fact indicates that a protein is a network of interacting residues, and the core nodes in this network determine the function of the protein [57]. Future point mutation experiments aiming at the eleven sites may have butterfly effect on the protein structure and function. psbA gene, along with three other chloroplast genes, encodes the core proteins complex in the chloroplasts of ferns. However, as we have concluded before the functional adaptation of this complex might be caused by the inter-coevolution among different proteins [32]. The current intra-molecular evolution analysis may have shortcomings on predicting the complicated coevolutionary networks among the entire functional core complex. Further conclusion could not be declared without overall analysis based on the four related genes: psbA, psbB, psbC and psbD.

Conclusion
In the current research we present evidence that point to a complex adaptive process mediating the functional innovation of the D1 protein. This process involves a multiply checking of the structural and functional consequences of the fixation of functionally novel mutations and the amelioration of the effects by such mutations may have through compensatory replacement events. A serial amino acid sites are identified as co-evolution positions while significant positive selection signals are detected in the single copy psbA gene from the fern species. One hypothetical scenario is put forward: i) single copy psbA gene fern species may adapt to the newly formed living circumstance by the modification of amino acids in D1 protein; ii) by the meantime, the dosage effects of D1 protein are the possible strategy against the rising of angiosperm in psbA duplicated ones. The selective relaxation in duplicated psbA genes could be observed throughout different models. However, no evidence stands for the functional divergence in the duplicated psbA genes. Future investigation will shed new lights on this question. Although this research covers only a little of the diversity of species in fern, our sampling included all 11 orders. Our research however opens exciting new avenues that will hopefully lead to a more complete understanding of the functional novelties and dosage effects in the D1 protein among ferns.

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 [58]. 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 A 260 /A 280 ratio and A 260 absorption value.
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 [59].

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 [60]. 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 [61]. 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 [62], 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 [63]. 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 [64] and 4) Polypodiaceae clade [25].
According to authors' suggestion [37], 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 [38]. A uniform distribution is applied over the estimating of the absolute ages via the MCMC process. For each MSA, BEAST runs 6 × 10 7 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 postburn-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][43][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 N /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 [66].
Besides, another three models for detecting codons under selection are implemented on Datamonkey website [43]: SLAC, FEL and REL. SLAC originated from the Suzuki-Gojobori counting approach [67] 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 [68]. 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 [69].

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 [70], which compares the correlated variance of the evolutionary rates at 2 sites corrected by the time since the divergence of the 2 sequences [46]. 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 stepdown permutation procedure was employed to correct for multiple tests and non-independence of data [71]. 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 [72].
The mathematical models for detecting the coevolving 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 [47], 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 [73].

Additional files
Additional file 1: Table S1. Universal primers for fern psbA gene. S2 Species-specific primers for psbA gene from eight fern species. S3 Six datasets of different species and fragments. S4 Species and accession number of the retrieved data. S5 LRTs of the random-site models in PAML version 4.1*. S6 Positively selected codons determined via REL model.

Additional file 2:
Currently determined 20 sequences. The sequences started at the ATG codon of psbA, ended before the ATG codon of trnH. They included the stop codon of psbA and the intergenic region between psbA and trnH. They were applied to reconstruct the phylogenetic trees.