Sequence data from Qi et al. (2018) and Shen et al. (2018) were downloaded from the NCBI SRA. The quality of raw sequences was assessed with FastQC¹ and reads were trimmed of adaptors and the first five bp with Trimmomatic ver. 0.39 (Bolger et al., 2014). Transcriptomes were assembled following the methods in One Thousand Plant Transcriptomes Initiative (2019); trimmed reads were assembled with SOAPdenovoTrans ver. 1.03 (Xie et al., 2014) with a kmer size of 25 bp. To remove any plastid sequence data that may have ended up in the nuclear assemblies, available fern plastome sequences were downloaded from NCBI (accessed July, 2020). BLASTN (Altschul et al., 1990) was used to compare assemblies against the plastome database and scaffolds or contigs with significant hits (e-value < 1 × 10^–4, overlap > 300 bp, bitscore > 50) were removed from the assembly. Sequences with greater than or equal to 98% similarity were clustered with CD-HIT ver. 4.6.8 (Fu et al., 2012) to reduce transcript redundancy. Assemblies from 1KP were downloaded from the Cyverse repository for downstream analyses. Transcriptome completeness was assessed using BUSCO ver. 3.02 (Simão et al., 2015) with eukaryota odb9 lineage database (see Supplementary Appendix 1). Peptide and coding sequence (CDS) files were generated for each transcriptome using TransDecoder ver. 5.5.0 (Grabherr et al., 2011). We used Kruskal–Wallis tests by rank to determine whether assemblies from different studies were significantly different from one another for several metrics of interest, followed by pairwise comparisons using Wilcoxon rank sum tests with corrections for multiple testing using the Holm method (Holm, 1979).

Peptide and corresponding CDS files for the outgroups Amborella trichopoda, Arabidopsis thaliana, Physcomitrella patens, and Selaginella moellendorffii were downloaded from Ensembl Plants 51 (Howe et al., 2020). To ensure that all major lineages of land plants were represented, peptide and CDS files for the gymnosperm Ginkgo biloba were also downloaded from Guan et al. (2019). Peptide files from outgroups and ferns with >55% BUSCO completeness (92.6%, 239 of 258) were passed to OrthoFinder ver. 2.3.11 (Emms and Kelly, 2015, 2019) to identify orthogroups (OGs). OGs identified with OrthoFinder were filtered using the custom bash script get_orthogroups.sh to generate datasets for single- and multi-copy datasets. OGs that were single-copy and contained at least 60, 75, and 85% of the input transcriptomes were used to generate the “SCO60,” “SCO75,” and “SCO85” datasets. Multi-copy OGs that contained all transcriptomes were used to generate the “MCO” dataset.

The custom python script extract_cds.py (modified from Kasey K. Pham, pers. comm.) was used to obtain the corresponding coding sequences for each taxon in the filtered OGs. Sequences were aligned with the codon-aware alignment program MACSE ver. 2.0.4 (Ranwez et al., 2011) and gappy sites were removed with trimAl ver. 1.2 (Capella-Gutiérrez et al., 2009) by retaining sites that contained at least 50% of tips. Maximum likelihood gene trees were constructed from both the nucleotide and peptide alignments with IQTREE2 ver. 2.1.2 (Minh et al., 2020), with ModelFinder (Kalyaanamoorthy et al., 2017) and 1000 ultrafast bootstraps (Hoang et al., 2018) on 2 CPU threads where possible based on the recommendations from Shen et al. (2020). Note that trees constructed from the concatenated matrices required additional RAM and could not be run on only two CPU threads; these deviations are specified in our code. The optimal maximum likelihood gene trees for each locus were passed to Astral ver. 5.7.7 (Zhang et al., 2018) for the SCO datasets and Astral Pro ver. 1.1.2 (Zhang et al., 2020) for the MCO dataset to generate species trees under the multispecies coalescent (MSC). Concatenated alignments were used to generate a maximum likelihood species tree for the SCO datasets following the methods above, partitioned by locus (Chernomor et al., 2016). Discordance in the data was visualized using DiscoVista (Sayyari et al., 2018) and generalized Robinson-Foulds (GRF) distances between estimated species trees were calculated using the R package phangorn ver. 2.5.5 (Schliep, 2011).

Given heterogeneity among lineages (see Section “Discussion”), we compared the results from traditional models of sequence evolution to GHOST models (Crotty et al., 2020). GHOST is a mixture-model that takes a user-supplied number of components (k) and was developed to account for heterotachous evolution in datasets. For each locus in the SCO60 NT dataset, we ran IQTREE2 using the best-fitting model with classes k = 2–6 in both the linked and unlinked implementations of GHOST. The best number of classes was assessed using AIC. The optimal trees for the best-fitting number of classes under the GHOST model were passed to ASTRAL ver. 5.7.7 (Zhang et al., 2018) to construct the species tree as above. We also used a custom python script (extract_codons.py) written with the BioPython ver.1.79 module (Cock et al., 2009) to extract the first and second codon positions and third codon position from each untrimmed locus alignment in the SCO60 NT dataset. We then used trimAl on the first and second and third codon position alignments to remove sites with less than 50% tip occupancy. Gene trees were generated as above with IQTREE2 which were used to construct a species tree with ASTRAL.

