A total of 47 Schizaeaceae specimens from 27 species (~ 77% of the species diversity of the family; PPG I, 2016) were sampled, including seven Microschizaea spp., seven Actinostachys spp., and 13 Schizaea spp. (Supplementary Table 1). This sampling covered all bioregions in each genus/clade (Table 1). For outgroups, we sampled the other two Schizaeales families: Anemia phylliditis from Anemiaceae and Lygodium japonicum from Lygodiaceae (PPG I, 2016). A modified CTAB protocol was used for DNA extractions (Kuo, 2015). The rpoC2, rbcL, and trnL-L-F (trnL gene + trnL-F intergenic spacer) were sequenced for our three-plastid-region dataset. PCRs were performed in 15 μl reactions each with 20 ng of genomic DNA, 0.5 μM of each primer, and 1 × SuperRed PCR Master Mix RED (TOOLS, New Taipei City, Taiwan). The resulting PCR products were purified and sequenced using the standard Sanger method with Applied Biosystems 3730XL (Thermo Fisher Scientific, Waltham, MA, United States of America) at Genomics Corp. (New Taipei City, Taiwan).

For the plastome phylogeny, we selected both Schizaeales outgroups and a total of 11 Schizaeaceae representatives. These included species from every genus and major clade within Schizaeaceae (two Microschizaea spp., two Actinostachys spp., and seven Schizaea spp.; Supplementary Table 1). Among these, published plastome sequences were already available for L. japonicum and three of the 11 species of Schizaeaceae (Gao et al., 2013; Labiak and Karol, 2017), and these were used for all downstream analyses. To assemble the remaining plastomes, we used the Illumina reads generated by either the GoFlag project (Breinholt et al., 2021) or the genome skimming of this study. When the sequencing depth of a GoFlag sample was insufficient to yield a circular plastome using NOVOplasty (see below for details), we designed PCR primers and closed the gaps between the contigs. The PCR recipe was the same as described earlier. For the genome skimming, we first sheared the genomic DNA into 400 ~ 500 bp fragments using a Covaris S2 ultrasonicator (Covaris, Woburn, MA, United States of America), which were then input into the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, Ipswich, MA, United States of America). Sequencing was done by HiSeq X Ten (Illumina, San Diego, CA, United States of America) with 150 bp PE and ~ 3Gb per sample. Fastp (Chen et al., 2018) was used to trim the reads using the default settings. NOVOplasty (Dierckxsens et al., 2017) was used to assemble the plastomes with the setting of “Kmer = 39,” and conspecific rbcL sequences were used as the input seeds. These plastome assemblies were annotated using Geneious (Kearse et al., 2012) with the published Schizaeaceae plastomes (GenBank accessions: KU764518, KX258660-61) as references. We manually inspected every protein-coding gene annotation, and adjusted the coordinates if necessary.

Details about the PCR primers are provided in Supplementary Table 2. Details about the voucher information, GenBank accessions, and NCBI Sequence Read Archive accessions can be found in Supplementary Table 1.

For the three-plastid-region dataset (i.e., “rpoC2 + rbcL + trnL-L-F”), the nucleotide sequences were first aligned with MUSCLE (Edgar, 2004), as implemented in AliView (Larsson, 2014), and then concatenated into a single matrix. This three-plastid-region matrix was partitioned by gene, by intergenic spacer (IGS), and by codon position in order to find the best partition scheme and substitution models using ModelFinder (Kalyaanamoorthy et al., 2017) with the Bayesian information criterion (BIC; Schwarz, 1978). Based on the inferred partition scheme and the substitution models, IQtree 1.6.8 (Nguyen et al., 2015) was used to construct the maximum likelihood (ML) phylogeny with 1,000 ultrafast bootstrap replicates (UFBS). The Bayesian phylogeny was inferred using MrBayes 3.2.7 (Ronquist et al., 2012). Two simultaneous runs were carried out with four chains (10 million generations each). Each chain was sampled every 1,000 generations. Log likelihoods of MCMC runs were inspected in Tracer 1.6 (Rambaut and Drummond, 2013) and RWTY (Warren et al., 2017) to confirm their convergence. The first 25% of the sample was discarded as burn-in, and the remaining was used to calculate the maximum clade credibility consensus tree with TreeAnnotator (Rambaut and Drummond, 2013).

