Genetic and morphological identification of a recurrent Dicksonia tree fern hybrid in New Zealand

Hybridization is common in many ferns and has been a significant factor in fern evolution and speciation. However, hybrids are rare between the approximately 30 species of Dicksonia tree ferns world-wide, and none are well documented. In this study we examine the relationship of a newly-discovered Dicksonia tree fern from Whirinaki, New Zealand, which does not fit the current taxonomy of the three species currently recognized in New Zealand. Our microsatellite genotyping and ddRAD-seq data indicate these plants are F1 hybrids that have formed multiple times between D. fibrosa and D. lanata subsp. lanata. The Whirinaki plants have intermediate morphology between D. fibrosa and D. lanata subsp. lanata and their malformed spores are consistent with a hybrid origin. The Whirinaki plants–Dicksonia fibrosa × D. lanata subsp. lanata–are an example of hybridization between distantly related fern lineages, with the two parent species estimated to have diverged 55–25 mya. Our chloroplast sequencing indicates asymmetric chloroplast inheritance in the Whirinaki morphology with D. lanata subsp. lanata always contributing the chloroplast genome.


Introduction
Tree ferns (Cyatheales) are often a locally dominant feature of wet temperate Southern Hemisphere forests. Within New Zealand's forests tree ferns often represent a significant proportion of the forest community based on the number of individuals, biomass and basal area [1]. Two genera are represented in New Zealand's tree ferns, Cyathea Sm. (Cyatheaceae) and Dicksonia L'Hér (Dicksoniceae) [2,3]. The taxonomy of the New Zealand species in both of these genera is thought to be well understood, with no new species described since 1910 [2].
Dicksonia contains around 30 species that are found in Central and South America, Southeast Asia, eastern Australia, New Zealand and the Pacific [2,4]. Three indigenous Dicksonia species are recognized in New Zealand; all are endemic [2]. Dicksonia fibrosa Colenso (whekīponga) and D. squarrosa (G.Forst) Sw. (whekī) have trunked rhizomes up to 6m tall and 7m tall, respectively. Both are found throughout the North Island in lowland to montane forest and are largely confined to lowland and coastal areas in the South Island [2]. Two allopatric a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 1. the Australian species D. antarctica Labill., which is sometimes confused with D. fibrosa [2], including sometimes claimed to be present in New Zealand; 2. an undescribed species; 3. a homoploid hybrid between D. fibrosa and D. lanata subsp. lanata that has formed once and has spread vegetatively (both D. squarrosa and D. lanata produce underground stolons from their rhizomes [2]); 4. a recurrent homoploid hybrid between D. fibrosa and D. lanata subsp. lanata; 5. a fertile allopolyploid between D. fibrosa and D. lanata subsp. lanata.
If the Whirinaki morphology is either D. antarctica (hypothesis 1) or an undescribed species (hypothesis 2) then, given its restricted distribution, it would likely be considered threatened under the New Zealand Threat Classification System [7]. Hypotheses (3), (4) and (5) all involve hybridization, which is rare in Dicksonia, with only unconfirmed reports of wild Dicksonia hybrids [4,8] and a single recorded hybridisation event in cultivation between D. arborescens L'Hér. from St Helena and D. antarctica [9].

Sampling and DNA extraction
A total of 43 Dicksonia specimens were collected for this study (Table 1), with four to 14 individuals sampled for each of the New Zealand taxa, and one specimen of D. antarctica. Samples were collected under the Department of Conservation permit number CA-31615-OTH. Plants of the Whirinaki morphology were collected from two sites: c. 38˚48' S, 176˚37' E, and c. 38˚49' S, 176˚39' E. Fresh frond tissue was collected into silica gel. DNA was isolated from the dried tissue using a modified-CTAB extraction method (steps 1, 3-7 from Table 1 in [10]).

Chloroplast sequencing and analyses
Six chloroplast loci were tested for variation in a subset of the New Zealand Dicksonia samples.
Sequence files were edited with Sequencher version 5.2 (Gene Codes Corp., Ann Arbor, MI, USA). The trnL-trnF region demonstrated the most variation and was sequenced for a total of 33 Dicksonia specimens, as described above. The trnL-trnF sequences contained only two unambiguous indels and these were aligned by eye. A median-joining network [17] for the trnL-trnF sequences was produced using PopArt [18] with indels (insertions and deletions) recoded as single events.

