Worldwide relationships in the fern genus Pteridium (bracken) based on nuclear genome markers.

PREMISE: Spore-bearing plants are capable of dispersing very long distances. However, it is not known if gene flow can prevent genetic divergence in widely distributed taxa. Here we address this issue, and examine systematic relationships at a global geographic scale for the fern genus Pteridium.
METHODS: We sampled plants from 100 localities worldwide, and generated nucleotide data from four nuclear genes and two plastid regions. We also examined 2801 single nucleotide polymorphisms detected by a restriction site-associated DNA approach.
RESULTS: We found evidence for two distinct diploid species and two allotetraploids between them. The "northern" species (Pteridium aquilinum) has distinct groups at the continental scale (Europe, Asia, Africa, and North America). The northern European subspecies pinetorum appears to involve admixture among all of these. A sample from the Hawaiian Islands contained elements of both North American and Asian P. aquilinum. The "southern" species, P. esculentum, shows little genetic differentiation between South American and Australian samples. Components of African genotypes are detected on all continents.
CONCLUSIONS: We find evidence of distinct continental-scale genetic differentiation in Pteridium. However, on top of this is a clear signal of recent hybridization. Thus, spore-bearing plants are clearly capable of extensive long-distance gene flow; yet appear to have differentiated genetically at the continental scale. Either gene flow in the past was at a reduced level, or vicariance is possible even in the face of long-distance gene flow.

Spore‐bearing plants are dispersed by microscopic haploid propagules that can travel vast distances in wind currents. Thus, many species of ferns, lycophytes, fungi, and other such groups are theoretically capable of extensive long‐distance dispersal, some of which may lead to establishment. How much intercontinental gene flow occurs? Is this sufficient to prevent genetic divergence? And at what scale are we able to detect these processes in spore‐bearing plants? Here we address these questions in one of the most widely distributed genera on the planet: Pteridium (bracken). Pteridium Gled. ex Scop. (Dennstaedtiaceae) is a fern genus found in most temperate to tropical regions from almost 70 degrees north (Alaska, Scandinavia, Siberia) to almost 50 degrees south (Argentina, Tasmania, New Zealand) and from sea level to above 3000 m (GBIF Secretariat, 2018). Bracken is a noxious weed in many parts of the world, often invading grazing land and outcompeting forage species (Holm et al., 1997). Furthermore, Pteridium is toxic to humans and other animals; consumption can result in teratogenesis, carcinogenesis, thiamine deficiency, retinal degeneration, and acute poisoning (Alonso‐Amelot and Avendano, 2002; Gil da Costa et al., 2012). Nevertheless, bracken has been used as fuel, thatch, litter, compost, food, medicine, and in the production of potash for the glass, soap, and bleaching industries (Rymer, 1976; McGlone et al., 2005).

Following a landmark revision by Tryon (1941), Pteridium has generally been treated as monospecific, although previously many more species had been described at regional levels (e.g., Ching, 1940). In a later review, Page (1976) noted that more systematic work is needed at a global scale to understand the number of distinct entities within the genus worldwide, and this is still the case. The majority of brackens are diploid (2n = 104: summarized by Page, 1976; Marrs and Watt, 2006). Current treatments have converged to recognize two diploid species (Thomson, 2012; Zhou et al., 2014; Wolf et al., 2015; Schwartsburd et al., 2018). Pteridium aquilinum (L.) Kuhn is predominant in the northern hemisphere whereas P. esculentum (G. Forst.) Cockayne is predominantly southern hemisphere in distribution. Additional to the two widespread diploid species, we recognize two allotetraploid taxa resulting from hybridization of these two species: (1) P. caudatum (L.) Maxon from Central and South America (and parts of North America); and (2) P. semihastatum (Wall. ex J. Agardh) S.B. Andrews from southeastern Asia and northern Australia. Other allopolyploids and/or introgressants may also exist locally where the northern and southern diploid species are, or have been, in contact (Brownsey, 1989; Wolf et al., 2015).

