1. Introduction

Cryptic species, which have long puzzled biologists, are often misidentified due to morphological similarity to closely related species (Sáez and Lozano, 2005; Bickford et al., 2007; Winker, 2005; De, 2007). These species can be derived from isolation differentiation (Rieseberg and Willis, 2007), hybridization, and polyploidization (Mallet, 2007; Rieseberg and Willis, 2007). Convergent evolution (Peter et al., 2009; Bravo et al., 2015) is often considered necessary for cryptic speciation in special habitats like arctic, aquatic, desert, rocky, and alpine regions. However, over the last three decades, an increasing number of cryptic species (both animals to plants) have been identified in numerous habitats (from marine to freshwater) (Vrijenhoek et al., 1994; Hebert et al., 2004; Feulner et al., 2006; Grundt et al., 2006; Lombard et al., 2010; Pillon et al., 2010; Li et al., 2015; Liu et al., 2015; Hollingsworth et al., 2016; Hinojosa et al., 2019). More recently, DNA (deoxyribonucleic acid) barcoding has been used to identify cryptic species of lycophyte and fern genera, such as Selaginella (Wang et al., 2022), Ceratopteris (Zhang et al., 2020b; Yu et al., 2022), Hymenasplenium (Xu et al., 2019a, 2019b), and Didymochlaena (Shang et al., 2020). As biodiversity conservation based on species diversity has become a hot topic globally, accurate identification of cryptic species is critical for the conservation of biodiversity (Bickford et al., 2007; Nygren, 2014; Fier et al., 2017; Thompson et al., 2020), especially because many cryptic species are at risk of extinction prior to their discovery (Liu et al., 2022). However, the origin of species remains a fascinating scientific problem, and much of the evolution of organisms, as well as the total number of species on earth, still baffles scientists (May, 1988, 1990; Prance et al., 2000; Dirzo and Raven, 2003; Joppa et al., 2011; Pimm et al., 2014). The known number of species may be much smaller than the actual amount because significant morphological changes do not accompany all speciation, and genetic diversity within a species also implicates underappreciated mechanisms of morphologically static species (Bickford et al., 2007).

The stability and homogeneity of aquatic environments have led to morphological similarities between related aquatic organisms (Bickford et al., 2007; Fier et al., 2017). One such group of organisms is the quillworts (Isoëtes spp.), an ancient aquatic group of species. Isoëtes is the only surviving genus of the Isoëtaceae family. These heterosporous lycophytes have fascinated botanists and paleobotanists alike because of their distinctive morphology and long evolutionary history (Gifford and Foster, 1989; Pigg, 2001; Cúneo, 2009; Pereira et al., 2021a). Although the origin of Isoëtes can be dated back to the Devonian (Pigg, 1992), most extant diversity is relatively recent (Wood et al., 2020). Isoëtes are widely distributed from the tropics to subarctic regions and survive in habitats that range from submerged, evergreen aquatic environments of cold, clear water lakes to upland terrestrial ephemerals with seasonally parched soils over bedrock (Smolders et al., 2002; Troia et al., 2016).

Although more than 350 species of Isoëtes have been reported, fewer than 200 have been universally accepted by the scientific community (Troia et al., 2016). Species number varies greatly from country to country; for example, more than 50 species have been reported in the United States (Schafran, 2019), whereas only ten species have been recorded in China (Li et al., 2019; Lu et al., 2021; Shu et al., 2022). These Chinese Isoëtes vary in ploidy level, and include diploid species (e.g., I. hypsophila Hand.-Mazz.; I. shangrilaensis Xiang Li, Yuqian Huang, Xiaokang Dai & Xing; I. yunguiensis Q.F. Wang & W.C. Taylor; and I. taiwanensis De Vol), tetraploid species (e.g., I. sinensis Palmer) (Palmer, 1927; Liu et al., 2002), and hexaploid species (e.g., I. orientalis H. Liu & Q.F. Wang) (Liu et al., 2005). Interestingly, ploidy level appears to be related to geographical distribution. For example, diploids have only been reported in the Qinghai-Tibet Plateau, Southwest Sichuan, and Taiwan (Handel-Mazzetti, 1923; Devol, 1972; Wang et al., 2002; Li et al., 2019; Liu et al., 2002); in contrast, hexaploid species have only been found in Zhejiang Province. Previous studies of Isoëtes in China have ignored geographical variation in ploidy level among known Isoëtes species (Pang et al., 2003; Ye and Li, 2003; Chen et al., 2004; Taylor et al., 2004; Kang et al., 2005; Wang et al., 2006; Liu et al., 2008; Xie et al., 2019; Dai et al., 2020). Thus, nearly all the populations from Anhui, Hunan, Jiangsu, and Zhejiang have historically been treated as I. sinensis.