Given that the phylogenetic position of two samples (Onoclea sensibilis ONSE and Plagiogyria japonica PLJA) suggested misidentifications, we extracted the longest plastid contig or scaffold from the initial assembly and used BLASTN (Altschul et al., 1990) against the nr database to determine their appropriate identification. The Onoclea sensiblis ONSE sample matched the Matteuccia struthiopteris chloroplast genome (KY427353.1) with a 98.7% identity compared to 91.26% identity to the O. sensibilis chloroplast genome (KY427354.1). The P. japonica PLJA sample matched the P. subadnata chloroplast genome (MN623362.1) with 99.25% identity. These samples have been named according to their respective hits for downstream analyses.

In a supplemental analysis, we downloaded sequence data from Dong et al. (2019) for Sphaeropteris lepifera and assembled the transcriptome as above. Coding sequences of Cystodium sorbifolium and Saccoloma campylurum from Qi et al. (2018) were provided by the author as the raw sequence files were corrupted and not available on the NCBI SRA. For these three transcriptomes, we used BLASTP (Altschul et al., 1990) to identify significant hits to the 2884 OGs in the SCO60 dataset. The contig or scaffold with the highest bit-score was extracted and nucleotide alignments, gene trees, and species trees were constructed as above.

We used SortaDate (Smith et al., 2018) to calculate the root-to-tip variance and bipartition support for each locus in the SCO60 NT dataset relative to the SCO60 NT MSC species tree. We then selected loci that were above the 90th percentile for bipartition support and below the 15th percentile for root-to-tip variance. This filtered dataset contained 99 single copy loci and was 136,137 bp in length. We generated a maximum likelihood tree with IQTREE2 as above using the SCO60 NT MSC species tree as a topological constraint to generate branch lengths relative to the number of substitutions per site. The resulting tree was then used to generate a dated phylogeny under a penalized likelihood method with treePL (Smith and O’Meara, 2012) with fifteen fossils (Supplementary Appendix 2) as calibration points along the phylogeny. To account for uncertainty in the dataset, we generated 100 bootstrap alignments from the 99 locus matrix and constructed trees from these bootstrap alignments using RAxML ver. 8.2.12 (Stamatakis, 2014). We then ran treePL on each of the bootstrap trees with three cross-validation runs and summarized them with treeAnnotator in BEAST ver. 2.5.0 (Suchard et al., 2018). We used the same methodology to date the supplemental tree with the inclusion of C. sorbifolium, S. campylurum, and S. lepifera with two additional fossil calibrations: we placed Cyathocaulis fossils (Tidwell and Nishida, 1993; Lantz et al., 1999) at 113–145 Ma for the stem Cyatheaceae as in Schuettpelz and Pryer (2009) and Du et al. (2021) and we used a stem age of 100.5–113 Ma for Lindsaeaceae based on the Regaladgo et al. (2017) as in Du et al. (2021). Although we did attempt to use MCMCTree (Yang, 2007) to estimate divergence times, due to the size of the dataset (e.g., number of loci, number of tips), Bayesian dating analyses were unfeasible.

The program wgd ver. 1.1.1 (Zwaenepoel and Van de Peer, 2019) was used to generate paralog age distributions (K_S plots) for each transcriptome. Normal mixture models were fit to the K_S distributions using the R package mixtools ver. 1.2.0 (Benaglia et al., 2009); AIC values of models with one component were compared to those of models with more than one component to determine if fits were statistically different. The R package SiZer ver. 0.1–7 (Chaudhuri and Marron, 1999; Sonderegger, 2020) was used to detect significant changes in slope (at α = 0.05). Briefly, SiZer identifies changes in slope based on changes in the first derivative in the curve (Chaudhuri and Marron, 1999). With respect to K_S plots, changes, particularly increases, would represent deviations from the background paralog distribution. Uncorrected interspecific K_S values were calculated for select species using wgd’s one_v_one feature, which generates a K_S distribution of one-to-one orthologs for two taxa of interest. We compared the uncorrected K_S rates and subsequent WGD hypotheses with rates corrected using ksrates ver. 1.1.1 (Sensalari et al., 2022). Using a set of transcriptomes, ksrates accounts for differences in synonymous divergence attributable to lineage-specific rate heterogeneity using two focal taxa with an inferred WGD based on an initial K_S analysis and three to four outgroup taxa to correct K_S values. We selected taxa that represent major lineages and used BUSCO scores to inform our selection process.

Sets of transcriptomes were carefully selected based on hypotheses about the relative placement of WGDs from K_S plots and previous studies (Huang et al., 2019; One Thousand Plant Transcriptomes Initiative, 2019). Taxa were picked based on their phylogenetic position to represent major clades and we used transcriptome completeness (i.e., BUSCO scores) to aid in our decisions of which taxa to include. For taxa with multiple samples (e.g., Psilotum nudum, Dicksonia antarctica), we selected the assembly with the highest BUSCO score for our analyses. We extracted OGs where at least 45% of these transcriptomes were present at least once. OGs were aligned and gene trees were generated as above. Using subsets of the species tree and these gene trees, we performed several MAPS analyses (Li et al., 2015; Li Z. et al., 2018) to estimate the placement of WGDs in a phylogenomic context. Average gene birth (λ) and death (μ) rates for each dataset were estimated following Tiley et al. (2016). Briefly, we generated OGs using OrthoFinder for just the taxa of interest and filtered gene families to only include those with at least one copy in the outgroup and at least one copy in any other taxon in the dataset and removed any gene families with greater than 100 members in any taxon. We generated ten random subsets of 500 gene families to estimate λ and μ with WGDgc ver. 1.3 (Rabier et al., 2014) using the geometric mean of gene family size as the root prior and the “oneInBothClades” likelihood conditioning. To create a null distribution, we generated three sets of 1,000 simulated trees without a WGD using GuestTreeGen from PhyloGenData (Sjöstrand et al., 2013): 1,000 with the maximum likelihood values of λ and μ; 1,000 with λ and μ at half the ML values; and 1,000 with λ and μ at three times the ML values following (Li F. W. et al., 2018, Li Z. et al., 2018). We then generated 3,000 trees as above with WGDs in the midpoint of branches of interest (e.g., leading to nodes with higher subtree duplications relative to other nodes, nodes corresponding to K_S peaks) to create a positive distribution with a retention rate of 20%. We compared the MAPS analysis with the empirical data against these simulations using Fisher’s exact tests.

