In this study, we sampled a total 88 material, including 21 individuals of suspected hybrid M. matthewii and its potential parents (19 individuals of M. hancei and 44 individuals of M. marginata) based on the hypotheses outlined in the introduction section of this study and findings from previous research (Wang, 2016; Luo et al., 2018). M. marginata includes the original variant, var. calvescens (sometimes treated as M. calvescens in some literature, e.g. Liu et al., 2016; Wei and Zhang, 2016; Fraser-Jenkins et al., 2017), var. bipinnata, with an additional four individuals of M. strigosa as outgroups based on previous study (Luo et al., 2018). Most samples (60 individuals, 68.2% of the total) were collected at the population level, while the remaining samples were obtained from previous field collections stored in the National Wild Plant Germplasm Resource Center for Shanghai Chenshan Botanical Garden (ZWGX2202). The sampling locations are represented on the map ( Figure 1I ), and further details regarding the sampling are listed in Supplementary Table S1 .

These specimens were field-collected and dried using silica gel. Voucher specimens were curated at Chenshan Herbarium (CSH). DNA extraction from silica gel-dried leaves followed the modified CTAB method, with subsequent determination of DNA concentration and total amount. Samples passing the test proceeded to double digest Genotyping-by-sequencing (dd-GBS) library construction, involving double-enzymatic (Msel-Taqal) digestion and adapter attachment to the cleavage site ends of digested fragments. Primers designed based on adapter sequences were used for fragment amplification. After passing quality control, bipartite sequencing was conducted on the Illumina HiSeq X Ten platform (Illumina, USA), capturing approximately 150 bases of sequence on either end. Raw data underwent adapter trimming and quality screening to obtain clean data. Reads of low quality, including those with more than 40% of nucleotides having a quality value lower than 15, more than 10% of N nucleotides, or a length less than 30bp, were discarded.

The chloroplast and nuclear SNP datasets were obtained separately. To call chloroplast SNPs, we used the published plastid genome of M. marginata (GenBank assembly accession: MT130649; Du et al., 2021) as the reference. The clean data of 88 samples were then mapped to this reference using BWA v.0.7.12-r1039 (Li and Durbin, 2009) with maximal exact matches. The resulting BAM files were directly used in the standard GATK Best Practices workflow to call variants via the GATK HaplotypeCaller algorithm. Default GATK program settings were applied to filter the variants, and InDel variants were excluded from our study. In the end, we obtained 4649 chloroplast SNPs.

Nuclear SNPs were called by assembling the retained clean reads using the de novo assembly pipeline Stacks v2.1 (Rochette et al., 2019), and before calling, reads that could be mapped to the chloroplast reference were excluded. We required at least six identical reads (option –m 6) to create a stack for each individual using ustacks. A catalog of all loci across populations was then constructed using cstacks. After matching each sample against the catalog with sstacks, tsv2bam and gstacks were executed to incorporate paired-end reads, identify, and phase the SNPs. The variant call format (VCF) of called variants was exported using populations. The final dataset for subsequent analysis was filtered by vcftools with a minor allele frequency of 0.01, a maximum missing rate of 0.5, and a minimum depth of five. In the end, we obtained 1274 nuclear SNPs.

Phylogenetic trees were constructed from a nuclear SNP dataset using maximum likelihood (ML) concatenation analysis. The best-fit nucleotide substitution models were selected based on the Bayesian information criterion (BIC) using ModelFinder (Kalyaanamoorthy et al., 2017) implemented in IQTREE v2 (Nguyen et al., 2015). IQTREE was executed with the ascertainment bias correction (+ASC) model (Lewis, 2001) for the nuclear SNP dataset, and 1000 ultrafast bootstrap replicates were performed.

Population structure analysis was performed using ADMIXTURE v1.3.0 (Alexander et al., 2009). The number of pre-defined genetic clusters ranged from K = 1 to K = 7, determined based on the number of species and presumed genetic components, and adjusted according to the final results. The best-fitting model was determined using cross-validation (CV) error.