For the phyloplastomic datasets, we only included coding sequences, introns, and intergenic spacers (IGS) that were consistently located in the large single copy (LSC) region of the Schizaeaceae plastomes. The restriction to LSC genes aimed to reduce phylogenetic artifacts resulting from substitution rate heterotachy, because the movement of genes between the inverted repeat (IR) and single copy (SC) regions would lead to changes in their substitution rates (Li et al., 2016). For alignment of every individual DNA region, we used MAFFT v7.450 (Katoh and Standley, 2013) and MACSE v2.03 (Ranwez et al., 2011). In total, we compiled five different matrices that consist of either the coding genes, noncoding regions (i.e., “IGS/INTRON”), or both (i.e., “67CDS + IGS/INTRON”). In the three coding gene matrices, 67 loci from the LSC were included with sequences of either (1) the first two codon positions (i.e., “codon1 + 2”), (2) the third codon position (i.e., “codon3”), and (3) all three codon positions (i.e., “67CDS”). For each of the five matrices, we performed four analyses with different models and partitions, as detailed in Table 2. With all the matrix/partition/model combinations, we conducted a total of 20 (= 5 × 4; Table 2) ML phylogenetic analyses in IQtree 1.6.8 each with 1,000 UFBS replicates and a series of tree topology tests (such as KH, SH, ELW, and AU; detailed in Minh et al., 2020) with a RELL replicate number of 10,000.

To examine the potential effects of homoplasy driven by rapidly evolving sites in CDS, particularly the third codon positions, we also analyzed two degenerate-coded matrices: one “67CDS” with “Degen”-coding (Zwick et al., 2012) and one “codon3” with RY-coding of third codon positions, and compared these results with those generated from original non-degenerated ones. For these two matrices, the same region partitions and model selection were used as those inferred by the original non-degenerated matrices, and the ML phylogenetic analyses was then performed with 1,000 UFBS replicates in IQtree 1.6.8. In addition, to examine the evolution of Schizaeales chl genes in more detail, we combined their sequences with those used in Kuo et al. (2018a), which has a complete order-level sampling. These sequences were aligned with MAFFT, and used to infer an ML phylogeny in IQtree 1.6.8 with the model GTR + FO*H4 and 1,000 UFBS replicates.

We examined herbarium collections from MO, P, TAIF, UC and WELT to conduct morphological comparisons between the S. dichotoma complex from Africa/Malagasy and Asia/Oceania. The Asian/Oceanian S. dichotoma complex contains several species—S. asperula, S. bifida, S. biroi, S. dichotoma and S. forsteri—whereas the African/Malagasy complex was represented by a single species—S. medusa sp. nov. We also measured several quantitative traits from these specimens, such as rhizome thickness, frond sizes, and stipe sizes (Supplementary Table 3).