Species identification in Isoëtes is challenging due to morphological convergence between species that has resulted from widespread interspecific hybridization and allopolyploid speciation (Hickey, 1986; Taylor and Hickey, 1992; Hoot and Taylor, 2001). Several features have long been used to distinguish Isoëtes species, including habitat (Pfeiffer, 1922), megaspore and microspore texture (Hickey, 1986; Lellinger and Taylor, 1997; Holmes et al., 2005), and chromosome number (Takamiya et al., 1994; Liu et al., 2002). More recently, molecular barcoding methods using complete chloroplast genome data have been used to distinguish species within this genus (Choi et al., 2018; Schafran et al., 2018; Pereira et al., 2021a, 2021b). However, chloroplast barcoding methods have still failed to distinguish problematic Isoëtes species from species complexes (Romero and Real, 2015; Kim and Choi, 2016; Dai et al., 2020; Suissa et al., 2020). Nevertheless, several new Isoëtes species have been recognized and described around the world (Moraolivo et al., 2016; Pereira et al., 2016, 2017, 2019a; Schafran et al., 2016; Li et al., 2019; Lu et al., 2021; Shu et al., 2022), indicating that Isoëtes species diversity may be grossly underestimated, with many cryptic species unidentified (Schafran, 2019). One approach to identifying additional Isoëtes species is the use of larger data sets such as genomic data.

Genomic-scale data sets have been used to resolve complex evolutionary relationships, including the branching order within rapid radiations (Wei et al., 2017; Wei and Zhang, 2019). Plastid genome (plastome) DNA sequences have been extensively used in recent plant molecular systematics because of their mode of uniparental inheritance and high content of informative loci (Givnish et al., 2010, 2015; Jansen et al., 2011; Ruhfel et al., 2014; Lu et al., 2015; Ross et al., 2016; Schafran et al., 2018). Over the past few years, several studies have demonstrated that the complete plastome has the potential to provide strong support for deep relationships in the plant tree of life, even when dealing with explosions of diversity over short periods of time (Barrett et al., 2016; Welch et al., 2016; Wei et al., 2017; Wei and Zhang, 2019; Foster et al., 2018; Nie et al., 2020; Pereira et al., 2021a; Larsén et al., 2022). Furthermore, plastome sequences provide new opportunities to explore the relationships among lineages with ambiguous or unresolved phylogenies (Wei et al., 2017). Previous studies of Isoëtes species based on plastome data have been reported, and some questions regarding their evolutionary history and phylogeny have been discussed (Karol et al., 2010; Schafran et al., 2018; Santos et al., 2020; Suissa et al., 2020; Pereira et al., 2021a). However, we believe that the species diversity of quillworts might be seriously underestimated, and many cryptic species might be misidentified in China. This belief led us to search for cryptic speciation in Chinese Isoëtes species by re-identifying population ploidy levels, re-observing spores, analyzing haplotypes, and reconstructing phylogenetic trees of almost all Chinese Isoëtes populations.

2. Materials and methods 2.1. Taxon sampling