Plink v2.0 (Purcell et al., 2007) was applied to conduct principal component analysis (PCA) on the nuclear SNPs. Each sample was plotted based on the first two principal components (PCs), and the resulting figure was generated using R 4.1.3 with the ggplot2 package (https://ggplot2.tidyverse.org).

To assess historical introgression between species in the nuclear genome, we employed Patterson’s D-statistic (ABBA-BABA statistics; Durand et al., 2011) in the Dsuite program (Malinsky et al., 2021) and gene frequency covariance in the treemix software (Pickrell and Pritchard, 2012).

For Treemix analysis, the vcf file was obtained using vcftools and converted to the Treemix format using a Python script (https://github.com/wk8910/bio_tools/tree/master/03.treemix, accessed on 7 May 2023). Then, the analysis was conducted incrementally, considering pre-defined numbers of migration events (1-5) using the software Treemix. The increase in explained variation was compared for each number of migration events to infer the best result. The tree for the best result was visualized in R 4.1.3.

In Dsuite analysis, we assessed admixture across the M. hancei to M. marginata complex utilizing Patterson’s D (Green et al., 2010; Durand et al., 2011) and the f-branch (fb) statistic (Malinsky et al., 2018) with Dsuite v0.4r38 (Malinsky et al., 2021). To determine Dmin, representing the minimum allele sharing for all trios of ingroup lineages (n = 5), we employed Dtrios from Dsuite with the SNP dataset and the treemix phylogeny, irrespective of any assumptions about the tree topology. The fb statistic with Fbranch from Dsuite summarized rates of introgression. The generated “.tree” file by Dtrios and the treemix phylogeny were employed, and resulting fb statistics were plotted on the treemix phylogeny using “dtools.py”. M. strigosa served as the out-group.

Phylogenetic analyses of the chloroplast SNP dataset were conducted employing three different methods: maximum parsimony (MP), maximum likelihood (ML), and Bayesian Inference (BI). These analyses were performed using PAUP * 4.0b10 (Swofford, 2003), RAxML-HPC (Stamatakis, 2006), and MrBayes v3.2.5 (Huelsenbeck and Ronquist, 2001; Ronquist et al., 2012), respectively.

The ML analyses were performed using RAxML-HPC (Stamatakis, 2006) with ML tree searches and bootstrapping. The default model of -m GTRCAT was applied, along with 1000 rapid bootstrap analyses. A one-time search for the best-scoring tree was conducted (Stamatakis et al., 2006).

The MP analyses were carried out in PAUP* ver. 4.0b10 (Swofford, 2003) with 1000 tree-bisection-reconnection (TBR) searches. One thousand replicates were performed, each with 10 TBR searches, and a maximum of 100 trees were held per TBR search.

The best-fitting likelihood model for Bayesian analyses was chosen based on the Bayesian information criterion using jModeltest2 (Darriba et al., 2012). Bayesian inference (BI) was performed using MrBayes 3.2.5 (Huelsenbeck and Ronquist, 2001; Ronquist et al., 2012). We executed two independent runs, each employing four chains – one cold and three heated. The temperature parameter was set at 0.2, and the transition/transversion rate ratio was designated as beta, while the priors remained at their default settings. At the outset of each run, a random tree was initialized, and subsequently, every 1000 generations, a single tree was sampled, totaling 10,000,000 generations. The evaluation of convergence and stationarity was conducted using Tracer version 1.4 (Rambaut and Drummond, 2007), and to ensure the convergence of runs, the initial 25% of trees were discarded as burn-in. Subsequently, the remaining trees were employed to compute posterior probabilities (PP) for the majority-rule consensus topology.

The G-test was used to test the parentage ratio of the 21 individuals of the suspected hybrid M. matthewii (whether the ratio of M. calvescens: hancei significantly deviation from the expected 1:1 ratio). This analysis was conducted using the RVAideMemoire package (Herve, 2023) in R 4.1.3.

Furthermore, Bayesian Species Delimitation analyses (BFD) were conducted using the chloroplast SNP dataset to assess whether M. calvescens is a distinct species compared to other M. marginata variants, without considering hybridization. The nexus file was used to create the BFD input xml files in BEAUti v.2.4.5 (Bouckaert et al., 2014). We tested five different species assignment models, either combining the variants of M. marginata differently or treating them as distinct units. The BFD analyses were executed in SNAPP v.1.2.5 (Bryant et al., 2012), utilizing 12 initialization steps and chain-lengths of one million for each clade model. The remaining settings followed the guidelines set forth by Brandrud et al. (2020).

On the nuclear SNP phylogenetic tree ( Figure 2A ), the samples of M. hancei, M. marginata (excluding M. calvescens) and M. calvescens formed monophyletic groups separately (therefore, in the following text, any reference to M. marginata will exclude M. calvescens). However, M. matthewii did not form a monophyletic clade.

The PCA analysis demonstrated a clear separation of the 84 individuals into four lineages on the first two axes (PCs), which accounted for 73.9% of the total variation (see Supplementary Figure S1A ). These lineages corresponded to M. hancei, M. marginata, M. calvescens, and M. matthewii (positioned between M. calvescens and M. hancei) ( Figure 2B ).

The admixture analysis confirmed the patterns observed in PCA ( Figure 2C ), especially at K = 3, effectively distinguishing M. hancei, M. marginata, M. calvescens, and M. matthewii as a mixture of M. hancei and M. calvescens. The result at K = 4 showed a division of M. marginata into two genetic components. However, based on the CV error, K = 5 was determined as the optimal value (see Supplementary Figure S1B ), leading to the further division of M. calvescens into two genetic components. Conversely, all M. marginata variations except M. calvescens grouped together in the K = 2 and K = 3 results.

The ABBA-BABA tests revealed that three out of the ten tested four-taxon phylogenies, considering triplets involving M. matthewii and M. hancei, showed a significant signal of introgression (P < 0.05, non-zero D-statistics) ( Table 1 ). D-statistics ranged from 0.88 to 0.92 for all the significant tests. Additionally, one gene flow from M. hancei to M. matthewii was identified on the f-branch (fb) statistic plot ( Figure 3A ).

The general topology of the Treemix ML tree ( Figure 3B ) was consistent with the phylogenetic relationship recovered from the other phylogenetic analysis. One migration event between the M. matthewii and M. hancei improved with the explained variation from the initial 96% (with zero migration events allowed) to 99.96% (with one migration event allowed). The variation of further migration events (two to four) are 99.98%, 99.99% and 100%. As a max explained variation improving, the scenario with only one migration event was chosen to be the best model.

By employing the chloroplast genome of M. marginata as a reference, a comprehensive set of 4,649 SNP loci was acquired. Subsequent tree construction using various methodologies, including maximum parsimony, maximum likelihood, and Bayesian methods, consistently revealed similar topologies ( Figure 4 ). Excluding outgroups, three well-supported monophyletic lineages were identified within the inner group. M. calvescens constituted a distinct monophyletic group, while the other variants of M. marginata formed a sister group without M. calvescens, representing a separate monophyletic lineage. The third monophyletic lineage comprised M. hancei.

The two most optimal models identified in the Bayesian Species Delimitation (BFD) analyses both supported the topology in which M. calvescens is considered a distinct species, separate from other M. marginata variants, as indicated by the maximum marginal likelihood estimates (MLE) shown in Table 2 . The MLEs of the first and second optimal models were very close (Δ<10), and significantly higher (Δ>1000) than that of the third-ranked model.

Within these lineages, four individuals in M. matthewii clustered in the clade of M. calvescens, while the remaining 17 individuals clustered in the M. hancei clade, resulting in a significant deviation from the expected 1:1 ratio (G test: G= 8.6618, P= 0.003249).

The three cluster analyses conducted in this study—phylogeny, structural analyses, and PCA—clearly depict a mixture pattern for M. matthewii, suggesting the possibility of it being a hybrid species. Further details are discussed below:

The result of admixture analyses for K = 2 didn’t agree with the phylogenetic tree, as the M. marginata complex separated first, and M. hancei, M. calvescens and M. matthewii clustered together in another group. In many studies, the first run of structure at K = 2 tends to be the earliest of the two major branches to be separated or not genetically mixed (e.g., Yi et al., 2023) and often exists as the optimal K (e.g., Chattopadhyay et al., 2019; Xue et al., 2021; Encinas-Viso et al., 2022) because there are generally much longer branch lengths between the two main clades than within them. However, in this study, M. calvescens did not cluster with its closest relative on the phylogenetic tree, M. marginata, but instead with M. hancei at K = 2. This may be attributed to strong gene flow between M. calvescens and M. hancei, causing the two, along with their hybrid offspring, to cluster into a single genetic component first. Additionally, the CV error of K = 2 is very high and not optimal, predicting the implausibility of the structure.

The outcome from K = 3 indicates that M. calvescens forms a unique evolutionary cluster on the structure diagram, suggesting it should be classified as a distinct species. The remaining variants of M. marginata, including the original variant, var. bipinnata, and so forth, appear identical. The phylogenetic tree verifies the separation of M. calvescens and the other M. marginata varieties as two independent sister groups, both forming monophyletic clusters. Additionally, the principal component analysis (PCA) results indicate that M. calvescens and other M. marginata variants constitute two distinct groups.

The results of K = 4 for M. marginata reveal two distinct genetic groups, implying further differentiation within the species M. marginata or among these variants. However, the genetic groups do not correspond to the variants, suggesting that the differentiation is incomplete or ongoing and there is no evidence of new species arising at the present stage. And the results of PCA also did not support differentiation, with all individuals clustered together closely. Although some individuals in M. marginata are slightly separated on the PCA plot, consistent with the results of K = 4, they are so close that M. marginata is still considered a clustered entity.

The results of K = 5 for M. calvescens also reveal two distinct genetic groups, but this is not supported by the PCA, where all the plots clustered together.

Although the optimal K value is five with the lowest CV error value, the CV error values for K values ranging from three to six are very close to each other ( Supplementary Figure S1 , 0.28398, 0.27836, 0.26708, and 0.27742). Based on the results of the PCA and phylogenetic tree, K = 3 appears to be the most appropriate choice. Most importantly, the findings from K = 3 to K = 5 analyses reveal two significant points: 1) M. calvescens constitutes a distinct genetic component separate from other M. marginata variants, consistent with both the phylogeny and BFD analyses; 2) The presence of a mixture of M. calvescens and M. hancei genetic components in M. matthewii suggests a potential hybridization event between these two species. This speculation is also supported by the fact that M. matthewii lies between M. calvescens and M. hancei on the PCA plot.