To confirm whether the nuclear-encoded ndh complex related genes was present or not in the Schizaeales, we downloaded their homologue sequences from the well-annotated case of Apostasia odorata (detailed in the Supplementary Material S2 of Lin et al., 2017), and blast-searched them against six Schizaeales transcriptomes (Qi et al., 2018; Shen et al., 2018; One Thousand Plant Transcriptomes (OTPT) Initiative, 2019) and two fern whole genomes (Li et al., 2018). For PnsB1, the homologue sequence from Marchantia polymorpha (Ueda et al., 2012) was used instead, because the annotated sequence from A. odorata (GenBank accession: KX156894) appeared to be misidentified. We also included the ndhV from Zostera marina (Ma et al., 2022), because the sequence of A. odorata appeared to be unusually diverged from other plants. The transcriptomes were obtained from previous studies which used foliar, photosynthetic tissues for RNAseq, and included two Lygodiaceae (Lygodium japonicum and Lygodium flexuosum), two Anemiaceae (Anemia tomentosa and Anemia phyllitidis), and two Schizaeaceae (Actinostachys digitata and Schizaea dichotoma). We first used the sequences from A. odorata, Z. marina, and M. polymorpha as the query for tBLASTx searches against the Azolla filiculoides and Salvinia cucullata genomes via FernBase (https://www.fernbase.org/). Based on the criterion of a percentage of identity, we obtained the best-matched sequences from the fern genomes. To confirm whether these sequences are homologues or not, we conducted tBLASTx searches against the nucleotide collections of NCBI, and inferred their gene trees based on the blast results. After obtaining the fern homologous sequences from these genomes, we then added these sequences into the tBLASTx searches against the transcriptomic assemblies using blast-2.10.0 + (Camacho et al., 2009). Only the transcripts having > 60% sequence identity and consistent hits to all query homologues (i.e., those from the two ferns, A. odorata, Z. marina, and M. polymorpha) were kept.

To examine the possibility that the missing protein-coding genes in the plastome have been transferred to and are now transcribed in the nuclear genome, we also conducted local tBLASTx searches for these genes as mentioned above. The query sequences were derived from plastomes of the closest relatives that still retain those genes (e.g., using Anemiaceae’s rps16 sequence to blast against the Schizaeaceae’s transcriptomes).

In total, nine plastomes were assembled in this study, including one Anemia and eight Schizaeaceae. They all assembled into a single circular contig, except for the plastome of Microschizaea tenella in which the gap between psbK and trnQ could not be closed by PCR and was thus separated into two contigs. All members of Schizaeaceae share several plastome features, including the expansions of the IR and losses of all ndh genes (ndhA–J), rps16, ycf94, and several tRNA genes (Figure 1). IR/SC boundaries of the Schizaeaceae plastomes are almost identical, and their small single copy region (SSC) are highly reduced (Figure 1). In the most extreme case, the SSC of S. pusilla is only 1,091 bp in length and contains only one gene, trnN. Actinostachys, on the other hand, appears to be an outlier where the SSC expanded to include ccsA, rps15, and ycf1. The chl genes (chlB, chlN, and chlL) are absent in the plastomes of two clades—S. pusilla and Actinostachys, but this is unlikely to have resulted from a single loss event (Figure 2; see below in “Discussion”). The loss of psaM gene was shared by Schizaeaceae and its sister family—Anemiaceae (Figure 2).

From our three-plastid-region and plastome datasets (Figures 2, 3) four well-supported clades were recovered (UFBS = 100, Bayesian inference = 1.00) within Schizaeaceae—Schizaea, S. pusilla, Actinostachys, and Microschizaea (Figure 3). These clades closely matched their generic definition by Reed (1947), except for S. pusilla which formed a rather distinct clade from other Microschizaea spp. However, the inferred inter-clade relationships differed among our phylogenetic analyses (Table 2). Generally, topology I received the highest branch supports, and another two were rejected by most of our topological tests (Table 2). Only the “codon3” matrix with a gene-partitioned GTR + F + R10 model highly supported topology II (Table 2), but this support was declined or nearly unchanged when applying the degenerate-coded matrices with analyses (Supplementary Table 4).

In the Schizaea clade, the S. dichotoma complex was non-monophyletic with two separate lineages, one from Africa/Madagascar, and the second from Asia/Oceania (Figure 3). In the Microschizaea clade, M. robusta and M. australis were nested within M. fistulosa, rendering M. fistulosa polyphyletic (Figure 3).

We detected the nuclear-encoded ndh-related genes in the transcriptomes of other Schizaeales (Anemiaceae and Lygodiaceae) and the genomes of Salvinia and Azolla, except for pnsL1–4, psnB2, ndhV, CRR2–4, and CRR21 (Figure 4). For CRR41, we identified it in Salvinia but not in Azolla. This gene likely existed in other Schizaeales (Anemiaceae and Lygodiaceae), whose CRR41-matched transcripts showed consistent blast-hits to the homologous queries of CRR41 although relatively low in the identity ranging from 40–54%. In contrast, all plastid and most nuclear-encoded ndh-related genes were missing in the Schizaeaceae transcriptomes (Figure 4). We were also unable to detect any chl genes in the transcriptome of Actinostachys, whose plastomes also lack these genes. The plastid rps16 and ycf94 genes were also absent from the transcriptomes of Schizaeaceae, and psaM was absent from both Schizaeaceae and Anemiaceae.

Our phylogenetic analyses based on different datasets, models, and partitions resulted in three general topologies (Figure 2). While data from plastome structure and gene content failed to provide additional evidence, we have more confidence in topology I (Table 2). This topology received, on average, higher branch supports (i.e., nodes A and B), and was not rejected in most of the topology tests (Table 2). On the other hand, although topology II also received high branch supports in some analyses, it was mostly derived from matrices comprised of the third codon positions. It is likely that the support for topology II is a result of high homoplasy at the third codon position (Supplementary Table 4), which is known to have higher substitution rates and therefore be more likely to reach saturation (e.g., Breinholt and Kawahara, 2013).

Notably, our study is the first to include Microschizaea species other than S. pusilla. All the inferred phylogenies suggested that Microschizaea sensu Reed (1947) is not monophyletic and fell into two distinct, non-sister clades—Microschizaea and S. pusilla (Figure 2). The Microschizaea clade comprises all Microschizaea species except S. pusilla (Figures 2, 3). The position of S. pusilla varied across the three inferred topologies, but was never found to be sister to Microschizaea. In the most highly supported topology (topology I), S. pusilla was placed sister to the Schizaea clade similar to what was reported by Wikström et al. (2002). The close relationship between the S. pusilla and the Schizaea clades is morphologically supported by the shared presence of multicellular hairs intermixing with sporangia, a trait that is absent in the Microschizaea clade (Table 1). Schizaea pusilla also differs from Microschizaea by its more ellipsoidal (i.e., bilateral) spores with a ratio of polar to equatorial length (P:E) about 0.66, compared to the subglobose or ovoid spores in the Microschizaea clade with P:E around 0.72–0.76 (Reed, 1947; Table 1). Finally, spore perines of S. pusilla are alveolate with shallow pits but those of the Microschizaea clade are granular or smooth (Table 1).

Despite bearing morphological and phylogenetic differences, Microschizaea was retained in synonymy of the genus Schizaea in the latest phylogenetic classification of PPG I (2016), and Actinostachys was there recognized as the second genus of Schizaeaceae. Here, we show Schizaea as circumscribed by PPG I is very likely paraphyletic. We therefore propose a three-genus framework for Schizaeaceae that recognizes Microschizaea, Actinostachys, and Schizaea. This classification is similar but different from Reed (1947) in that we placed S. pusilla in Schizaea rather than in Microschizaea. It can be argued that given its unresolved placement and distinct morphology, S. pusilla could be erected as a separate genus. We are however hesitant to do so based on the current data. Regarding the infra-generic system of Reed (1947), many of his subgenera and (sub)sections are non-monophyletic; further studies are needed to provide a clear picture for these infra-generic schemes.

The simple foliar structure in Schizaeaceae presents taxonomic challenges, particularly in the Microschizaea fistulosa and Schizaea dichotoma species complexes. The M. fistulosa complex contains the polyphyletic “M. fistulosa” and most of its congeneric members (Figure 3), and has a broad distribution across Southeastern Asia, Pacific Islands, Oceania, and South America (Lash, 1966; Brownsey and Perrie, 2014; Giacosa et al., 2015). Morphologically variable but continuous forms make taxonomy of this complex still unsettled. There are more than ten names within this complex but these names are always treated under a single name—M. fistulosa or with another—M. australis (Lash, 1966; Brownsey and Perrie, 2013, 2014). In addition, two cytotypes/ploidy levels have been discovered. The specimens in New Zealand can be separated into M. australis with a lower chromosome number (n = 94) and smaller overall plant sizes, and M. fistulosa with a higher chromosome number (n = c.150, 190) and larger individual sizes (Brownsey and Perrie, 2013). From our plastid tree, Pacific samples of this complex appear to form a group, while circum-Antarctic samples form another (Figure 3). However, because plastid sequences can track only the maternal lineage in ferns (reviewed in Kuo et al., 2018b), the current phylogeny is still insufficient to shed light into any polyploidization and reticulation history in the M. fistulosa complex. Future systematics studies need to incorporate cytological information, details of microcharacters, and analyses of nuclear markers.

The S. dichotoma complex likewise exhibits wide morphological and cytological variations (Brownsey and Perrie, 2014). Some forms have been formally named to reflect their distinct morphology, such as S. asperula showing sterile-fertile dimorphic fronds and S. bifida showing only twice-to-thrice bifurcated lamina. However, most members of this complex are poorly characterized, and have been collectively lumped into “S. dichotoma”, which is found to be polyphyletic in our phylogeny (Figure 3). We discovered that the African and Malagasy specimens form a clade that is phylogenetically distinct from the rest of the S. dichotoma members (Figure 3). Because the types of all the previously named species in this complex are based on Asian and Oceanian materials, a new species is warranted, which we named Schizaea medusa. By quantifying several key morphological characters, we showed that S. medusa can be distinguished from S. dichotoma (Supplementary Figure S1; also see in “Taxonomic Treatment”). There are likely additional cryptic species in the S. dichotoma complex and more studies are clearly needed.

Gametophytes of Schizaeaceae are peculiar among extant ferns because of their strongly mycotrophic lifestyle, whether they are achlorophyllous (like Schizaea and Actinostachys) or not (like S. pusilla and Microschizaea). Only one exceptional case, reporting a surface-living and chlorophyllous gametophyte in the Schizaea clade, is from the observation of S. bifida by Thomas (1902). However, this record needs to be confirmed because S. bifida usually co-occurs with Microschizaea (Brownsey and Perrie, 2014), and fern gametophytes in the field could be easily misidentified without genetic evidence (Nitta and Chambers, 2022; Wu et al., 2022). Despite being chlorophyllous, gametophytes of S. pusilla and Microschizaea still show several properties that are not found in other chlorophyllous fern gametophytes. Their gametophytes are filamentous and partially or sometimes completely underground, with their subterranean parts being achlorophyllous (Britton and Taylor, 1901; Lash, 1966; Bierhorst, 1968). In addition, previous studies failed to regenerate sporophytes from these gametophytes under axenic conditions (Lash, 1966; Swatzell et al., 1996), implying that microbial symbionts might be required for sexual reproduction (Britton and Taylor, 1901; Lash, 1966; Swatzell et al., 1996). Finally, at least in S. pusilla, the gametophytes unusually exhibit a negative phototropism during spore germination (Kiss, 1994). Taken together, Schizaeaceae are clearly heterotrophic or mixotrophic at the gametophyte stage, and have intimate associations with fungal symbionts.

Heterotrophic or mixotrophic plants often display certain diagnostic genomic signatures (e.g., Vogel et al., 2018; Su et al., 2019; Xu et al., 2021), such as plastomes with extensive gene losses (Graham et al., 2017; Hadariová et al., 2018; Wicke and Naumann, 2018). Here we confirmed that the loss of plastid ndh genes is a synapomorphy of extant Schizaeaceae (Figure 2), and we also could not detect expression of most nuclear ndh-related genes from Schizaeaceae transcriptomes. In other words, there might be a concerted loss of ndh-related genes in both genomic compartments (Figure 4), a feature that has been suggested to be the initial (but irreversible) step toward mycoheterotrophy in flowering plants (Graham et al., 2017; Lin et al., 2017). Interestingly, a separate fern lineage Stromatopteris (Gleicheniaceae), which also produces mycoheterotrophic gametophytes (Bierhorst, 1971b), has likewise experienced ndh gene losses in the plastome (Du et al., 2022).

Loss of the plastid chl genes is frequently found in plants with a heterotrophic and/or achlorophyllous nature for their gametophyte generation (e.g., all flowering plants), and is thus considered as one genomic indicator for heterotrophs (Ueda et al., 2014). These plastid chl genes encode all subunits for light-independent protochlorophyllide oxidoreductase (DPOR), which plays an important role for chlorophyll synthesis under dark, particularly during the gametophyte stage (Suzuki et al., 1998; Ueda et al., 2014). Complete lack of chl genes had been documented in two fern lineages, Psilotaceae and Actinostachys in Schizaeaceae (Grewe et al., 2013; Zhong et al., 2014; Labiak and Karol, 2017; Kuo et al., 2018a), whose gametophytes are also achlorophyllous (Bierhorst, 1971b). Here we found that in Schizaeaceae, there were actually two independent losses of chl genes (Figures 1, 2; Supplementary Figure 2). One of them, however, is found in S. pusilla, a species that produces chlorophyllous gametophytes (Figures 1, 2). The link between chl gene loss and the achlorophyllous nature of gametophytes might therefore seem not obligate. The gametophytes of S. pusilla likely rely on nuclear-encoded light-dependent NADPH-protochlorophyllide oxidoreductase (LPOR) instead, the alternative pathway for chlorophyll synthesis, but advanced genetic evidence is required to test such a hypothesis. On the other hand, several fern and lycophyte lineages producing achlorophyllous gametophytes retain functional chl genes in their plastomes, such as Lycopodiaceae, Ophioglossaceae, and the Schizaea clade in the present case (reviewed in Kuo et al., 2018a). These retentions imply that DPOR remains important for chlorophyll synthesis in the autotrophic (or mixotrophic) sporophytes of these ferns. Alternatively, DPOR could be involved in different physiological functions, resembling the situation in the (potentially) fully mycoheterotrophic bryophytes, which also have chl-retained plastomes (Bell et al., 2020).

In addition to ndh and chl genes, several other coding and tRNA genes have disappeared from Schizaeaceae plastomes (Figures 1, 2). From our blast results against the transcriptomes, the coding genes are unlikely to have been transferred to the nuclear genome. However, these genes do not seem to be associated with the specialization of gametophytic lifestyle in Schizaeaceae, and are mostly considered to have minor functions in the plastid. For instance, rps16 has been lost several times in ferns (Zhang et al., 2014; Kuo et al., 2018a; Du et al., 2022). Ycf94, which was recently identified in plastomes of seed-free plants and whose function is still unknown (Song et al., 2018), seems to be absent only in Schizaeaceae (Du et al., 2022).

Schizaeaceae is one of the most understudied families of ferns. With the most comprehensive sampling to date, our phylogenomic analyses resolved important relationships in this family, and provided the most robust infrafamiliar backbone for Schizaeaceae. We proposed a new phylogenetic classification modified from the Reed’s (1947) system, in which a portion of Schizaea sensu PPG I (2016) is moved to a third genus—Microschizaea. In addition, our species-level phylogeny illustrated species complexes within this family. One cryptic species hidden in the S. dichotoma complex was identified and described as a new species, Schizaea medusa. Using our new phylogenetic framework, we were able to trace the evolution of plastome features as well as the gametophytic lifestyle in Schizaeaceae. Specifically, we provided a better picture of the relationship between plastid gene loss and mycoheterotrophy in fern gametophytes.

Importantly, this phylogenetic study sets the stage for future work delving into the unique biological features of Schizaeaceae and in Schizaeales, such as Mesozoic biogeography (Skog, 2001), epiphytism on tree ferns (Amoroso et al., 2020), leaf simplification and lamina reformation (Vasco et al., 2013), biosynthesis of silica bodies (Ribeiro et al., 2007; Tzu-Tong Kao personal communications), and recruitment of symbiotic microbiome (e.g., Chen et al., 2022). Lastly, the pheromone-mediated (i.e., antheridiogen) mating system has been well-studied in the other two Schizaeales families (Yamane, 1991, 1998; Tanaka et al., 2014) but not yet for Schizaeaceae. It would be very interesting to explore whether gametophytes of Schizaeaceae rely on the similar pheromone system for their “underground” mating.