A total of 46 samples of Isoëtes were used in this study: 26 were collected from the field (Table 1), four were cultivated in the greenhouse (Isoëtes sp.8 was in Lushan Botanical Garden, Jiangxi; I. yunguiensis 2 in Harbin Normal University, Heilongjiang; Isoëtes sp.4 in Guangxi Normal University, Guangxi; Isoëtes sp.13 in Zhongshan Botanical Garden, Jiangsu), and three were obtained from herbarium specimens [Isoëtes sp.1 was from South China Botanical Garden Herbarium (IBSC), and I. japonica A. Br. and Isoëtes sp.12 were from China National Herbarium (PE)]. A total of 13 complete plastome sequences from Isoëtes species were downloaded from GenBank (http://www.ncbi.nlm.nih.gov/genbank/), including three accessions (I. sinensis, I. yunguiensis 3, and I. taiwanensis) from China and 10 accessions from South America (Table 1).

2.2. Spore observation

The morphology of spores was examined under a scanning electron microscope (Quanta 250 FEI, USA). Untreated plant spores were collected after maturation and dried under natural conditions. Spores were mounted on a double-sided adhesive tape attached to metal stubs, sputter-coated with gold, and viewed under a Quanta 250 system at an accelerating voltage of 25 kV (kV) for megaspores and 30 kV for microspores. Previous studies were used as references to describe spore ornamentation (Hickey, 1986; Ranker, 1993; Watanabe et al., 1996; Wang and Yu, 2003; Wang and Dai, 2010; Holmes et al., 2005; Liu et al., 2008). To determine the average size of spores, we measured the length of the equatorial axis of 10–30 megaspores and microspores. Photographs are shown in Figs. 1 and 2, and the spore descriptions are presented in Table 2.

2.3. Chromosome number

Chromosomes were counted using the methodology of Zhang et al. (2019) but with some procedural modifications (see Table 3). Young root tips of the sporophytes were pretreated in saturated p-dichlorobenzene aqueous solution for 5 h and were then fixed in Carnoy's solution for 1 h at 4 ℃. They were dissociated in a mixture containing a ratio of 3% cellulase: 2.5% pectinase (1:1) for 15 min. The chromosomes from 20 samples were counted, and photographs were taken in the laboratory (Figs. 2 and 3). Three Isoëtes xiangfei Y.H. Yan, Y.F. Gu & J.P. Shu (Shu et al., 2022) samples from the same collection site had equivalent chromosome numbers. The chromosome number was determined using a line drawing on each photograph.

2.4. DNA extraction, library preparation, and plastome assembly

Isoëtes plastome library construction was done at Majorbio Corporation (Shanghai, China). Briefly, DNA was extracted from 33 leaf samples (fresh or silica-dried) and then sequenced for plastome library construction using the TruSeq Nano DNA HT Sample Preparation Kit (Illumina, San Diego, CA) following the manufacturer's recommendations. Specifically, each DNA sample was indexed with a unique marker for downstream sequence identification. The samples were divided into 350-base-pair (bp) fragments, which were end polished, A-tailed, and ligated with full-length adapters for subsequent Illumina sequencing and PCR amplification. The libraries were sequenced using the Illumina HiSeq X Ten platform (Illumina), and 150-bp paired-end reads were generated with insert sizes of approximately 400 bp. After enrichment of each sample, between 2 and 4 Gigabases (Gb) of sequence data were obtained.

To assemble the reads, raw reads were filtered using the default settings (-L: 5 -p: 0.5-N: 0.1) in ng_QC v.2.0 (Majorbio Corporation, Shanghai). Then the scripts and libraries in the GetOrganelle toolkit (Jin et al., 2020) were used to recruit target organelle reads from the sequence data, with Isoëtes yunguiensis used as the reference. The total target-associated reads underwent de novo sequencing assembly into a FASTG file using the St. Petersburg genome assembler (SPAdes v.3.13.0). The slimmed FASTG file was visualized using a Bioinformatics Application for Navigating De novo Assembly Graphs Easily (Bandage) (Wick et al., 2015) to finalize the complete target genome manually. Geneious Prime 2021.0.3 (https://www.geneious.com) was also used for data visualization and sequence correction.

2.5. Phylogenetic analyses for the complete plastid genome

Phylogenetic analyses of Isoëtes were conducted by both Maximum likelihood (ML) and Bayesian inference (BI). First, we aligned the complete plastome sequences of 46 Isoëtes samples using Mauve Aligner v.2.4.0 implemented in Geneious Prime. We used an advanced algorithm and assumed collinearity, excluding poorly aligned regions from the complete plastome data set, using Gblocks v.0.91b (Gerard and Jose, 2007) implemented in PhyloSuite v.1.2.2 (Zhang et al., 2020a) with DNA as the type of sequence (-t) up to half gap positions enabled (b5: half), and other parameters within the default settings. The resulting alignment was first trimmed for phylogenomic analysis to exclude misaligned regions and positions with > 90% missing data. ML analyses were performed using IQ-TREE v.1.6.12 (Nguyen et al., 2015) with 10, 000 bootstrap replicates. The best fitting model was selected by ModelFinder (Kalyaanamoorthy et al., 2017) and implemented in IQ-TREE. We visualized ML trees and ML bootstrap values (MLBS) using FigTree v.1.4.4 (Andrew, 2018); trees were rooted using Isoëtes nuttallii as the reference sequence in accordance with the methodology of Dai et al. (2020). BI analysis was performed using MrBayes v.3.2.6 (Ronquist et al., 2012) in PhyloSuite (Zhang et al., 2020a) under a GTR + I + G model selected by the Bayesian information criterion in ModelFinder. Analyses were performed using four parallel Markov Chain Monte Carlo (MCMC) runs for 20, 000, 000 generations, sampling every 1000 generations with a 25% burn-in. We combined simultaneous, independent runs to obtain majority rule consensus trees and to calculate Bayesian inference posterior probabilities (BIPP).

2.6. Variant calling with sequence data

All sequences, which we treated as low-depth, were mapped to the Isoëtes xiangfei reference genome (Shu et al., unpublished data) using the Burrows-Wheeler Alignment Tool (BWA v.0.6) with default settings (Li et al., 2009). Picard-tools v.1.92 (https://picard.sourceforge.net) was used to assign read group information that contained the sample library number, lane number, and sample identity. BWA was also used to sort the Sequence Alignment/Map (SAM)-format files, remove reads marked as duplicates, and generate Binary Alignment Map (BAM)-format files. Raw single nucleotide polymorphisms (SNPs) were filtered with VCFtools (Danecek et al., 2011) (using the parameters -minMapQ 30 -minQ 20 -max-missing 0.98) to obtain high-quality SNPs. To infer the individual ancestry proportions, we used ADMIXTURE v.1.23 (Alexander and Lange, 2011) with cross-validation for numbers of genetic clusters (K) from 1 to 20. The optimum number of clusters (K) was determined using the cross-validation errors (CV-error).

2.7. Haplotype and network

DNA Sequence Polymorphism (DnaSP v.6) was used to identify the haplotypes in the plastid sequences (Julio et al., 2017), producing an RDF format (.rdf) file. Network 10.2 was employed to analyze the.rdf file and draw the network (Bandelt and Dress, 1992).

2.8. Estimation of divergence times

The divergence times were estimated using the BEAST v.2.6.6 program (Bouckaert et al., 2019). The GTR + F + I model for the combined data set was applied with a Yule speciation tree and a relaxed clock log normal clock model after. We generated two calibrated notes by inputting the results from Dai et al. (2020) into the TimeTree online platform (http://www.timetree.org/). The age of the root was set to 'normal distribution, ' and the treeModel.rootHeight parameter was set to normal. BEAST trees were run for 200, 000, 000 generations and were sampled every 1000 generations. Tracer v.1.7.1 (http://beast.bio.ed.ac.uk/Tracer) was used to check whether the effective sample size (ESS) parameter values were greater than 200. The final annotation was completed using TreeAnnotator v.1.8.0 after discarding the first 10% of generations.

3. Results 3.1. Spore characteristics

We examined spores from 24 Isoëtes samples, including megaspores and microspores. The megaspores of Isoëtes samples were globose and trilete, with an equatorial ridge dividing the proximal and distal hemispheres and three radial ridges converging at the pole of the proximal hemisphere, where the surface is laevigate or has various types of ornamentation (Figs. 1 and 2; Table 2). Microspores were ellipsoid, monolete, with a smooth or textured surface, ornamented with spines, granules, or tubercles (Figs. 1 and 2; Table 2).

Megaspores can be divided into eight types based on ornamentation characteristics. Megaspores are (ⅰ) laevigate in Isoëtes hypsophila; (ⅱ) reticulate in I. fengii (Gu et al., unpublished data) and Isoëtes sp.3; (ⅲ) echinate in Isoëtes sp.2 and Isoëtes sp.9; (ⅳ) cristate in Isoëtes sp.11; (ⅴ) echinate-cristate in Isoëtes sp.4, Isoëtes sp.6, Isoëtes sp.8, I. baodongii, and I. sinensis; (ⅵ) cristate-reticulate in Isoëtes sp.5, I. yunguiensis, Isoëtes sp.7, I. orientalis, Isoëtes sp.10, I. xiangfei, and Isoëtes sp.12; (ⅶ) tuberculate-cristate in I. longpingii Y.H. Yan, Y.F. Gu & J.P. Shu (Shu et al., 2022) and Isoëtes sp.13; and (ⅷ) rugulate in I. taiwanensis. Microspores can be divided into five types. Microspores are (ⅰ) echinate in I. fengii, Isoëtes sp.4, Isoëtes sp.6, Isoëtes sp.8, I. hypsophila, Isoëtes sp.9, I. longpingii, I. orientalis, I. taiwanensis, I. baodongii, Isoëtes sp.13, and I. sinensis; (ⅱ) levigate-granulate in Isoëtes sp.3, Isoëtes sp.5, Isoëtes sp.10, and I. yunguiensis; (ⅲ) tuberculate-echinate in Isoëtes sp.7; and (ⅳ) laevigate in I. xiangfei. We assessed other studies for any missing information. The sizes of both megaspores and microspores are shown in Table 2.

3.2. Chromosome count

Chromosomes from 18 samples were observed, photographed, and counted with line drawings (Fig. 3). Four samples—Isoëtes yunguiensis, Isoëtes sp.10, I. baodongii, and Isoëtes sp.13—have 22 chromosomes; I. fengii, I. orientalis, and Isoëtes sp.3 have 66 chromosomes, and ten samples have 44 chromosomes. For samples in which chromosomes were uncounted in this study, we supplemented our observations with data from previous studies (see Table 2).

3.3. Genomic structure and phylogenetic analyses

We obtained new plastome sequences for 33 samples. Complete plastome lengths for these Isoëtes samples ranged from 145, 474 to 146, 363 bp (Table 1). The GC content of the plastomes was approximately 38.0% (Table 1). All 46 sequences used in this study presented strong collinearity; furthermore, Mauve Aligner identified only one locally collinear block without any rearrangement (Fig. 4).

Phylogenetic trees based on complete plastome sequences of 46 Isoëtes samples were reconstructed by ML (Fig. 5a) and BI (Fig. 5b). ML and BI phylogenetic trees were highly consistent except for some variation in the positions of four samples. The ML and BI analyses provided a well-supported phylogenetic structure, with most branches supported by > 90% MLBS values and 1.0 BIPP values (Fig. 5). All 35 samples from China and one from Japan clustered into a large monophyletic branch, in which I. hypsophila formed a monophyletic clade. Isoëtes fengii, together with I. japonica, also formed a monophyletic clade. The tetraploids did not form a monophyletic clade, as I. taiwanensis and I. orientalis were embedded in the branch formed by 14 tetraploid samples. Isoëtes baodongii and Isoëtes sp.13 were clustered with two other tetraploid samples (Isoëtes sp.2 and Isoëtes sp.6).

3.4. Genetic structure analysis

A total of 36, 378 SNPs was obtained and employed to estimate the optimal clustering number. The minimum CV-error was found when K was 19 (Fig. 6A). Almost every sample presented a single structure when the genetic structure was calculated with SNPs (Fig. 6B). Introgression occurred only in Isoëtes longpingii and two I. hypsophila samples. I. xiangfei and Isoëtes sp.4 clustered into the same structure, and Isoëtes sp.2 presented the same structure with Isoëtes sp.9 but not with Isoëtes sp.11, which clustered with I. baodongii. The genetic structure differed between I. yunguiensis distributed in Yunnan Province and I. yunguiensis distributed in Guizhou Province. Similarly, the two samples collected in Yunnan (Isoëtes sp.10 and Isoëtes sp.12) had different genetic structures. Samples with different ploidy on the phylogenetic tree were well separated (Fig. 6C). I. hypsophila formed a monophyletic clade; I. japonica is the sole member of one clade; hexaploid species I. fengii formed a monophyletic clade together with I. yunguiensis and another two Yunnan samples (Isoëtes sp.10 and Isoëtes sp.12); another two hexaploidy species were clustered at a different location in a large clade that also contains two diploid samples (I. baodongii and Isoëtes sp.13) and all tetraploid samples.

3.5. Haplotype and network analysis

Six Isoëtes hypsophila samples shared one haplotype (Hap22); two I. fengii samples shared one haplotype (Hap20); three I. xiangfei samples shared one haplotype (Hap1); three Isoëtes sp.7 samples shared one haplotype (Hap6); four I. yunguiensis samples shared one haplotype (Hap13); two Isoëtes sp.9 shared haplotype 2 (Hap2); each of the other samples possessed one haplotype.

When we used these shared haplotypes in network analysis, we found no net structures (Fig. 7). However, a total of 10 median vectors (mv) were produced. Isoëtes hypsophila and I. fengii connected to mva; Isoëtes sp.10, Isoëtes sp.11, and I. yunguiensis connected to mvb; Isoëtes sp.6, Isoëtes sp.12, Isoëtes sp.2, and I. baodongii shared mvc then connected to mvb; Isoëtes sp.1, Isoëtes sp.5, and I. longpingii shared mvf then connected to mve (which only Isoëtes sp.8 connected); no haplotype connected to mvd directly; Isoëtes sp.3, Isoëtes sp.7, I. xiangfei, and I. orientalis shared mvg; Isoëtes sp.4 and Isoëtes sp.9 shared mvh which connected to mvg; I. sinensis connected to mvi; I. taiwanensis and I. japonica shared mvj. A long line segment between Isoëtes hypsophila and the median vector was identified, and I. fengii, Isoëtes sp.1, and I. japonica also showed a long line segment to the median vector they connected with.

3.6. Expansion and contraction of the inverted repeat (IR) region

High gene arrangement synteny was observed in 35 Isoëtes plastomes from China and Japan (Fig. 4C). Differences in plastome size of Isoëtes can be primarily attributed to the extent of the IR (inverted repeat) region for each individual species (Fig. 4C). All Isoëtes plastomes within this study contained the SSC (small single copy region)/IRB boundary within the ycf2 gene and the SSC/IRA boundary between ycf2 and ycf1 genes. Two types of plastomes were characterized by IR length variation. The type I plastome was identified in I. hypsophila samples, which were characterized by an expanded IR region (reaching a length of 13, 691 bp), and the other 30 samples were identified as a type Ⅱ plastome (ranging from 13, 195 to 13, 207 bp).

3.7. Molecular dating

Divergence time for each sample was estimated using the BEAST v.2.6.6 program based on plastome data (Fig. 8). The divergence of Isoëtes distributed in South America and the samples from this current study likely occurred in the Middle Paleocene (~57.79 Ma; 95% HPD = 43.21–71.76 Ma). The first branch to diverge was Isoëtes hypsophila, at about 48.05 Ma (95% HPD = 30.78–66.31 Ma), in the early Eocene. The divergence times for the remaining 30 samples were estimated to occur between 3 and 20 Ma.

4. Discussion 4.1. Phylogenetic tree reveals complex relationships between different ploidy species

Our phylogenetic analysis results showed that the Chinese samples (not including I. japonica) formed a clade, largely consistent with the results of Choi et al. (2018). All the Isoëtes distributing in East Asia clustered a monophyletic clade together with samples from New Guinea and East Australia. Several clades received low support values and some taxa are resolved differently in the ML and BI analyses, but most clades are well supported especially at the first lineage of the two phylogenetic trees. These findings further indicate that our inferred phylogeny can serve as a reliable framework for exploration in the classification of the Chinese Isoëtes diploid–polyploid complex. The species diversity is much higher than we know and cryptic species in this genus are common in China. In our sampling we treated all samples (except for the type species from the original locality) as unknown species, thus, we can avoid preconceived identification of species. Samples used in this study have the widest coverage of Isoëtes and they can represent the taxa in China. The phylogenetic framework can be used to assess the species taxonomy and to reevaluate the delimitation of this family.

Historical events known as ancient hybridization (Peters et al., 2007) and paleopolyploidization (Barker et al., 2008) may have occurred before the existence of presently known Isoëtes species. This may be the reason why we got the chaotic phylogenetic tree, samples with different chromosome number are clustered together. When diploid samples were separated out to reconstructed phylogeny relationship (Fig. 9), we aquired almost the same topology as the original phylogeny tree employed with all the samples and mainly four clades were well supported of Chinese samples. As ancestor group, these diploids played important roles in the speciation of polyploid species in China. Of course, hybridization between diploids cannot be ignored either. Isoëtes sinensis had been known as formed by I. taiwanensis and I. yunguiensis (Liu et al., 2004; Taylor et al., 2004; Dai et al., 2020) but been queried by Kim et al. (2010), with new tetraploid (Shu et al., 2022) and diploid species (Lu et al., 2021) were discovered, our phylogenetic tree indicates that the speciation of this tetraploid should be reconsidered, I. baodongii may also be involved in the formation of tetraploid species including I. sinensis. Besides allopolyploid, there should be autopolyploids in this genus, the unique genetic structure of each sample in the genetic structure analysis may suggest this possibility. Hybridization, migration and isolation between diploids have led to their present genetic diversity and species diversity, all known and unknown diploids formed other polyploids by hybridization and polyploidization, so that all these samples present an intricate phylogenetic structure.

4.2. Multiple lines of evidence indicate high species diversity in Chinese Isoëtes samples

Despite the advent and acceleration of the genomics area, defining what makes a species genetically and morphologically distinct remains a problem in many fields of study. Canonically, a speciation event is defined as when an individual organism (ⅰ) forms a monophyletic branch when reconstructing the phylogenetic tree (Maddison, 1997; Maddison and Lacey, 2006; Liu and Pearl, 2007), (ⅱ) presents a different ploidy (Troia et al., 2016) causing reproductive isolation from other species (Fontaneto et al., 2015; Fier et al., 2017), and (ⅲ) possesses a morphological feature markedly different from other species.

Previous studies of Isoëtes in China more or less ignored the ploidy of the samples collected from different places (Pang et al., 2003; Ye and Li, 2003; Chen et al., 2004; Taylor et al., 2004; Kang et al., 2005; Wang et al., 2006; Liu et al., 2008; Xie et al., 2019; Dai et al., 2020). Thus, nearly all the populations from Anhui, Hunan, Jiangsu, and Zhejiang have historically been treated as I. sinensis. The total number of chromosomes has been found to vary among species (Hickey, 1984; Bhu and Goswami, 1990; Liu et al., 2002), and the physical appearances of the mature megaspores were first emphasized by Edward (1922) and are still considered an important taxonomic feature (Schafran, 2019).

The hexaploid species was reported later but is by no means a simple group. Besides Isoëtes orientalis, a hexaploid species, two other hexaploid samples were collected from Taining County (Isoëtes sp.3) and Chuxiong City (I. fengii). Although Isoëtes sp.3 is a sister clade to I. orientalis in the plastome tree, its position is different in the phylogenetic tree based on SNPs (Figs. 6C and 10), and the megaspore morphological features presented dissimilarly to its sister clade, revealing that Isoëtes sp.3 should be considered a new species. Isoëtes fengii formed a solely distinct clade both in the phylogenetic tree based on plastome data and SNP data, confirming that it is a new species. The tetraploid I. sinensis is a historically complex species in China since its megaspore ornamentation can be divided into six types rather than just one (Table 2). The spore sizes of tetraploids often vary, the megaspore size of Isoëtes sp.9 is 400–460 μm (mean = 435) and is the largest, the megaspore size of I. longpingii is 310–350 μm (mean = 325) and is the smallest of the tetraploids. Moreover, all tetraploids appear in nine clades on the phylogenetic tree instead of forming one monophyletic clade (Figs. 4 and 5). By integrating spore morphology and molecular phylogenetic evidence, all tetraploids could be identified as at least six species. The SNP data supported this result; even though some samples' positions changed, they were still separated into five clades.

Our analysis found that diploids did not necessarily cluster into one clade. Some were embedded in different branches of the phylogenetic tree, whereas others clustered with other tetraploid or hexaploidy samples. Four diploids of Isoëtes yunguiensis, I. taiwanensis, I. hypsophila, and I. baodongii have been previously reported in China (Handel-Mazzetti, 1923; Devol, 1972; Wang et al., 2002; Li et al., 2019; Lu et al., 2021); however, here we identified a new diploid sample of Isoëtes sp.13. Morphologically, megaspores of Isoëtes sp.13 were smaller than those of other diploids and had different surface ornamentation (Table 2; Fig. 1). However, Isoëtes sp.13 had the same genetic structure and clustered with I. baodongii. Therefore, we treated Isoëtes sp.13 as I. baodongii. Although four I. yunguiensis samples from Guizhou clustered into one clade, they did not form a monophyletic clade unless together with two Yunnan samples (Isoëtes sp.10 and Isoëtes sp.12). Even more notably, the two Yunnan samples (Isoëtes sp.10 and Isoëtes sp.12) possessed different genetic structures. According to these results, it may be inaccurate to identify the quillworts distributed in Yunnan Province all as I. yunguiensis.

There was some conflict between the phylogenetic tree based on plastome data and the phylogenetic tree based on SNP data (Figs. 8 and 9), which supposed to conflicts between the phylogenies of plastome and nuclear data. These conflicts between different data sets may have been caused by incomplete lineage sorting and hybridization/introgression (Funk and Omland, 2003; Yi et al., 2015; Wei and Zhang, 2022). The network result does not present any net structure, which may be caused by the single parent (maternal) genes. The structure of the network is almost consistent with phylogenetic tree reconstructed based on plastome sequences, besides Isoëtes japonica connected the same mv with I. taiwanensis with so long genetic distance between them. Even though all the samples have short branch lengths on the phylogenetic tree, their haplotypes showed different genetic distances. This result suggests a high genetic diversity in these samples.

As more and more new distribution sites are found in China, we suspect that Isoëtes distribution in China should be universal. Just as Troia et al. (2016) reported, several regions in Asia appear to be without Isoëtes, mostly because of insufficient exploration rather than an unsuitable climate for growth. This region has a potentially high level of species diversity because cryptic species in these unexplored regions were easily ignored, especially for morphologically simple and similar plants inhabiting comparable environments. Few species within this region have been reported and identified because of insufficient field studies and a lack of appropriate research.

4.3. Cryptic speciation of quillworts in China

Historically, we have defined and considered three main pathways for speciation—allopatric speciation, neighborhood speciation, and sympatric speciation. Congruently, isolation, differentiation, and genetic variation have been the main mechanisms we have identified that catalyze the formation of new species (Abbott et al., 2013). Recent studies have found that natural hybridization and polyploidy are the main mechanisms of speciation (Mallet, 2007; Rieseberg and Willis, 2007; Abbott et al., 2013). As Darwin stated in On the Origin of Species (Darwin, 1951), all species living in different places have a common ancestor, who migrated to different places and adapted to new geographical conditions, went through natural selection, accumulated mutations, and evolved into new species. Simpson (1961) challenged this viewpoint and postulated that almost any taxon can appear to be a center of origin if it steadily spreads in all directions until encountering impenetrable barriers and gradually begins to contract (often splitting into small disjunctive areas), creating new species. Geographical isolation can prevent gene exchange between populations, and long-term geographical isolation leads to genetic isolation among individual populations (Rieseberg and Willis, 2007).

Isoëtaceae are likely to have evolved from Annalepis (Pigg, 1992; Meng, 1998), with the family becoming widespread during the Triassic period because global sea-level movements enabled them to travel long distances with ocean currents. With the Tethys Sea flooding from west to east through Yunnan and Sichuan, Isoëtes expanded their range to the entire Yangtze River basin (Meng et al., 2000). Finally, their distribution range has gradually shrunk from east to west because of the periodic regression from east to west in the Yangtze region (Meng et al., 2000). Previous studies have identified. The Asian quillwort is thought to have originated from Annalepis in the Yangtze River Valley before the uplift of the Qinghai-Tibet Plateau, which occurred about 50 Ma (Eocene) (Harrison et al., 1992; An et al., 2001). This study estimates the node age showed that the Chinese samples and I. japonica were divided from the South American samples around 57.79 Ma (Paleocene) before I. hypsophila occurred. This period may be when the Isoëtes species originated in East Asia. The uplift of the Qinghai-Tibet Plateau resulted in western China being higher than eastern China, and the Yangtze River flowed from higher regions to lower regions; thus, the spread of Isoëtes species appeared to flow along with water to low elevation regions. Plants inhabiting high-elevation areas are subjected to low temperatures and high-intensity ultraviolet rays, which create species variants. For example, the expanded IR region of I. hypsophila resulted in a longer chloroplast genome with longer non-coding regions than other related species.

4.4. River system isolation promoted quillworts speciation

With the uplift of the Qinghai-Tibet Plateau, a change in topography could split a continuous water system into isolated systems, which limited gene flow among populations (Liu et al., 2004). When mapping the sites of Isoëtes distribution in China to the water system, each species sample corresponded to a independent water system mainly concentrated in the Yangtze River basin (Fig. 11). Water flow is the primary mode of plant dispersal for Isoëtes (Gentili et al., 2010; Troia et al., 2016) because their sporophylls can be carried by flowing water. This means the species dispersal occurs from higher to lower elevations in a single direction unless other media help to spread the sporophylls. Liu et al. (2004) reported that when isolated spores from populations located in high elevation regions or the headstreams of mountains came together in low elevation regions via the downstream movement of water, there would be a greater opportunity for speciation in low elevation regions. On the upon points, we think that when a new species occurs in a low elevation area or lower reaches of the river, it will be difficult found in a high elevation area or upper reaches of the river. Consequently, polyploid species occurring in Hunan, Jiangxi, Guangxi, Zhejiang, Anhui, and Jiangsu cannot be found in the Qinghai-Tibet Plateau and Yunnan-Guizhou Plateau. Almost all the samples in this study grown in the upper reaches of a river or near the end of tributaries, but it's difficult to find them in the downstream. We observed that these quillworts inhabited locations with gently flowing water so they could not be washed away, and these locations with slow-moving water also facilitated reproduction.

The distribution of species with different ploidy is affected by the geographical environment. For Isoëtes, diploids might have a greater tolerance for environmental extremes than polyploids (Liu et al., 2004); therefore, we can find diploids in the Qinghai-Tibet Plateau and Yunnan-Guizhou Plateau, where they have stronger ultraviolet rays and a lower temperature. However, polyploids could be more competitive within their geographical range (Liu et al., 2004), which may be why we find polyploid species in almost all areas of low elevation along the Yangtze River. It appears that with the uplift of the Qinghai-Tibet Plateau, sporophylls of diploids were swept away by floods from high elevation areas, aggregating in low elevation areas. However, it is difficult to explain the occurrence of hexaploid I. fengii found in Yunnan-Guizhou Plateau and diploid I. baodongii inhabiting low elevation areas. Isoëtes fengii and I. japonica formed a monophyletic clade, which was before I. yunguiensis and I. taiwanensis, suggesting that I. yunguiensis was not involved in the formation of I. fengii, and that I. hypsophila should be its mother. As I. fengii is a hexaploid species, an unknown tetraploid species should be involved in its speciation. At low elevation area, Zhejiang Province may be a central area for the speciation of quillworts living in the third terrain step geographically, where the speciation of tetraploids and hexaploids is related to the diploids, such as I. taiwanensis and I. baodongii. With this in mind, future studies are needed to explore these notable knowledge gaps.

4.5. Disappearing cryptic species calls for a new biodiversity conservation strategy

Cryptic species are at great risk of extinction, with several cryptic species reported extinct each year. Cryptic species, whose prevalence impairs biodiversity estimations, often challenge taxonomists tasked with their identification. Cryptic species seems common in fern taxa, in the last two decades, many cryptic species have been reported, such as Aspleniaceae (Jiang et al., 2011; Xu et al., 2019a, 2019b, 2022), Hypodematiaceae (Li et al., 2018), and Dicranopteris (Wei et al., 2021). Especially some endangered groups of quillworts (Pereira et al., 2019b; Schafran, 2019; Lu et al., 2021; Shu et al., 2022), Ceratopteris (Zhang et al., 2020b; Yu et al., 2022), Ottelia spp. (Ito et al., 2019; Li et al., 2020). However, cryptic species substantially add to our overall understanding of biodiversity, and, thus, require increased conservation efforts (Krug et al., 2013; Nygren, 2014; Pante et al., 2015; Struck et al., 2018). Because dispersal and reproduction in Isoëtes plants depend on water (Gentili et al., 2010; Troia et al., 2016), human activities that lead to poor water quality, habitat loss, and invasion by exotic species have increased the number of Isoëtes species listed as threatened worldwide (Chen et al., 2005, 2008; Kang et al., 2005; Kim et al., 2008; Gentili et al., 2010). Furthermore, human activities have decreased the distribution ranges and sizes of Isoëtes populations in China, and have led to the disappearance of some Isoëtes species from their known habitat locations (e.g., I. sinensis, I. yunguiensis, I. longpingii, and unnamed species Isoëtes sp. 1 and 12). Recent research has also predicted that human activities will continue to decrease potential suitable habitat for Isoëtes plant species (Yang et al., 2022).

Extinction is inevitable, but present extinction rates are exceptionally high because of anthropogenic activities (Pimm et al., 2014). It has been posited that approximately 1 million out of 8.7 million species are thought to be under threat of extinction (Tollefson, 2019). This number is still considered an underestimation because it is based on extinction risk estimates for all known species and does not consider whether undiscovered species may have higher extinction risks than species that have already been described (Liu et al., 2022). Among the anthropogenic changes on earth, habitat conversion represents the leading cause of species extinction, and agricultural activities are the drivers of the destruction of natural environments (Gonalves-Souza et al., 2020). According to the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, 1 in 9 species are threatened by extinction (Tollefson, 2019). Additionally, the global warming crisis has led to increasing biodiversity concerns, and some studies suggest that 12% of species face extinction (Malcolm et al., 2010). Other studies estimate that 7–24% of plant species face extinction-level events (Vuuren et al., 2006).

How do we reduce biodiversity loss and conserve wildlife? These are significant challenges humanity faces that require global participation. Current conservation efforts often focus on species discovered centuries ago or species identified as having significant economic value, excluding newly described and cryptic species. Many unknown, cryptic species become extinct before they are even discovered (Scheffers et al., 2012; Lees and Pimm, 2015). The Chinese government is exploring a set of feasible plans for the protection of cryptic species (40 individual taxa), including Huperzia, Ceratopteris, Cyatheaceae, Angiopteris, Cycas, Podocarpus, Amentotaxus, and Taxus, in which all species, known or unknown, have been listed in the checklist of nationally protected plants. Although the number of protected plants falls short of recommendations of the Convention on Biological Diversity (CBD)—at least 10% of species in each ecological area will be protected (Coates, 2016), contributing to global biodiversity conservation. It is easier and more effective to protect a taxonomic group than just one species, particularly for taxa that include cryptic species. Wildlife conservation should be a long-term and continuous process, requiring the participation of all countries to reduce the global extinction rate and increase biodiversity conservation efforts.

Acknowledgments

We want to thank Jian-Ying Xiang (Southwest Forestry University), Jian-Yun Wang (Baoshan University), Xue-Liang Hou (Xiamen University), and Shi-Wei Yao (Zhongshan Botanical Garden) for providing Isoëtes materials. We would also like to thank Wen Shao (Shanghai Chenshan Botanical Garden) and Guo-Hua Zhao (Fairy Lake Botanical Garden) for their help scanning spores. We also would like to thank TopEdit (www.topeditsci.com) for linguistic assistance during manuscript preparation. This study was supported by the Key Laboratory of National Forestry and Grassland Administration for Orchid Conservation and Utilization (grant number OC202103), the Harbin Normal University Postgraduate Innovation Project (grant number HSDBSCX2021-01), the National Natural Science Foundation of China General Projects (grant number 32170216), and the Hangzhou Science and Technology Development Project (grant number 20201203B113).

Authors contributions

Yu-Feng Gu: Data curation, Formal analysis, Investigation, Methodology, Wrote the original manuscript draft. Jiang-Ping Shu: Methodology, Formal analysis, Wrote the original manuscript draft. Yi-Jun Lu and Hui Shen: Investigation, Funding acquisition. Wen Shao: Scanning spores. Yan Zhou and Qi-Meng Sun: Investigation, Resources. Jian-Bing Chen, Bao-Dong Liu, and Yue-Hong Yan: Funding acquisition, Methodology, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.