We further compiled haploid (gametophytic) chromosome counts (n) and 1C genome sizes from Fujiwara et al. (2021), the Chromosome Count Database (Rice et al., 2015), and the Kew Plant C-Values Database (Leitch et al., 2019). For more recent inferred WGDs and where there was a direct sister clade without an inferred WGD, we compared chromosome number and genome size of the inferred polyploid taxa to their sister clade.

Transcriptomes were blasted against the Araport11 A. thaliana protein dataset (Berardini et al., 2004). Gene ontology (GO) terms were sorted and tallied into GO Slim categories and visualized using custom perl scripts. Paralogs falling within ±1 SD of the mean of the putative WGD peak(s) and whose posterior probability was highest for the mixture component corresponding to the inferred WGD were also annotated using this pipeline. We used Chi-squared tests to determine if the GO composition of paralogs were significantly different from their respective full transcriptomes. Following Barker et al. (2008) and Shi et al. (2010), GO Slim categories with residuals >2 are considered to be overrepresented in paralog sequences compared to the full transcriptomes, while GO Slim categories with residuals <-2 are considered to be underrepresented. We further compiled GO annotations from all transcriptomes and paralogs for each WGD event identified along the phylogeny and compared the GO Slim composition of paralogs retained in each of these events.

We compiled a total of 261 fern transcriptomes that represent 230 species; after filtering for transcriptome completeness, there were a total of 242 transcriptomes representing 211 species from 43 families (89.6% of 48 total families) in all 11 orders. We found that BUSCO completeness was significantly different between all three studies (Supplementary Figure 1A, χ² = 74.565, df = 2, P < 0.001), as well as the total number of scaffolds and contigs (Supplementary Figure 1B, χ² = 154.57, df = 2, P < 0.001). There were additional significant differences in total transcriptome length (Supplementary Figure 1C, χ² = 75.388, df = 2, P < 0.001) and number of predicted genes (Supplementary Figure 1D, χ² = 9.5959, df = 2, P < 0.01) in which the 1KP assemblies were significantly shorter and had fewer predicted genes compared to the assemblies from the Qi et al. (2018) and Shen et al. (2018) studies, but the other assemblies were not different from one another. Assembly statistics for transcriptomes from each of the studies are given in Table 1 and Supplementary Appendix 1. Note that these analyses and statistics do not include the assemblies for C. sorbifolium and S. campylurum since we only had access to the CDS files.

OrthoFinder identified a total of 6,507,715 genes from 244 transcriptomes (239 ferns and 5 outgroups). Nearly all genes (94.9%) were assigned to one of 126,747 shared OGs. The mean and median OG size was 48.7 and 3.0 genes, respectively, with 1796 OGs that were represented at least once per transcriptome. After filtering and trimming gappy sites (see Section “Materials and Methods”), 2884, 1161, and 135 single-copy OGs were retained in the SCO60, SCO75, and SCO85 datasets, respectively (Table 2). A total of 1585 multi-copy OGs were retained in the MCO dataset (Table 2).

Tree topology was generally consistent across analyses with most nodes having full or high support values (Figures 2, 3 and Supplementary Figure 2). All orders, except Gleicheniales, were monophyletic; all families were monophyletic, except for Tectariaceae and Athyriaceae in some analyses (see Section “Discussion”). In nearly all cases, genera that contained multiple samples were monophyletic, though some were not (e.g., Cheilanthes). GRF distances between trees were relatively low, with generalized scores all equal to or less than 0.06 (GRF ≤ 0.06; Supplementary Figure 3). Analyses were clustered first by the type of analysis [maximum likelihood on concatenated data matrix (ML) vs. multispecies coalescent (MSC) in ASTRAL], then generally by data type (nucleotide vs. amino acid), and finally by the dataset used (single copy, multi-copy, and percent of missing taxa). The maximum likelihood trees were the most similar to each other, having GRF ≤ 0.03, with MSC analyses differing by up to GRF = 0.06. There were clusters of low GRF values within the MSC analyses, where SCO75, SCO60, and MCO MSC analyses were clustered with GRF ≤ 0.02, with clusters for nucleotide and amino acid data types. The tree that differed the most from the other analyses was the SCO85 MSC analysis on amino acids, which had GRF values from 0.04 to 0.06 compared to other trees. Interestingly, the topology of both linked and unliked GHOST species trees were identical to the SCO60 NT MSC tree (GRF = 0, Supplementary Figures 2, 3). There were some topological differences when the first and second (CP12) and third (CP3) codon positions were analyzed separately in the SCO60 NT dataset. Tree topologies are compared in further detail for specific clades of interest in the Discussion.