Microsatellite analyses
Eleven microsatellite loci isolated from Dicksonia sellowiana Hook. [19] were tested in New Zealand Dicksonia. Primers to a further ten microsatellite loci developed from sequences obtained from a preliminary ddRAD-seq run on four Dicksonia (one each of D. fibrosa, D. squarrosa, D. lanata subsp. lanata and the Whirinaki morphology) were also screened. Primers were designed for these extra loci (S1 Table) using Primer3 [20].
An M13 tag (TGTAAAACGACGGCCAGT) was added to the 5' end of the forward primer of each locus. Loci were amplified individually in 10 μL PCR reactions that contained 1 μL of diluted template DNA, 0.2 μM forward primer, 0.8 μM reverse primer, 0.8 μM M13 primer (labelled with FAM, NED, PET or HEX) and 1× MyTaq mix (Bioline). PCR thermocycling conditions were an initial denaturation of 94˚C for 5 min followed by 35 cycles of 94˚C for 20 s, 52˚C for 20 s, and 72˚C for 20 s and a final extension at 72˚C for 15 min. For loci that amplified, genotyping was performed on an ABI 3730 DNA sequencer at the Massey Genome Service (Massey University, Palmerston North, New Zealand). Alleles were sized using the internal size standard GeneScan 500 LIZ (Applied Biosystems) and scored using Geneious version 10.2.6 (Biomatters Ltd., Auckland, New Zealand).

ddRADseq library preparation and sequencing
Double-digest restriction site associated sequencing (ddRAD-Seq) libraries [21] were prepared for 39 individuals with 8 of these processed in duplicate as technical replicates. For each sample, 300 ng of DNA was digested with two restriction enzymes (AvaII and MspI, New England Biolabs Inc), following manufacturer's instructions. These two enzymes were selected based on the recommendation of Yang et al [22], who found that this enzyme pair produced a high number of fragments for many plant species. Adaptors containing sample-specific barcodes and Illumina indices were ligated to each sample. Samples were pooled into three index pools and a size-selection performed on each pool for 300-500 bp fragments using excision from an agarose gel, followed by extraction with a Qiaquick gel extraction kit (Qiagen). Pooled size-selected samples were PCR-amplified to add Illumina indices using Phusion flash high fidelity PCR master mix (Thermo Scientific). Each pooled sample was purified and concentrated with a MinElute kit (Qiagen), quantified with a Qubit dsDNA HS (high sensitivity) assay kit (Thermo Fisher Scientific) and combined in equimolar amounts. A detailed protocol has been deposited in protocols.io (dx.doi.org/10.17504/protocols.io.zgyf3xw). The library was sequenced across a third of a lane using the Illumina HiSeq 2500 to generate 2 × 125 bp reads.

ddRADseq assembly
Raw paired-end reads were demultiplexed with ipyrad v0.7.28 [23]. All Fastq sequence files are available from GenBank at the National Center for Biotechnology Information short-read archive database (accession number: PRJNA522301). Four samples had low numbers of reads and were removed prior to assembly. For the remaining samples, reads were demultiplexed, had their adaptors removed, and were merged and assembled into de novo loci using ipyrad (the params file we used in iyprad is provided in S1 Text). Reads were clustered at 90% similarity, with a minimum depth of coverage of six. Samples were treated as diploid allowing two alleles per site. Only loci present in at least 50% of the samples were retained.
A preliminary NeighborNet network (see below) was constructed from this initial assembly to confirm that the eight technical replicates clustered with their duplicate. The two duplicated sets of reads were then pooled and assembled as described above. The D. antarctica specimen produced a high number of raw reads but these were assembled into few loci compared with the other samples so this sample was excluded from the final dataset used for the analyses described below.