Hybridization, involving gene flow between species, can lead to a mixed genetic structure. However, it is not the only cause. Two other mechanisms, namely ghost admixture and recent bottleneck effects, can also result in a mixed pattern (Lawson et al., 2018). Therefore, classical methods such as phylogeny, structural analyses, and PCA are currently insufficient for fully identifying instances of species hybridization. Gene flow assays are necessary tools that must be employed to augment these methods.

Therefore, we employed two methods to calculate gene flow to test whether the mixed pattern of M. matthewii is caused by hybridization. One of them, treemix, is based on gene frequency covariance. The presence of gene flow is indicated if the actual value is less than the estimated value, as gene flow can reduce gene frequency covariance. In the treemix analysis, gene flow from M. hancei to M. calvescens was detected, and M. matthewii and M. calvescens were clustered as sister groups on the maximum likelihood tree. This same gene flow was confirmed by the Dsuit results, using the ABBA-BABA principle, finally confirming that M. matthewii is a hybridization between M. hancei and M. calvescens.

By employing phylogenetic analysis using chloroplast SNPs to determine the parentage of M. matthewii, we discovered evidence of bidirectional and asymmetrical hybridization. The chloroplast phylogenetic analysis showed that most hybrids clustered with M. hancei, while a few clustered with M. calvescens. This suggests that both M. hancei and M. calvescens could be maternal or potential paternal parents of M. matthewii. This different from the numerous fern hybridizations with defined parents proved in previous studies, such as A. meishanianum, a hybrid with A. menglianense as father and A. malesianum as mother (Zhang et al., 2014; Shang et al., 2016).Therefore, it is likely that M. matthewii is a product of bi-directional hybridization involving M. hancei and M. calvescens.

Furthermore, it was observed that despite the minimal contribution of plastid genetics, offspring resulting from the two types of hybridization still displayed morphological differences. For instance, M. matthewii, with M. calvescens as its mother, exhibited longer petioles (Luo, 2018). This bidirectional hybridization may have further enriched biodiversity. Similar to M. matthewii, samples of M. krameri are positioned on the branches of both M. hancei and M. marginata in the chloroplast phylogenetic tree (Luo et al., 2018). This positioning suggests a potential bi-directional hybridization between these two species. However, further verification is required through the collection of additional material.

Out of the 21 hybrid offspring detected, 17 had M. hancei as the mother, while only 4 had M. calvescens. This resulted in a ratio of 4.25:1, significantly deviating from the expected 1:1 ratio. This type of asymmetrical hybridization has been previously reported in hybridization studies of two North American Dryopteris hybrids. For example, a 7.1:1 ratio of carthusiana: intermedia genome inheritance was found for D. × triploidea, and a 3.6:1 cristata: intermedia ratio was detected for D. × boottii. These deviations were driven by an array of reproductive traits, such as archegonial neck canals and spermatozoids (Testo et al., 2015).