Our penalized likelihood dating analysis resulted in family stem ages similar to those in previous studies (see Supplementary Appendix 3 and Supplementary Figures 4–6). We estimated that ferns originated around 346.7 Ma (range 342.5–346.7 Ma), leptosporangiate ferns around 300.9 Ma (range 299.3–312 Ma), Polypodiales around 185.0 Ma (range 172.2–222.0 Ma), and eupolypods around 114.0 Ma (range 79.2–123.6 Ma). With the three additional taxa and 17 fossil calibrations (see Section “Materials and Methods”) we estimated that ferns originated around 346.5 Ma (range 345.7–346.7 Ma), leptosporangiate ferns around 299.1 Ma (range 299–300.3 Ma), Polypodiales 166.9 Ma (range 160.6–181.9 Ma), and eupolypods around 95.3 Ma (range 86.2–111.6 Ma). Divergence times of specific families and clades are further compared in the Section “Discussion.”

Of the 239 fern transcriptomes analyzed (not including the three additional taxa in the supplemental analyses), 193 had evidence of at least one peak in their K_S plots (80.75%); of these, 163 transcriptomes had one peak, and 30 had two peaks (Supplementary Appendix 4 and Supplementary Figure 7). Median peak K_S values ranged from 0.104 to 2.191, with nearly all inferred WGDs (25 out of 27) having median peak K_S < 2, suggesting that these analyses do not suffer from saturation effects. Briefly, at high synonymous substitution values (typically K_S > 2), a build-up of synonymous mutations not related to a duplication event may appear as a “saturation peak” in K_S plots (Vanneste et al., 2013), resulting in a false-positive inference of a WGD. The two events with median K_S > 2 were supported using MAPS (see below). In general, SiZer identified significant increases in slope consistent with K_S plots (Supplementary Figure 7), although there were several cases where significant increases were detected when no discernable K_S peaks were found (Supplementary Figure 7). Furthermore, SiZer failed to detect any significant increases in slope of K_S distributions of transcriptomes that had more than one peak, such as Alsophila spp. (Supplementary Figure 7). We conducted a total of 18 MAPS analyses that utilized 147,273 empirical gene trees with an average of 8663 gene trees per analysis (range 8043–11334), and 102,000 simulated gene trees (Supplementary Appendix 5). While WGD inferences from corrected and uncorrected K_S values were generally consistent, there were some instances where they conflicted (Figure 4). Using a combination of evidences, we inferred 27 large- and small-scale WGDs throughout the phylogeny (Figure 2, Supplementary Figures 8, 9, Supplementary Appendix 5); 18 large-scale WGDs (present in more than one sampled taxon) and 9 small-scale WGDs (present in only one sampled taxon). Sixteen of these WGDs were supported in both the MAPS and K_S analyses.

There were several independent, large-scale inferred WGD events in the extant eusporangiate ferns, with WGDs separately shared by Equisetales (EQUI), Sceptridium (OPHIO.1), Ophioglossum (OPHIO.2), Psilotales (PSIL.1), Tmesipteris (PSIL.2), and Marattiales (MARA) (Figure 2). A small-scale, unshared duplication in Ophioderma pendula (OPHIO.3) was also identified. In these lineages it was difficult to assess an additional lines of chromosome number or genome size evidence as there were no direct sister groups lacking an inferred WGD. For example, all Equisetum species sampled have n = 108 with a base chromosome number of x = 108, suggesting that they are all diploid (or ancient polyploids). Furthermore, repeated rounds of neopolyploidy or multiple cytotypes within a genus or species may obscure signals of paleopolyploidy events. For example, Sceptridium, Ophioglossum, and Ophioderma should have the same haploid chromosome number if they each underwent independent WGDs. However, the haploid chromosome number in Sceptridium dissectum is n = 45, while Ophioglossum vulgatum has n = 240–1320 and O. pendula has n = 370. In these instances, we could not use chromosome number as a line of evidence to support or reject our hypothesized WGD events. We do find chromosomal support for the PSIL.2 event, with the most common haploid chromosome number of P. nudum of n = 52 (although several counts have been reported ranging from n = 46–210), compared to n = 104 for Tmesipteris tannensis. Furthermore, the genome size estimated for P. nudum (1C = 32.8 pg) is just about half that reported for T. tannensis (1C = 74.84 pg), supporting our inference of the PSIL.2 WGD event.

While there was not a significant increase in duplicated gene trees relative to a null distribution in the putative WGD in Marattiales, there was a significant increase at the MRCA of Angiopteris fokiensis and Ptisana pellucida which is consistent with a WGD (Supplementary Appendix 5). However, both corrected and uncorrected K_S plots suggests a WGD shared by all Marattiales (Figure 4). Furthermore, the base chromosome numbers of genera in Marattiales do not show an increase at the node identified by the MAPS analysis alone (at the MRCA of A. fokiensis and P. pellucida). There are a range of haploid chromosome counts in Marattiales; A. fokiensis has n = 40, Christensenia aesculifolia has n = 80, and the base haploid chromosome number of Ptisana is n = 39, with no discernable pattern to confirm the within-Marattiales WGD event suggested by MAPS alone. While we cannot directly compare the karyotypes of Marattiales to a sister group that has not undergone a WGD, treating Osmundaceae as a diploid outgroup with x = 22 would further support the inference of a WGD at the base of Marattiales.