ddRADseq data analyses
To examine whether there was conflicting phylogenetic signal in the ddRAD data, which may indicate a history of hybridization, a network was constructed with the NeighborNet algorithm [24], implemented in SplitsTree4 v4.14.6 [25] using uncorrected p-distances and the equal angle algorithm. Support with assessed with 1000 bootstrap replicates.
We used STRUCTURE v2.3.4 [26,27] to examine genetic structuring without a priori inferences. A single snp was randomly selected from each ddRAD locus and the number of genetic clusters (K) was set between 1 and 5, with 10 permutations for each. We used the admixture model with correlated allele frequencies and ran STRUCTURE with a burn-in of 100,000 generations followed by 500,000 Markov Chain Monte Carlo (MCMC) iterations. The optimal number of genetic clusters (K) was determined by calculating the ΔK statistic [28] in STRUC-TURE HARVESTER web v.0.6.94 [29]. We also examined all clustering results that warranted biological interpretation, following [30]. CLUMPP v.1.1.2 [31] was used to average iterative runs of K and the results visualized graphically with STRUCTURE PLOT [32].

Morphology
The morphology of the Whirinaki plants was compared to D. fibrosa and D. lanata, based on published descriptions [2,5]. The spores of 10 plants with the Whirinaki morphology were examined with a compound microscope. Spores were also examined from plants of D. fibrosa and D. lanata subsp. lanata.

Chloroplast sequencing
The trnL-trnF alignment was 972 bp in length and contained 17 substitutions and two unambiguous indel events. The relationships between the trnL-trnF sequences, with the indels coded, are shown in the median-joining network (Fig 1). The most common haplotype detected was shared by most of the D. lanata subsp. lanata, D. lanata subsp. hispida and Whirinaki morphology samples. Two additional haplotypes were found in samples with the Whirinaki morphology; each differed from the most common D. lanata haplotype by a single mutational change. Two haplotypes were restricted to D. lanata, one in each subspecies, and they differed from the common D. lanata haplotype by one or two substitutions. The two D. squarrosa samples had a haplotype that differed from the most common D. lanata haplotype by one mutation (a 10 bp insertion). The haplotypes detected in D. antarctica and D. fibrosa differed from the D. lanata haplotypes by a 6 bp insertion plus at least 12 substitutions. The D. antarctica trnL-trnF sequence was identical to one of the D. fibrosa sequences with the other two D. fibrosa sharing a haplotype that differed by one substitution.

Microsatellites
Three microsatellite loci amplified and were variable within and/or between New Zealand Dicksonia species (Table 1). Neither D. fibrosa (n = 6) nor D. antarctica (n = 1) showed any variation at the three loci. Dicksonia lanata (n = 14) only varied at the DicMic01 locus, with six alleles detected. A number of D. lanata subsp. lanata samples failed to amplify at this locus suggesting that null alleles may be present in this taxon. Dicksonia squarrosa (n = 4) was invariant at the DicMic109 locus, had two alleles at the DicMic104 locus and failed to amplify at the DicMic01 locus. At locus DicMic109 all the individuals of the Whirinaki morphology (n = 14) were heterozygous for the 131 bp and 125 bp alleles. The 131 bp allele was the only allele

ddRAD-Seq
Illumina sequencing of Dicksonia ddRAD libraries for all 38 individuals plus 8 duplicates resulted in 303.48M reads after initial quality filtering. Excluding the four low coverage samples, 668K to 7.9M reads per individual were obtained (mean±SD = 3.43M±2.10M reads per sample) ( Table 2).
The preliminary NeighborNet including the duplicate samples and the D. antarctica sample is shown in S1 Fig. It demonstrates that the duplicates largely cluster together and that the D. antarctica sample is closely related to D. fibrosa. The duplicates were combined for subsequent analyses and the D. antarctica sample was excluded owing to low coverage.
The final dataset comprised 1244 loci across 33 specimens. The NeighborNet analysis of this dataset clustered each species together with strong support (Fig 2). There was little variation within D. fibrosa and D. squarrosa. In contrast both D. lanata subsp. lanata and D. lanata subsp. hispida exhibited greater variation between individuals and while there was not strong evidence for the separation of these two subspecies, they were somewhat separated in the network.