Studies on any taxon at a worldwide scale face the challenge of sampling populations in such a way as to represent diversity effectively, and this is inherently difficult at a global scale. In Pteridium, the problem is exacerbated by inconsistencies in morphological variation and what appears to be phenotypic plasticity for many traits (Page, 1976). In P. aquilinum (Boodle, 1903; Bright, 1928; Tryon, 1941) and in P. esculentum (O'Brien, 1963; Schwartsburd et al., 2014, 2018), environmentally labile characters include features regarded as critical for taxonomic assignment, as well as occasional paedomorphy in the latter species. Furthermore, because of uncertainties in taxonomic delimitations, and because fern spores are so mobile, it is not unexpected to find intermediate phenotypes and apparent hybrids. Plants collected for the present study were subsequently grown at a single site in Sydney, Australia under standardized garden culture in raised tubs preventing rhizome cross‐contamination, to facilitate taxonomic diagnosis and morphometric comparisons. This included 72 accessions (Thomson, 2000) used in both the present study and earlier work by Der et al. (2009) for a global phylogeny based on chloroplast DNA.

The number and rankings of sub‐specific taxa within the diploid brackens are still contentious. In Pteridium aquilinum, both morphometric analysis and molecular evidence from phenetic cluster analysis using the arbitrarily‐primed polymerase chain reaction (A‐P PCR) and inter simple sequence repeats (ISSR) suggest four subspecies in North America, two in Europe and North Africa, two in sub‐Saharan Africa, two in Asia, and one in the Hawaiian Islands (Thomson et al., 2005, 2008; Thomson, 2008). All eleven subspecies of P. aquilinum are represented in the present analysis. Strikingly, each of these taxa falls into one or the other of two major phenotypic groups separated on a set of morphological, phenological, and physiological criteria delimiting the “latiusculum” and “aquilinum” morphotypes respectively (Thomson et al., 2005, 2008). These are features likely to affect growth rate, frond size, laminal subdivision, light penetration, and litter breakdown in relation to adaptation to local climate (Thomson et al., 2005; Thomson, 2008). All four North American diploid brackens have the “latiusculum” morphotype while one subspecies of each morphotype is represented in Europe, sub‐Saharan Africa, and Asia respectively (Thomson et al., 2005). Intermediates consistent with hybridization have been observed in zones of contact or parapatry between the two morphotypes (e.g., Rumsey et al., 1991; Verdcourt, 2000; Karlsson, 2001). Nuclear genomic comparisons using A‐P PCR showed that sequence similarity for contrasting morphotypes from the same geographical region is generally higher than for morphotypes of the same kind from other geographical regions, providing strong evidence of local gene flow between them (Thomson, 2008). The genetic basis of these morphotypes remains unclear pending detailed comparative nuclear genomic analysis, although simple allelic differences do correlate with smaller‐scale polymorphisms in other ferns (Thomson, 2008).

Many of the main groupings discussed above were identified in a global scale analysis of molecular variation focused on chloroplast genes (Der et al., 2009). The results again highlighted the divergence between northern hemisphere (Pteridium aquilinum) and southern hemisphere (P. esculentum) clades, a pattern further corroborated in the work of Zhou and colleagues (Zhou et al., 2014). At the subspecies level, two results obtained by the latter authors are particularly noteworthy. As already noted by Der et al. (2009), chloroplast sequences from P. aquilinum subsp. japonicum (Nakai) A. Love and D. Love and P. aquilinum subsp. pinetorum (C.N. Page and R.R. Mill) J.A. Thomson were not distinguishable, suggesting their synonymy. In contrast, genomic comparison using AP‐PCR and ISSR sequence data clearly separate these two forms (see Fig. 6 in Thomson et al., 2005). Secondly, the chloroplast sequence analysis presented by Zhou et al. (2014) did not separate the two sub‐Saharan forms from each other or from P. aquilinum subsp. aquilinum from western Europe. Additional nuclear genome studies are particularly desirable as chloroplast DNA is maternally inherited in ferns that have been examined (Gastony and Yatskievych, 1992; Vogel et al., 2014; Kuo et al., 2018). Thus, hybrid and polyploid samples usually carry only one chloroplast haplotype and so only the maternal parent can be inferred. Conversely, recent hybrids usually retain nuclear DNA markers from both parents. Here, we addressed the issues of Pteridium systematics and hybridization at a global scale, using data from two distinct approaches. We used Polymerase Chain Reaction (PCR) to amplify chloroplast and nuclear genomic regions. We also prepared a genomic library with double‐digest restriction site‐associated DNA markers (ddRADseq). The objectives of our study are: (1) to infer distinct genetic entities in the genus and their geographic distributions; (2) to infer the number of diploid species; (3) to infer origins of tetraploid samples; and (4) to explore cases of genetic admixture resulting from dispersal between geographic regions.

MATERIALS AND METHODS

We sampled Pteridium from most named taxa (Fig. 1, Appendix 1) except for some recently described forms of P. esculentum from South and Central America (Schwartsburd et al., 2018). We used 101 samples for ddRADseq and 96 for amplicon analysis. Seven samples in ddRADseq samples were not in the amplicon samples (141PWLT, 321KISA, 293CPHW, MALW, CP73, 040CONT, AZ01) and two samples in the amplicon set were not used for ddRADseq set (Wolf_1006 and Wolf_1007). Total genomic DNA was extracted from tissue that had been silica‐dried, freeze‐dried, or pickled in CTAB/NaCl/ascorbate following Thomson (2002), or we extracted from silica‐dried material using the DNeasy Plant Mini kit (QIAGEN, Germantown, Maryland, USA), following the manufacturer's protocol.

Figure 1

Map showing location of samples used in this study. See also Appendix 1.

Amplicon sequencing

We selected two plastid genes (rpl16 and rps4) and four nuclear encoded genes (ApPEFP_C, CRY2, GapCp, and SQD) based on previous studies on Pteridium (Der et al., 2009; Wolf et al., 2015) and other ferns (Schuettpelz et al., 2008; Rothfels et al., 2013). PCR primer sequences were modified slightly (Table 1). Note that most of the plastid gene sequences are available on GenBank from an earlier study (Der et al., 2009); these were included for quality checking. Because sequence reads are shorter than the gene regions, we divided the genes into pieces (except for the targeted portions of CRY2 and GapCp which are <600 bp; Fig. 2) for a total of 12 amplicons (Table 1). Because amplicons were to be processed for sequencing on the Illumina platform, primers also included sequences that corresponded to the Illumina sequencing primers. Thus, all forward primers also included TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG at the 5′ end and reverse primers also included GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG at the 5′ end. We used a high‐fidelity DNA polymerase (Amplitaq Gold, Thermofisher Scientific, Waltham, Massachusetts, USA) to amplify the above 12 regions. Amplicons from all primer sets were then pooled for each individual and indexed during the standard library preparation. DNA sequencing was performed on a single run of the MiSeq platform with 2 × 300 bp paired‐end reads.

Table 1

Primers used for amplicon sequencing of Pteridium samples, based on previous studies (Schuettpelz et al., 2008; Der et al., 2009; Wolf et al., 2015). All forward primers included TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG at the 5′ end and all reverse primers included GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG at the 5′ end.

Marker name	Forward primer sequence	Reverse primer sequence
rpl16a	ACGCTTAGTGTGCGACTCGTT	GTTCAGATATTACTTGGTGC
rpl16b	GGTTAACCACTCAACTCCCG	TCCSCNATGTTGYTTACGAAAT
rps4a	ATGTCSCGTTAYCGAGGACCT	GGCTAGCAATGCGATTCACT
rps4b	TCTCTTTTGGGGGAAAAGATACC	AGAAGAGCGAAAGGGTTC
rps4c	CAGATTACTGAAAAACTAGC	TTACCGAGGGTTCGAATCCCTC
SQDa	GCAAGGGTACHAAGGTHATGATCATAGG	GCGTGARTCRTGCACTTTGCTRAGATG
SQDb	TGTGGGTGATATATGTGATTTTGAG	CCTTTGCCATAAACTGTAAGGGGGTG
ApPEFPa	ATACAGCTGCGGAAAATGCT	GCAATGCTGAAACAACCACA
ApPEFPb	ACACACACTATTGGGCACTGA	GAAGTGCAAAATACTAGAGCAGAAAC
ApPEFPc	ACTAGTTTGACAGAAACTCAGTACACA	TCAAAGGTGCGAGCTAACAA
CRY	AGGATGARYTGGAGAAAGGYAGCAATG	GTRTCCCAGAAATAYTTCATACCCC
GAPC	ATYCCAAGYTCAACTGGTGCTGC	GTAACCCAGAATGCCCTTCA

Figure 2

Relative locations of amplicons for each gene.

Amplicon analysis

Paired‐end sequences were first demultiplexed using the sample‐specific index. Sequences within samples were then paired and sorted by amplicon using PANDAseq with a strict quality threshold of 0.9 (Masella et al., 2012). We next used VSEARCH version 2.4.2 (Rognes et al., 2016) to cluster identical sequences within samples, while accounting for variation in read length. To remove possible sequence errors and PCR chimeras we used the usearch unoise algorithm (Edgar, 2016 [preprint]) to produce zero‐radius OTUs (zotus). This resulted in 595 amplicons with sequences. The remaining amplicons were clustered at 99%, or 98% identity using VSEARCH to retrieve sequences from amplicon‐sample combinations not in the zotus set. We estimated haplotype networks using the TCS method (Templeton et al., 1992), implemented in POPART (Leigh and Bryant, 2015). We attempted to phase haplotypes (Patterson et al., 2015) across overlapping amplicons, but could not achieve this for most samples due to a lack of single nucleotide polymorphisms (SNPs) within the regions of overlap, and a higher than expected number of haplotypes in some samples (see results).

ddRADseq

We generated a genomic library following a double‐digestion restriction site‐associated DNA sequencing (ddRADseq) protocol (Gompert et al., 2012; Parchman et al., 2012). We used EcoR1 and Mse1 to generate fragments of genomic DNA that were ligated to indexed oligonucleotides (with an index unique to each individual) on the EcoRI ends of the fragments. A standard, non‐barcoded oligonucleotide was ligated to the MseI ends of the fragmented DNA. Samples were then PCR‐amplified using iproof high‐fidelity DNA polymerase (New England Biolabs, Ipswich, Massachusetts, USA) with primers that overlap the ligated oligonucleotides. To reduce stochastic variation in PCR amplification, all fragments were first mixed with one other individual, which were further amplified in duplicate. The library was then reduced to fragments in the size range of 350‐450 bp using a Blue Pippin (Sage Science, Beverly, Massachusetts, USA). Quality and quantity was verified using TapeStation 2200 (Agilent Technologies, Santa Clara, California, USA). The size‐selected, multiplexed samples were run on a single lane of Illumina HiSeq 4000 with 100 bp single‐end sequencing at Genomic Sequencing and Analysis Facility at the University of Texas at Austin, Texas, USA (GSAF).

ddRADseq data processing

Processing the raw Illumina data was conducted in roughly three steps: (1) creation of a pseudo‐reference genome; (2) alignment of reads and variant calling; and (3) ENTROPY admixture analysis. Raw Illumina reads were demultiplexed and split by individuals using custom Perl scripts (Gompert et al., 2014). Since there is no reference genome available for Pteridium, we constructed a pseudo‐reference using the diploid taxa. We chose to use only the diploid taxa for the reference because they contain the majority of the sequence variation present in the tetraploids (known to be hybrid taxa between the diploids). Our pseudo‐reference was created by first clustering similar sequences for each diploid taxon separately using VSEARCH version 2.4.2 (Rognes et al., 2016); clustering was done at 92% similarity to create centroids for further clustering. This first stage attempts to combine different alleles from a locus, based on high levels of sequence similarity. However, similar paralogs within a diploid genome also tend to get included. Potential paralogs are removed, by eliminating loci that combine at a second clustering phase of at 84% sequence similarity. We then combined loci from the diploid species by selecting the intersection of all sequences from each diploid. After combining, we clustered using VSEARCH again at 84% similarity; the resulting sequences were then used as our pseudo‐reference genome.

Alignment of reads & variant calling

Before calling single nucleotide polymorphism (SNP) variants for all of the included species, we first had to index the pseudo‐reference genome from the previous step, which gives the sequence position points for the alignment. This was done using the INDEX function of BWA (Burrow‐Wheeler Aligner, version 0.7.10; Li and Durbin, 2009). Next we used PicardTools (website http://broadinstitute.github.io/picard/) to create a sequence dictionary; SAMTOOLS (Li et al., 2009) version 1.5 INDEX function was used to create a FASTA index file. For the alignment of reads we used BWA ALN to align all individuals to the pseudo‐reference. We then used the SAMTOOLS functions VIEW, SORT, and INDEX. To call variants we used the GATK (Van der Auwera et al., 2013) HaplotypeCaller (version 3.8.0) because it allows the user to specify ploidy. Because of this option, we called variants separately on the diploids and tetraploids. We used a custom Perl script to filter the resulting variant call format (VCF) files. We filtered for four stringency variables: (1) minimum coverage; (2) minimum number of sequences with an alternative allele; (3) loci fixed for the alternative allele; and (4) minimum mapping quality. Before combining the variants for the diploids and tetraploids we used custom Perl scripts to find the intersection of variants and then combine those intersecting variants from both files. Raw sequence data and details of amplicon analyses are available at website https://digitalcommons.usu.edu/all_datasets/67/.

ENTROPY analysis

We used the population genetics program ENTROPY (Gompert et al., 2014), which is very similar to STRUCTURE (Pritchard et al., 2000). Both are Bayesian, model‐based approaches to examine population structure, and both assume each individual's genome is a mosaic of loci inherited from K source populations. A key difference between them is that STRUCTURE assumes that individual genotypes are known, whereas ENTROPY uses genotype likelihoods calculated from raw sequence data and quality estimates. These genotype likelihoods are then used as an input for ENTROPY, whereas STRUCTURE calculates them repeatedly at each Markov chain Monte Carlo (MCMC) step from prior assignments (Gompert et al., 2014). The output graphic displays the population structure of samples at a range of source population numbers (K), with each sample coded by the proportion of the genome from each source. The first step of this analysis was to convert our VCF variant file to a Genotype Likelihood (GL) file format using a custom Perl script (see below). A second Perl script was used to convert the GL file to a matrix for input to R version 3.4.4 (R Core Team, 2008. We used the R package ADEGENET version 2.1.1 (Jombart, 2008) to perform a discriminant analysis of principal components to find the most likely source population for each individual. This analysis is similar to ENTROPY, but lacks much of the complexity. It is used to seed ENTROPY with starting values, which helps eliminate label swapping, and enables the Bayesian model to converge on the posterior more quickly. We followed the DAPC vignette (website http://adegenet.r-forge.r-project.org/files/tutorial-dapc.pdf), and chose a starting K value of 2.

Using our starting values, we ran ENTROPY for K = 2 through K = 15. Runs were 50,000 iterations with a burn‐in of 10,000 iterations and three independent chains per run. We chose to investigate those K values because K = 2 would be our most parsimonious explanation for the two diploid taxa contributing to the tetraploids; K = 15 would represent each sub‐species and tetraploid as its own separate population. We used a custom R script and functions to visualize our ENTROPY output, and chose K = 5 as the most informative. An R Markdown file with all custom scripts is located at website https://github.com/sylviakinosian/Pteridium_DDRAD-SEQ_Pipeline.

RESULTS

Amplicon haplotypes

We detected 1028 of the possible 1200 amplicons (i.e., 100 samples, each with a maximum of 12 possible amplicons). On average, we detected 1.57 haplotypes per plastid amplicon per sample. For nuclear amplicons we detected fewer haplotypes for presumed diploid samples (2.93 for Pteridium aquilinum, 2.64 for P. esculentum) than tetraploids (7.89 for P. caudatum, and 11.46 for P. semihastatum). However, a few amplicons had very high read depth and a corresponding high number of haplotypes, the highest of which was 49 haplotypes for GapCp in a sample of the tetraploid P. semihastatum (278KMAL). Thus, it is likely that at least some of retained haplotypes included both PCR and sequencing error. Removing sequencing and PCR error requires a tradeoff in losing possible true variants. Our approach has been more to the variant retention end of the spectrum, because we want to avoid losing variants associated with polyploidy and admixture. Broad scale relationships among haplotypes are illustrated by ApPEFP_C and SQD (Figs. 3 and 4), with remaining amplicon haplotype networks in supplementary data. In general, for each amplicon there are two dominant haplotypes, i.e., one from P. esculentum and one from P. aquilinum. All minor haplotypes emerge by a few steps (nucleotide substitutions) from these. Most of the presumed allotetraploids P. semihastatum and P. caudatum combine haplotypes of the two diploid species, or haplotypes that are within one or two steps (Figs. 3 and 4); however, some of the haplotypes found in the allotetraploids were unique to individual samples, and these may represent either haplotypes present in diploids that we failed to sample, or possible PCR chimeras.

Figure 3

Haplotype network for SQDa.

Figure 4

Haplotype network for ApPEFPb.

ENTROPY analysis of ddRADseq genotypes

We independently examined relationships among samples based on ddRADseq data. We detected 2801 SNP loci across diploid samples. With a few notable exceptions, at K = 2 the two species segregate and most presumed allotetraploid samples include contributions from each of the two diploids (Fig. 5). The exceptions include 040CONT (Pteridium esculentum from Tasmania), which grouped with P. aquilinum, and Wolf 921 (P. aquilinum from Connecticut, USA), which had some apparent contribution from P. esculentum at K = 2 (although at K = 5 this sample grouped with remaining North American P. aquilinum). Four of the six P. caudatum samples from Central and South America had equal contributions from P. aquilinum and P. esculentum. One sample from Bolivia (ij_2048) had mostly P. aquilinum at K = 2 (but see below at K = 5). The sample Wolf 1002 (P. caudatum from Galapagos) appeared to be P. aquilinum at K = 2 but a mix of P. aquilinum and P. esculentum at K = 5. All three P. semihastatum samples were an equal mix of P. aquilinum and P. esculentum.

Figure 5

Entropy patterns of genetic cluster contributions for K = 2 through K = 5 .

We detected more geographic and taxonomic resolution at K = 5 (Fig. 6), which appears to have the optimal signal, i.e., K > 5 did not appear to add any signal to the analysis. At K = 5. Pteridium aquilinum splits approximately into continents (Europe, Asia, North America, and Africa). However, P. esculentum in Australia and South America remain part of the same group, not splitting, even at high values of K. This result is consistent with treatment of P. esculentum as one species (Thomson, 2012). Within P. aquilinum, some subspecies appeared to involve admixture. Pteridium aquilinum subsp. decompositum (Gaudich.) Lamoureux ex J.A. Thomson appeared to include contributions from both North America and Asia, consistent with its location between the two, on the Hawaiian Islands. Likewise, P. a. subsp. pinetorum appeared to be a complex admixture involving contributions from Europe, Asia, North America, and even Africa. This is also somewhat consistent with morphology, distribution, and taxonomic treatments (Thomson, 2004). The most widely scattered genomic elements appear to be those of African P. aquilinum (subspecies centrali‐africanum Hieron. ex R.E. Fr. and capense (Thunb.) C. Chr.), components of which were detected in Europe, Australia, and one sample of P. caudatum (ij_2048) from South America.

Figure 6

Entropy patterns for K = 5.

DISCUSSION

Choosing the best approach for studying a group of closely related species and their derivatives is often challenging. Here we chose to use two very different data sets in an attempt to detect complementary signals. Amplicon sequencing exploits previously known phylogenetic markers to examine relationships among haplotypes that differ in nucleotide sequence, whereas ddRADseq samples many more genomic regions, in the form of single nucleotide polymorphisms. Filtering PCR error is challenging because an incorrect nucleotide incorporated at an early PCR cycle can result in a haplotype with high read depth (Andrews and Luikart, 2014). Furthermore, because polyploids possess extra gene copies, selection for maintaining extra functional copies may be reduced. This could lead to divergence, adding substantial noise to the analysis. We have been conservative with our filtering scheme, opting to risk retaining false haplotypes rather than eliminate possible real ones. However, we suspect that our amplicon data contain significant errors. Therefore, we have not submitted our amplicon sequences to GenBank, but instead have made them available from Digital Commons: website https://doi.org/10.26078/p64t-9h61. Despite extensive exploration, we are unable to determine the source of additional unexpected haplotypes, which may result from PCR error, despite our attempts to avoid it. We do not think that amplicon sequencing is an altogether ineffective method, although we suggest that longer read lengths, which eliminate the need for assembly and haplotype phasing, may result in higher success (Rothfels et al., 2017). Furthermore, various forms of target sequence capture of nuclear genes have been used very successfully at higher evolutionary levels (e.g., Uribe‐Convers et al., 2016; Dauphin et al., 2018). Here, we detected patterns consistent with previous findings in Pteridium, and the highest depth chloroplast haplotypes were identical to those obtained in earlier studies (Der et al., 2009; Wolf et al., 2015) using Sanger sequencing approaches. Conversely, ddRADseq has not been used previously for studies on Pteridium. We found that the ENTROPY approach revealed patterns than enabled a simpler interpretation of relationships among entities. Overall, we detected consistent patterns between the two approaches. Although both approaches entail similar levels of data analysis, ddRADseq requires considerably less lab work than amplicon sequencing. However, the latter can be streamlined by using alternative approaches to obtaining the data, such as using a target sequence capture protocol (Wolf et al., 2018).

Pteridium aquilinum versus P. esculentum

The strongest signal of genetic subdivision observed in both data sets globally distinguish P. aquilinum and P. esculentum. The same pattern has been observed with chloroplast‐only markers (Der et al., 2009) and on a more focused study on South American Pteridium (Wolf et al., 2015). The paradox surrounding this observation is that given the apparent high levels of admixture at the global scale and the fact that allopolyploids between these taxa have formed, how was it possible for P. aquilinum and P. esculentum to diverge in the first place. The division between the main Pteridium taxa is roughly that of a Laurasian P. aquilinum and a Gondwanan P. esculentum. We posit that these northern and southern groups were able to diverge because of reduced gene flow between them. Prevailing winds in equatorial regions are mostly easterly and westerly (Oort and Yienger, 1996), whereas air flows across the equator are relatively rare (Oort and Yienger, 1996; Wright et al., 2001). For example, the inter‐tropical convergence zone could enable spores to cross the equator, resulting in hybridization and subsequent allopolyploid formation in South America and in Australia/Southeast Asia (Thomson and Alonso‐Amelot, 2002). Another possibility is related to the conditions needed for spore production in Pteridium. In general, Pteridium sporulates more readily in dry conditions. If the climate was moist in the area where Pteridum was growing during the period of separation, then this could have contributed to accelerated divergence (via reduced gene flow) between P. aquilinum and P. esculentum.

Subspecific taxa

Two major features are apparent from the ENTROPY analysis of the ddRADseq data (Figs. 5 and 6). First, genetic distinctions within Pteridium aquilinum that reflect geographic distribution on a global scale are clearly evident at the subspecies level. Second, in agreement with previous studies (summarized by Thomson, 2008 and see Fig. 6 in Thomson et al., 2008), contrasting morphotypes from each continent are more similar genetically to each other than to the same morphotype from other continents.

In Asia, taxa of “aquilinum” and “latiusculum” morphotype are represented respectively by Pteridium aquilinum subsp. wightianum (J. Agardh) W.C. Shieh and japonicum which are seen as genetically close (Figs. 5 and 6 at K = 5,6), as are the sub‐Saharan African forms subsp. capense and centrali‐africanum. In western Europe, subsp. aquilinum is, with one exception, uniform genetically, i.e., all but one sample share a distinct genetic source. In contrast, subsp. pinetorum (“latiusculum” morphotype) shows elements of genetic similarity with subspecies latiusculum (Desv.) Hultén sensu Thomson, aquilinum, and japonicum. The ddRAD‐SEQ analysis supports distinction between the two sub‐Saharan African brackens and subsp. aquilinum (contra Zhou et al., 2014). In North and Central America, subsp. feei (W. Schaffn. ex Fée) J.A. Thomson, Mickel & K. Mehltreter, latiusculum, pubescens (Underw.) J.A. Thomson, Mickel & K. Mehltreter, and pseudocaudatum (Clute) Hultén all have a “latiusculum” phenotype. The close genetic similarity of the subspecies in each continental zone appears consistent with past interbreeding leading to reports of intermediate forms (e.g., Tryon, 1941; Page, 1976; Speer et al., 1999; Mickel and Smith, 2004).

In contrast to the genetic subdivisions observed in Pteridium aquilinum, we detected no clear differences between Australian and South American P. esculentum. This finding is more consistent with the treatment as a single species (Thomson, 2012; Zhou et al., 2014; Schwartsburd et al., 2018) rather than as separate taxa with P. esculentum in Australia and P. arachnoideum in South America (see e.g., Schwartsburd et al., 2014).

Origins of tetraploid taxa

The two tetraploid species, Pteridium caudatum and P. semihastatum, clearly show signs of admixture from the two diploid species in both our amplicon and ddRAD‐SEQ analyses. The only exception was one sample of P. caudatum from Bolivia that has a genetic composition more similar to African Pteridium. Otherwise, the presumed tetraploids were composed of approximately equal contributions from each of the two diploids. Several questions remain regarding the allotetraploid taxa: (1) How many independent origins have occurred? (2) Is allotetraploid formation ongoing? (3) Are tetraploids generated via sterile triploids? and (4) Can such triploids be detected? Evidence from P. caudatum suggests several origins and even possible backcrosses to parents (Wolf et al., 2015). Multiple origins of allotetraploid species is not an uncommon observation (see e.g., Soltis and Soltis, 1991; Zou et al., 2015). Several lines of theory and data suggest that tetraploids often arise in nature via a sterile triploid bridge in seed plants (Harlan and deWet, 1975; Bretagnolle and Thompson, 1995) and in ferns (Gastony, 1986). Initially, an unreduced spore is produced and the resulting diploid gametophyte is most likely to combine with a more common haploid gametophyte. The resulting triploid would be sterile except when it produces an unreduced (3x) spore that again, is more likely to combine with a common haploid gametophyte resulting in a tetraploid, which restores fertility so long as the combining genomes are balanced. To date, only a presumed autotriploid Pteridium has been observed (Sheffield et al., 1993), but allotriploids are likely to be common. The problem is that they are generally difficult to detect without cytological analysis. Future studies would benefit from sampling with a goal of detecting variation in genome size by flow cytometry, from which ploidy may be inferred.

Current patterns of admixture

Genetic groupings in Pteridium, even within P. aquilinum, were distinct at the continental scale. Furthermore, with a few notable exceptions, most samples within named taxa were genetically similar. However, we did detect several cases of possible admixture. One such case was the single sample of P. aquilinum subsp. decompositum from the Hawaiian Islands. The ENTROPY analysis of this sample indicated approximately three‐quarters contribution from North America, and one quarter from Asia, consistent with the location of the islands between these source continents. A similar situation is seen in P. aquilinum subsp. pinetorum, which is distributed in northern parts of Europe. Samples of this taxon have apparent contributions from Western Europe, Asia, and North America. This taxon has always been taxonomically challenging (Thomson, 2004). However, the proximity of P. aquilinum subsp. pinetorum (if looking from above the arctic) to all three source populations is consistent with the genetic results.

One of the major exceptions to consistency in the genetic results was the distribution of sub‐Saharan African genotypes. Samples from Pteridium aquilinum subsp. capense and subsp. centrali‐africanum grouped as a single source in the ENTROPY results at K = 5 (blue samples in Fig. 5). However, this genetic component was also detected in South America, Europe, and Australia. One possibility is that the African genotypes represent a relict southern entity. However, if this were the case then we would expect that admixture would have reduced the distinctness of the genotypes. Alternatively, the occurrence of African genotypes outside Africa more likely represents recent dispersal. If plants in Africa produce spores more readily (because of a drier climate) then it is not unlikely for them to enter airflows both above and below the equator. From there, reaching South America and Australia is not unreasonable. A previous study has detected North American genotypes in western Europe (Rumsey et al., 1991), providing evidence for recent trans‐oceanic dispersal.

Our findings are consistent with those of a chloroplast‐only study of many of the same samples (Der et al., 2009). However, inclusion of nuclear DNA markers enables us to see a more complete picture, including possible patterns of introgression and the origins of allotetraploid taxa. One possible issue that can confound analyses and make interpretations more challenging is possible variation in ploidy. In this study we assumed that taxa were either diploid or allotetraploid. In addition to alternative cytotypes, it is known that recent polyploids can have disrupted meiosis (e.g., Lim et al., 2008; Szadkowski et al., 2010; Chester et al., 2012). This can result in variation in the number of gene copies, a possible explanation for the high number of haplotypes observed in the amplicon data from some samples. In general, the results of the ddRADseq study were easier to summarize and interpret. Nevertheless, combining long sequences with phylogenetic signal, with large numbers of SNP data enables more insight into the evolutionary history of this group of ferns. We find evidence for two distinct (presumably diploid) clades with allotetraploids between them. Samples of the northern Pteridium aquilinum group approximately by continent, whereas the southern clade, P. esculentum, appears to be more homogeneous. Overlaid on this distinct geographic pattern is evidence for several cases of admixture between genetic (and geographic) groups. This is a similar pattern to that seen in scaly tree ferns (Korall and Pryer, 2013). We think that the treatment of Pteridium as two diploid and two allotetraploid taxa is warranted. Future studies would benefit from denser sampling in key geographic areas, especially near the equator and on oceanic islands. Screening for possible triploids, and other cytotypes, using flow cytometry would also provide useful information on the evolution of this ubiquitous group of plants.