In the leptosporangiate ferns, we identified an event shared by all leptosporangiates excluding Osmundales (LEPTO, but see corrected K_S plots in Figure 4 and Supplementary Figure 9 for an alternative hypothesis). In the clade comprised of Hymenophyllales + Gleicheniales, we inferred four WGDs: events shared by Gleicheniaceae (GLEI), Dipteridaceae (DIPT.1), Dipteris (DIPT.2), and Trichomanoideae (HYMN). The placement of several of these events were ambiguous in both the corrected and uncorrected K_S analyses, while MAPS found strong support for the separate events. The DIPT.2 event requires further investigation, as the haploid chromosome number of Dipteris conjugata (n = 33) is nearly a quarter that of Cheiropleuria dicuspis (n = 116). Both Dipteris and Cheiropleuria have base chromosome numbers of x = 33, so it is possible the WGD event inferred here is a burst of gene duplications. HYMN is additionally supported by chromosome counts and genome sizes, although it is not identified in the corrected K_S analysis (Figure 4 and Supplementary Appendix 5). The sampled Hymenophyllum all have n = 21 or 22 and H. badium has a 1C genome size of 16.15 pg, values which are nearly half that of those taxa hypothesized to have undergone the HYMN WGD event in Tichomanoideae, which have n = 36 and 1C genomes sizes from 25.04 to 25.24 pg (Crepidomanes minutum and Cephalomanes javanicum, respectively). If we compare the chromosome numbers of taxa that underwent the proposed DIPT.1 and GLEI WGD events to Hymenophyllum species which do not have a lineage-specific WGD event, the DIPT.1 and GLEI taxa have greater haploid chromosome numbers (n = 33–116 in Dipteridaceae and n = 39 in Gleicheniaceae) compared to Hymenophyllum (n = 22). These shifts in chromosome number may coincide with the inferred WGDs. In Schizaeales, a single unshared event was inferred in Schizaea dichotoma (SCHIZ).

In Salvinales, we found evidence for the WGD identified by Li F. W. et al. (2018) and One Thousand Plant Transcriptomes Initiative (2019), shared by Azolla species (SALV). There is greater than double the number of chromosomes in a haploid Azolla pinnata genome (n = 22) compared to Salvinia natans (n = 9). Interestingly, S. natans has a much greater genome size (1C = 1.82 pg) compared to Azolla filiculoides (syn. A. cf caroliniana), though this is not the case in other Salvinia species (e.g., S. culcutta has 1C = 0.25 pg). Three major WGD events were inferred in Cyatheales, one at the base shared by the order (CYATH.1), one shared by Culcitaceae and Plagiogyriaceae (CYATH.2), and one shared by Cibotiaceae, Cyatheaceae, and Dicksoniaceae (CYATH.3). In the MAPS analysis of CYATH.2, the empirical proportion of gene trees duplicated at the corresponding node of interest was significantly greater than the null distribution, but less than the positive distribution (Supplementary Appendix 5). While there are clear peaks in uncorrected K_S plots and evidence from MAPS, chromosome numbers are relatively uniform in Cyatheales. For example, the taxa that are proposed to have undergone the CYATH.3 WGD event have haploid chromosome numbers n = 65–69, while Thrysopteris elegans has n = 78. Furthermore, genome sizes of Cibotium barometz and Alsophila spinulosa (1C = 4.48 and 7.36 pg, respectively) are lower than that of T. elegans (1C = 10.23 pg). A similar pattern is observed in the CYATH.2, though there is greater variation in chromosome number n = 68–130 in Plagiogyria and Culcita. When we evaluated corrected K_S values, we found that the most recent peaks thought to represent CYATH.2 and CYATH.3 were found to be older than the inferred position by MAPS. The corrected K_S analysis instead suggests a single duplication at the base of Cyatheales, with an older event shared by Polypodiales + Salviniales + Cyatheales. It is possible that these inferred events within Cyatheales are bursts of gene duplication; genome sequencing projects focused on the tree ferns are in progress (or have been recently published, Huang et al., 2022) and should aim to tackle this question. We inferred two independent events in Lindsaeineae, one event in Lonchitis hirsuta (LONCH) and one shared by Lindsaeaceae (LINDS).

In Polypodiales, there was one inferred WGD shared by the order excluding Lindsaeineae and Saccolomatineae (POLY). One shared event was identified in the vittaroids within Pteridaceae (PTER.3), along with two independent events in Adiantum raddianum (PTER.2) and Ceratopteris thalictroides (PTER.1). Relative to Adiantum, the vittaroid ferns had much higher haploid chromosome numbers (n = 120) and genome sizes (1C = 32.81 and 33.17 pg for Antrophym callifolium and Vittaria lineata, respectively), which are nearly quadruple the chromosome numbers (n = 29, 30) and six times greater the genome sizes (1C = 5.17–5.58 pg for A. caudatum and A. aleuticum) of diploid Adiantum. Within the eupolypods, several small-scale events were also identified in Asplenium loriceum (ASPL.2), Asplenium polyodon (ASPL.1), Stenochlaena palustris (BLECH), and Pyrrosia subfuracea (PYRO). Compared to the base chromosome numbers in their respective genera (Asplenium x = 36, Stenochlaena and Pyrrosia x = 37), these taxa have twice the number of haploid chromosomes except for Pyrrosia subfurfuracea (n = 37).