Individuals of the Whirinaki morphology grouped together in an intermediate position between D. fibrosa and D. lanata.
For the Structure analyses ΔK indicated that the optimal K was 2 (S2 Fig; Fig 3). At K = 2 D. lanata and D. squarrosa were assigned to a single cluster and D. fibrosa was assigned to a second cluster. Individuals of the Whirinaki morphology were assigned equally to each of these two clusters. At K = 3 D. squarrosa was assigned with high probability to a third cluster while the Whirinaki morphology individuals were assigned equally between the D. lanata and D. fibrosa clusters. At K > 3 individuals were partitioned as for K = 3 but additional clusters were partitioned across all samples.

Morphology
About 80 individuals of the Whirinaki morphology were observed in the field in the Te Kohu area of Whirinaki Forest at an altitude of 800-920m above sea level. Plants with the Whirinaki morphology have a wide adventitious root-covered trunk similar to that of D. fibrosa but never exceeding 2m in height, whereas D. fibrosa is regularly taller (Table 3). Dicksonia fibrosa and the Whirinaki morphology also share skirts of dead fronds but the skirt is less well developed in the latter (Fig 4). The Whirinaki morphology is most easily distinguished from D. fibrosa by its longer stipes, which are usually orange-brown (only one of the 80 observed plants had green stipes), its wider fronds and the obtuse apices of its lamina segments (Table 3, Fig 5). Also, unlike D. fibrosa, the Whirinaki morphology has sparse brown woolly hairs up to 5mm long in the undersides of the rachis, pinna midribs and costae, especially in clusters in the coastae junctions. Similar hairs also occur in D. lanata subsp. lanata but the hairs are much denser than in the Whirinaki morphology. The spores examined from each of the ten accessions of the Whirinaki morphology were malformed and variable in size but generally smaller than those of D. fibrosa and D. lanata subsp. lanata (Fig 6). They also had no obvious exospore compared to D. fibrosa and D. lanata subsp. lanata.

The origin of the Whirinaki morphology of Dicksonia involves hybridisation
Our genetic results indicate that the plants with the Whirinaki morphology are not Dicksonia antarctica (hypothesis 1) because they differ in their chloroplast sequences and do not cluster together in the NeighborNet of the ddRAD-Seq data. Instead our analyses indicate that the Whirinaki morphology plants are hybrids between D. fibrosa and D. lanata. The NeighborNet of the ddRAD-Seq data placed samples of the Whirinaki morphology at an intermediate position between these two species and the STRUCTURE analyses assigned them approximately equally to the clusters containing D. fibrosa and D. lanata. Two of the microsatellite loci (Dic-Mic01, DicMic109) also provided support that the Whirinaki morphology plants are hybrids between D. fibrosa and D. lanata. Only the third microsatellite locus (DicMic104) did not show a pattern consistent with hybridization because D. lanata and the Whirinaki morphology were fixed for the same single allele and D. fibrosa was fixed for a different allele. None of the Dicksonia individuals genotyped were heterozygous at this locus and it is possible that this locus derives from a uniparentally-inherited organelle.     The two subspecies of Dicksonia lanata are only weakly separated in the ddRAD-Seq analysis, and not at all in the STRUCTURE analysis. This may reflect the make-up of the sample set, particularly the inclusion of hybrids. The AFLP analysis of a larger sample set by Lewis cited in [5] recovered strong support for separation of the two subspecies. This, together with morphological and geographic evidence, means we think there are strong grounds for maintaining them as separate taxa [2]. Furthermore, while the genetic results do not implicate one subspecies over the other in the hybridisation, only D. lanata subsp. lanata occurs in close proximity, and morphologically the Whirinaki plants show characters of both D. lanata subsp. lanata (cf. subsp. hispida) and D. fibrosa ( Table 3).
The abnormal spores observed in the Whirinaki morphology indicate that these plants are infertile F1 hybrids and the genetic analyses largely also support this conclusion. Particularly informative was the DicMic109 microsatellite locus where all individuals of the Whirinaki morphology were heterozygous for the same two alleles, one of which was fixed in D. lanata and the other was fixed in D. fibrosa. This is the expected allele pattern for F1 hybrids (homozygotes, as well as heterozygotes, would be predicted for other hybrid classes such as F2 and backcrosses). Furthermore no more than two alleles were detected per locus at the Genetic and morphological identification of a recurrent Dicksonia tree fern hybrid in New Zealand microsatellite loci, as would be expected for diploid hybrids. The high level of variation observed in the Whirinaki morphology in the ddRAD NeighborNet and at microsatellite locus DicMic04 indicates that this form has not arisen from a single hybridization event (hypothesis 3) but from repeated crosses between genetically differentiated parent individuals (hypothesis 4). The three chloroplast haplotypes found in the Whirinaki morphology also supports recurrent hybridization.