By comparing the number of putative paralogs within ±1SD of the K_S peak mean to the number of predicted genes from TransDecoder, we estimated that gene retention in ferns is low, with 11.97% of genes remaining in duplicate following an inferred WGD event (range 4.10–20.35%). Over- and under-represented GO Slim categories were similar throughout taxa and events (see Section “Discussion,” Supplementary Appendix 6, Supplementary Figures 11, 12) although lineage-specific differences are present. Generally, processes involved in binding (protein binding, nucleic acid binding, DNA binding, RNA binding), responses to stimuli (response to endogenous stimulus, response to abiotic stimulus, response to chemical), and certain organelles (nucleus, ribosome, endoplasmic reticulum, Golgi apparatus) were over-represented in paralogs. Processes involved in transport (transporter activity, transport, nuclear envelope), organelles derived from endosymbiotic events (mitochondrion, chloroplast), and signaling (signaling receptor binding and activity) were under-represented in paralogs.

Our understanding of the fern phylogeny has improved substantially over the last several decades, as our field has transitioned from morphological to molecular to genomic-based phylogenetic methods. The phylogenies we reconstructed here are largely consistent with the topologies from most recent plastid and nuclear analyses of ferns. Below we discuss several areas of the fern phylogeny that have been and remain difficult to resolve; for each clade we discuss relationships, gene tree-species tree conflict, and inferred ages, focusing primarily on the SCO60 NT MSC tree/dataset which contains the largest number of loci, and highlight points of discordance between this, our other datasets, and the literature.

Despite recent advances in fern systematics and our ability to use thousands of markers throughout the genome to reconstruct phylogenetic relationships, challenges certainly remain in fully resolving the fern tree of life. In particular, we address the difficulties posed by lineage-specific rate heterogeneity, nuclear-plastid incongruence, and polyploidy.

Polyploidy has long been recognized as an important evolutionary mechanism in plants (Manton, 1950; Stebbins, 1950; Klekowski and Baker, 1966), although our view of the role of polyploidy in plant evolution has shifted dramatically throughout the past century (reviewed by Soltis et al., 2014b). One of the major debates in the polyploid literature is the role of genome duplications in diversification. Polyploid plants have been suggested to diversify at faster rates than diploids (Soltis et al., 2014a; Tank et al., 2015; Landis et al., 2018; Román-Palacios et al., 2020; but see Mayrose et al., 2011; Arrigo and Barker, 2012; Mayrose et al., 2015 for contrasting results). Although we do not address diversification here, it is of note that the two inferred WGDs shared by the largest number of taxa are at the base of the most diverse fern clades. By placing WGDs on the phylogeny constructed by Qi et al. (2018); Huang et al. (2019) found that clades of ferns that underwent two or more WGD events had higher diversification rates than other clades. They propose that the three rounds of fern radiations (Rothwell, 1987; Rothwell and Stockey, 2008; Schuettpelz and Pryer, 2009) correspond to WGDs at the base of the leptosporangiates, core leptosporangiates, and Polypodiales (Huang et al., 2019). While this is an enticing possibility, further testing will be required to explore these patterns, such as running these analyses with the varying hypotheses of the placement of WGDs in ferns posited by Huang et al. (2019), One Thousand Plant Transcriptomes Initiative (2019), and this study. Furthermore, generating a phylogeny with more comprehensive sampling or placing WGDs on an existing tree with extensive sampling may be more appropriate for these diversification analyses. Additional work will be required to determine if there is a relationship between polyploidy and diversification in ferns.

Paralog age distributions (K_S plots) have been the primary tool of choice for inferring WGDs; briefly, the synonymous distances between genes within a gene family (often defined by a clustering algorithm such as OrthoMCL, Li et al., 2003) are calculated and a histogram of the frequency or count of these distances is plotted. Peaks in the K_S distribution are interpreted as the result of large-scale duplication events (i.e., polyploidy; Lynch and Conery, 2000; Vanneste et al., 2013). By rigorously testing the applicability of K_S plots to infer ancient WGDs (paleopolyploidy events), Tiley et al. (2018) found that K_S plots should be used primarily as a hypothesis-building tool and should be supplemented with other lines of evidence. One of the biggest challenges in using paralog K_S values relative to speciation events (ortholog divergence) to infer WGDs is lineage-specific rate heterogeneity. Several methods have been developed to account for these rate differences (e.g., Barker et al., 2008; Sensalari et al., 2022) to yield corrected K_S values, which may reveal different phylogenetic placements of WGDs. Interestingly, when we compared WGD inferences from corrected and uncorrected K_S most of the WGD inferences had identical phylogenetic placements (12 identical phylogenetic placements out of 20, including LONCH; Figure 4). As expected, correcting K_S in lineages with shifts in molecular evolutionary rates such as the tree ferns and Osmundaceae resulted in different placements of the WGD; however, phylogenomic approaches generally supported the uncorrected K_S placements, rather than the corrected K_S placement. Even after correcting for rate heterogeneity, there were some cases where K_S between one-to-one orthologs were ambiguous (Figure 4 and Supplementary Figures 8, 9), likely due to rapid divergences (e.g., GLEI, DIPT.1) or differences in rates of molecular evolution (e.g., CYATH.1-3).

Recent advances in phylogenomic methods have been applied to placing WGDs on a species tree (e.g., Li et al., 2015; McKain et al., 2016; Li Z. et al., 2018). These methods examine the proportion of gene trees with a shared duplication event at the nodes of a species tree; many shared duplications support a large-scale duplication event. One of the limitations to MAPS is that the software requires a ladderized tree as input, which requires users to subsample their phylogeny. To test whether there was an effect of the sample choice, we ran MAPS using each P. nudum sample in the PSIL analysis (see Supplementary Appendix 5), which resulted in nearly identical proportions of gene duplicated gene trees at each node in the species tree (e.g., <1% difference). Similarly, using different lineages may impact the inference of WGDs. We tested whether this may impact our results by using different Equisetum species in the EQUI analysis (see Supplementary Appendix 5) and again found minimal changes to the proportion of duplicated gene trees at each node and no change in the inference of the WGD event. While we did not find that taxon selection had an impact on the inference of the phylogenetic placement of WGDs, further work should be considered to determine if and how this process can affect inferences. Both K_S plots and gene tree-based methods are accessible for data generated from transcriptomes. Full genome assemblies may further be interrogated to assess duplication events using synteny (i.e., gene order; Li F. W. et al., 2018; Krabbenhoft et al., 2021). To infer WGDs in this study, we used a combination of synonymous distance between paralogs within species and one-to-one orthologs between species, a phylogenomic method implemented in MAPS and rigorous statistical testing to place WGDs on the fern phylogeny. In general, these lines of evidence were consistent.

Overall, the placement of WGDs identified here was similar to previous estimates (Huang et al., 2019; One Thousand Plant Transcriptomes Initiative, 2019). Of the 18 large-scale WGDs we inferred, 12 had identical placements in the phylogeny produced by Huang et al. (2019). Below we discuss events which are novel or whose placement differed between our study and previous inferences, starting at the base of phylogeny. We emphasize that the inferred events should be treated as hypotheses and require further study, including whole genome analyses.

In contrast with One Thousand Plant Transcriptomes Initiative (2019), we did not find a duplication shared by all ferns, although the placement of our duplications suggest that nearly all ferns sampled have a polyploid history (Figure 2). The placement of the duplication at the base of Ophioglossales found in One Thousand Plant Transcriptomes Initiative (2019) was not recovered in our study, although a WGD shared by Ophioglossum was identified (OPHIO.2). Unlike Huang et al. (2019), we did not find that this duplication was shared by Ophioglossum and Ophioderma; rather we found that Ophioderma underwent a separate duplication event (OPHIO.3). These inferences are supported by several lines of evidence including corrected and uncorrected K_S values and MAPS. In Psilotales, we identified a novel duplication in Tmesipteris relative to Psilotum, with support from both K_S analyses, MAPS, and karyotypes. All three studies agree on the placement of a WGD at the base of the leptosporangiates excluding Osmundales (LEPTO). While Huang et al. (2019) and One Thousand Plant Transcriptomes Initiative (2019) found that Osmundales underwent a separate duplication event our analyses suggested that the uncorrected K_S inference was ambiguous, the corrected K_S suggested an event shared by all leptosporangiate ferns, and MAPS failed to identify a significant difference in the observed and null distribution of duplicated gene trees. While One Thousand Plant Transcriptomes Initiative (2019) did not identify duplication events in Hymenophyllales and Gleicheniales (only two taxa were included from these groups in their phylogeny), our findings are similar to Huang et al. (2019), although we found novel WGDs on the branches leading to Dipteris (DIPT2) and Trichomanoideae (HYMN), both of which had taxa that were sampled in Huang et al. (2019).

Unlike One Thousand Plant Transcriptomes Initiative (2019), we did not infer a WGD at the base of Schizaeales, and the inferred event may instead be representative of LEPTO. The uncorrected K_S plots for Azolla spp. were ambiguous but the corrected K_S clearly identified the inferred WGD shared by Azolla based on MAPS; syntenic analysis further support this placement (Li F. W. et al., 2018). Although both Huang et al. (2019) and One Thousand Plant Transcriptomes Initiative (2019) identified events on the backbone of the phylogeny either shared by Cyatheales + Polypodiales or Salviniales, Cytheales, and Polypodiales, respectively, we did not identify a WGD event occurring at either of these locations. Within Cyatheales, Huang et al. (2019), One Thousand Plant Transcriptomes Initiative (2019), and Huang et al. (2022) identified a WGD shared by Alsophila, although our analyses suggests instead it is shared by Cibotiaceae, Cyatheaceae, and Dicksoniaceae (CYATH.3). Although Huang et al. (2022) generated a chromosome-scale genome for A. spinulosa, they used the same approaches used here (K_S analyses and MAPS), although further syntenic evidence will be required to verify the placement of CYATH.3, especially given the uniformly high chromosome numbers in Cyatheales (see Section “Results”) and slower rates of molecular evolution in the tree ferns. Along the backbone of the phylogeny, we inferred a WGD shared by Polypodiales excluding Lindsaeineae and Saccolmatineae (POLY), the placement of which agrees with One Thousand Plant Transcriptomes Initiative (2019). Within Pteridaceae, we found one shared event at the base of the vittaroid ferns (PTER.3) that was not identified in One Thousand Plant Transcriptomes Initiative (2019). We did not find evidence of the WGD on the branch leading the eupolypods, which was recovered in Huang et al. (2019). Although we did identify a significantly greater proportion of duplicated gene trees than in the null distribution at the MRCA of the eupolypods, the observed distribution was not consistent with positive simulations of a WGD.