Asymmetric hybridization
The chloroplast haplotypes sequenced from the 15 hybrids were identical or very similar to the haplotypes detected from D. lanata, indicating that this species was the chloroplast donor. This 15:0 ratio is statistically different from the 1:1 expectation if D. lanata and D. fibrosa were equal contributors of chloroplasts to the hybrids (χ 2 = 15, d.f. = 1, p = 0.0001). Asymmetric hybridization has been noted previously in ferns [33][34][35][36] but the drivers remain poorly understood [35]. Factors such as the relative abundance of the parental species' gametophytes, differences between parent species' sperm size and dispersability and sizes of the archegonial neck canal, differential success of embryos from reciprocal crosses and whether parent species exhibit an antheridiogen system may contribute to asymmetric hybridization [35].

Hybridisation in Dicksonia
Our finding that the Whirinaki morphology comprises recurrent F1 hybrids was surprising for several reasons. Firstly, these Dicksonia fibrosa × D. lanata subsp. lanata plants are the first wild hybrids confirmed by genetics in the genus but they were locally common (although they represented only a small fraction of the tree fern community in the area). Despite being such a large plant, it is possible that Dicksonia fibrosa × D. lanata subsp. lanata plants occur elsewhere but have been overlooked. However, although the distributions of D. fibrosa and D. lanata broadly overlap over a wide area [2], the New Zealand National Vegetation Survey Databank suggests that the two species rarely grow in close proximity, which would limit their opportunities for hybridization. Of the 97555 permanent mainland plots in the New Zealand National Vegetation Survey Databank [37], D. fibrosa and D. lanata (either subspecies) are recorded from 1006 and 791 plots, respectively, but only occur together in 15 of these plots. Differences in ecological preferences may account for this lack of local overlap, with D. lanata subsp. lanata preferring hillsides and D. fibrosa more likely to be found on gully floors. The western Whirinaki area may be unusual in that the two species occur together over a large area, possibly because it is flatter, colder and drier forest than many other sites [38]. Examining other sites where D. fibrosa and D. lanata are sympatric may reveal further hybrids and provide insight into the environmental conditions that lead to hybridization between these two species.
Secondly, the hybridization has occurred between Dicksonia from chloroplast clades that have been estimated to have diverged 55-25 mya [6]. The Whirinaki morphology therefore provides another example of hybridization between deeply diverged fern lineages as has been found between Cystopteris and Gymnocarpium (~60 mya [39]) and Asplenium (45 mya [36]). Dicksonia squarrosa and D. lanata, which are widely sympatric and much more closely related to each other than either is to D. fibrosa [6], are not known to hybridise.
Lastly, although recurrently formed hybrids, Dicksonia fibrosa × D. lanata subsp. lanata exhibited remarkably conserved morphologies given the considerable differences between the morphologies of the parents.
Polyploidy is thought to be a significant process in the generation of fern diversity with an estimated 31% of fern speciation events accompanied by an increase in ploidy [40]. However, while many fern genera hybridize frequently and have highly duplicated genomes (e.g., Asplenium [41,42], others, such as Dicksonia, hybridize only rarely, leading to limited opportunities for polyploidy. Understanding the situations where the barriers preventing hybridization break down in rarely hybridizing genera such as Dicksonia may provide insight into the dynamics of fern evolution and a contrast to genera with frequent reticulate evolution.
We are not making a nothospecies name because the New Zealand practice for naming fern hybrids is to use a hybrid binomial i.e., Dicksonia fibrosa × D. lanata subsp. lanata.