As with most other plants, we found that the rate of duplicate gene retention is low following paleopolyploidy or large-scale gene duplications, with an average of around 11.97% of genes remaining in duplicate. On average, duplicate genes have a retention rate around 10% in plants (Tiley et al., 2016), but vary depending on the age of the duplication (Li et al., 2016), with up to 76.3% of genes duplicated in Glycine max (Tiley et al., 2016) which underwent a recent duplication event. Genomes which have undergone more recent duplications tend to show higher retention rates than older duplications, with the proportion of the gene families duplicated rapidly decaying with the age of events (Li et al., 2016). Given that the gene families we analyzed were not restricted to “core gene families” as in Li et al. (2016) and that the signal of WGDs becomes diluted as synonymous divergence increases, we were not able to replicate these findings in ferns. While our method of identifying paralogs from WGDs using K_S plots has been described and used before (e.g., Li Z. et al., 2018), it is important to note that other types of duplication events other than WGDs may contribute paralogs to the distribution. In some cases, particularly in older WGDs, peaks may be difficult to distinguish from the background processes affecting duplicated genes. Although we assigned paralogs to WGD events based on their probability of belonging to certain components in mixture models (see Section “Materials and Methods”), we may be including other duplicated genes that may not have arisen from WGDs.

While most duplicated genes are rapidly lost (Lynch and Conery, 2000), the functions and types of genes retained in duplicate following independent WGD events are similar (e.g., Barker et al., 2008; Freeling, 2009; Li et al., 2016; Li Z. et al., 2018; Figure 5). According to the DBH (Papp et al., 2003), the loss of duplicated copies of some but not all partners in an interaction would alter the stoichiometry of the product/network and therefore the ultimate function may be disrupted. Compared to full transcriptomes, we found that genes retained in duplicate have functions that are over-represented in binding processes (DNA binding, transcription factor activity, nucleic acid binding, protein binding, RNA binding, nucleotide binding), responses to stimuli (response to chemical, response to endogenous stimulus, response to abiotic stimulus), and certain organelles involved in gene and protein production and processing (nucleus, ribosome, endoplasmic reticulum, Golgi apparatus) (Figure 5). Binding processes (DNA binding, RNA binding, nucleic acid binding) were also found to be over-represented in paralogs retained in hexapods (Li Z. et al., 2018) and intermediate-age duplicates in angiosperms (Li et al., 2016). Processes involved in transport (transporter activity, transport, nuclear envelope), organelles derived from endosymbiotic events (mitochondrion, chloroplast), and signaling (signaling receptor binding and activity) were under-represented in paralogs (Figure 5). Similar patterns of gene loss are again seen in hexapods (transport, mitochondrion, nuclear envelope, Li Z. et al., 2018), and angiosperms (transport, transporter activity, mitochondrion, Li et al., 2016). While lineage-specific variation in patterns of gene retention is present (see below), some patterns appear to be conserved over deep evolutionary time among the kingdoms of life.

While the overall pattern of gene functions is similar across independent duplication events, there were several instances that suggest that other factors may also drive lineage-specific differences in retained duplicates. For example, retained duplicates following the inferred Psilotales (PSIL.1) and Tmesipteris (PSIL.2) events were over-representative of genes involved with endosymbiotic organelles (chloroplast, mitochondrion) and signal receptor binding, while duplicates from most other events were under-represented in those categories (Figure 5). Similarly, genes with functions related to the ribosome and cytosol were under-represented in retained duplicates following the Dipteris event (DIPT.2) (Figure 5), although this may be due to low sample size. In Asteraceae, Barker et al. (2008) found several GO Slim categories that were under-represented in duplicates following independent WGDs which we found to be over-represented in fern duplicates (e.g., DNA or RNA binding, nucleus) whereas others were over-represented (e.g., cytosol, protein metabolic process) or under-represented (e.g., chloroplast) in both analyses. Furthermore, the age of the inferred WGD may impact the function of gene sets retained, as in angiosperms (Li et al., 2016); for instance, genes with functions related to translation and metabolic processes were over-represented in older WGDs, but under-represented in more recent events and single-copy gene families. Gene function may therefore be an important factor in the long-term survival of duplicated genes.

While genic diploidization (i.e., the process of removal and loss of genes by molecular mechanisms, Li et al., 2021) is clearly occurring in ferns, perhaps through pseudogenization/gene silencing as hypothesized by Haufler and Soltis (1986), cytological diploidization may be slow to follow. Unlike other plants, chromosome number and genome size are positively correlated in ferns (Clark et al., 2016) and, taken with the relative stasis in genome size across ferns, may suggest that chromosomes are retained following WGDs rather than lost during genomic reorganization in angiosperms. A similar pattern has been observed in the catostomid fish Myxocyprinus asiaticus, which shows remarkable genome subgenome stability and retained synteny over 50 MY following a WGD (Krabbenhoft et al., 2021). In contrast, one subgenome tends to dominate in polyploid plants (Alger and Edger, 2020) and rearrangements drastically alter gene order and retention of synteny (Zhao and Schranz, 2018; but see VanBuren et al., 2020 for an example of subgenome stability in plants following a WGD 1 Ma). Whether similar patterns of genome reorganization or stability are present in ferns remains uncertain and an active area of research.

Importantly, transcriptomes are temporal and spatial “snapshots” of gene expression. Many of the transcriptomes used in this study were derived from young leaf material, although some are from other tissue types (e.g., fertile pinnae, gametophytes). Not all genes in the genome will be expressed in every tissue and therefore transcriptomes from single tissues will likely not represent the entirety of gene-space in a genome. The presence or absence of genes in a transcriptome may not necessarily correlate to the presence or absence of that gene in the genome but could rather be a product of differences in expression between tissues and time. The analysis presented here is one of the first to tackle functional gene retention following WGDs in ferns, but additional analyses will be needed to explore whether these patterns of gene retention are specific to transcriptomic study. As new genome assemblies are becoming available (e.g., Adiantum, Alsophila, Ceratopteris) patterns of gene retention should be further explored in more complete gene spaces.