{'fraction_non_alphanumeric': 0.040309651413128124, 'fraction_numerical': 0.02766707892446521, 'mean_word_length': 5.105145413870246, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Target sequence capture of nuclear‐encoded genes for phylogenetic analysis in ferns

Premise of the Study Until recently, most phylogenetic studies of ferns were based on chloroplast genes. Evolutionary inferences based on these data can be incomplete because the characters are from a single linkage group and are uniparentally inherited. These limitations are particularly acute in studies of hybridization, which is prevalent in ferns; fern hybrids are common and ferns are able to hybridize across highly diverged lineages, up to 60 million years since divergence in one documented case. However, it not yet clear what effect such hybridization has on fern evolution, in part due to a paucity of available biparentally inherited (nuclear‐encoded) markers. Methods We designed oligonucleotide baits to capture 25 targeted, low‐copy nuclear markers from a sample of 24 species spanning extant fern diversity. Results Most loci were successfully sequenced from most accessions. Although the baits were designed from exon (transcript) data, we successfully captured intron sequences that should be useful for more focused phylogenetic studies. We present phylogenetic analyses of the new target sequence capture data and integrate these into a previous transcript‐based data set. Discussion We make our bait sequences available to the community as a resource for further studies of fern phylogeny.

capable of hybridization among lineages that are highly diverged. One such example is hybridization across fern genera that shared their last common ancestor approximately 60 million years ago (Rothfels et al., 2015a), and additional, potentially comparable hybridizations have been reported in other groups of ferns (e.g., Wagner et al., 1992;Wagner, 1993;Larsson, 2014b;Engels and Canestraro, 2017). It is hypothesized that organisms that undergo abiotically mediated fertilization may lack the types of pre-mating barriers that are found in groups such as insect-pollinated angiosperms; thus, ferns might evolve reproductive barriers more slowly than angiosperms (Ranker and Sundue, 2015;Rothfels et al., 2015a). However, several questions remain: How common are such wide hybridization events in ferns, relative to other plant groups? And are the hybrids evolutionary dead ends? Such hybrids are probably always sterile. However, it is not known whether fertility can be restored via allopolyploidy across such widely divergent genomes, and if so, is there an upper divergence limit to such events? These types of evolutionary processes would be entirely missed with plastid-based studies, but analyses based on a sufficient number of low-copy nuclear-encoded genes should help detect signatures of such complex evolutionary histories.
Published studies using nuclear genes in ferns either focus on one or a few genes (e.g., Ishikawa et al., 2002;Pryer et al., 2004;Ebihara et al., 2005;Adjie et al., 2007;James et al., 2008;Schuettpelz et al., 2008;Chen et al., 2012;Sessa et al., 2012;Schneider et al., 2013;Rothfels and Schuettpelz, 2014;Rothfels et al., 2017) or genes from transcriptome data Li et al., 2015;Rothfels et al., 2015b;Shen et al., 2017). These approaches either scale poorly (PCR amplification of nuclear markers) or are expensive and restrictive with respect to the material that can be used (transcriptomics). Here we present a cost-effective approach for acquiring multi-locus nuclear-encoded gene data in ferns. We apply a target sequence capture approach (Mertes et al., 2011) for 25 low-copy genes, designed to work across the fern tree of life. This type of approach has been used successfully on other groups of organisms (Lemmon et al., 2012;Hart et al., 2016;Heyduk et al., 2016;Leveille-Bourret et al., 2018) and can use genomic DNA as a starting material. During the preparation of a genomic library, templates are filtered to enrich for fragments that match certain "bait" oligonucleotides, which are designed from independent data, such as transcriptome sequences. The advantages of this method are that time-and labor-intensive PCR is circumvented, templates can be genomic DNA rather than RNA, template quality conditions are less restrictive than they are for PCR-based approaches, large contiguous information-rich homologous sequences (including introns) can be captured, and larger sample sizes can be used. Our goal is to provide the fern systematics community with a cost-effective tool for generating nucleotide data from nuclear genes. Our approach is focused to exploit data already available for ferns.

Bait design
We started with alignments of 25 genes from a previous transcriptome study (Rothfels et al., 2015b) and designed baits from 33 taxa, spanning fern phylogenetic diversity, for which coverage was good (the list of taxa is available on Digital Commons [https://doi.org/10.15142/ t3mg95]). We used RepeatMasker version 4.06 (Smit et al., 2013) to flag simple sequence repeats and low-complexity DNA. From the remaining sequences, we designed 120-nucleotide baits with ~2× flexible tiling density, for a total of 20,353 unfiltered baits. We then collapsed baits that were 95% identical or higher (with at least 75% sequence overlap) into 19,863 clusters, from which baits were chosen, one from each cluster. We screened these sequences against 38 published fern chloroplast sequences (from GenBank); no baits were flagged as being similar to any of these sequences. Baits were designed and synthesized by Arbor Biosciences (formerly Mycroarray; Ann Arbor, Michigan, USA) and are available on Digital Commons (https://doi.org/10.15142/t3mg95).

Samples, library prep, and sequencing
We chose 24 fern samples to evaluate the efficacy of our approach (Table 1). Taxa were selected from across fern diversity with a slight emphasis on Dennstaedtiaceae, Lindsaeaceae, and some other early diverging Polypodiales, which were underrepresented in the earlier transcriptome-based work (Rothfels et al., 2015b). Thirteen of the samples were from genera used for bait design, and four were from families absent from the Rothfels et al. (2015b) study (Dicksoniaceae, Cystodiaceae, Saccolomataceae, and Metaxyaceae; Table 1). Samples were extracted from silica-dried material, except for Osmundastrum C. Presl, which was from a herbarium specimen. We used a cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1990) to extract DNA, which was then sent to Arbor Biosciences for further processing. Half the volume of each sample was purified using QIAquick PCR Purification columns (QIAGEN, Germantown, Maryland, USA) to remove residual contamination and eluted in 50 μL of EB buffer. Each purified sample was quantified with Quant-iT PicoGreen dsDNA Assay Kit (Molecular Probes Inc., Eugene, Oregon, USA), and up to 4 μg of total genomic DNA was taken through a sonication and bead-based size-selection protocol to achieve modal fragment lengths of approximately 600 nucleotides. Between 99 and 509 ng of processed DNA were then converted to TruSeq-style libraries and index-amplified with unique dual eightnucleotide indexes for six cycles with KAPA HiFi polymerase (454 Life Sciences, a Roche Company, Branford, Connecticut, USA). Each target enrichment reaction contained two libraries of similar library yield, ranging from 100 to 400 ng of each library per pool. Enrichment reactions followed the standard myBaits version 3 (Arbor Biosciences) protocol (see manual provided on Digital Commons [https://doi.org/10.15142/t3mg95]) using 500-ng baits, 62°C for 20 h of hybridization, bead-bait binding, and wash steps. After resuspension of bead-bound enriched libraries in 30 μL of recommended buffer, 15 μL were amplified for 10 cycles using KAPA HiFi polymerase, and the reactions purified with Solid Phase Reversible Immobilization (SPRI) beads (Applied Biological Materials, Richmond, British Columbia, Canada). The final enriched libraries were then quantified using library qPCR, combined equimolar, and passed to a single lane of MiSeq PE250 (Illumina, San Diego, California, USA). FASTQ files were delivered following standard post-processing pipelines in BaseSpace Sequence Hub (Illumina). The Culcita C. Presl sample failed on first attempt and therefore was added to a later sequencing run on a HiSeq platform, with 125-bp paired-end reads.
adaptive quality trim with a target length of 10 nucleotides and strictness of 0.2, and discarded any reads trimmed to under 20 nucleotides. We then used HybPiper version 1.3 (Johnson et al., 2016) to extract and sort genes for downstream phylogenetic analyses. This suite of programs uses BLASTX (Altschul et al., 1990) or BWA (Burrows-Wheeler Alignment tool; Li and Durbin, 2009) to assign reads to target genes, and reads are then assembled separately for each target using SPAdes (Bankevich et al., 2012). Our target file was based on the alignments from Rothfels et al. (2015b), with some modifications (e.g., removal of outgroup sequences) to better facilitate the sorting of reads to the appropriate genes. We used 33 taxa, across the fern phylogeny, for which we had most genes. For gene families with duplication events within ferns, one paralog was arbitrarily selected to include the pre-duplication copies (ApPEFP_AC for ApPEFP_B-ApPEFP_C; gapCpSh for gapCpSh-gapCpLg; CRY1 for CRY1-CRY2-CRY5; and CRY3 for CRY3-CRY4 [see Rothfels et al., 2013]) and the reference sequences modified accordingly. We ran HybPiper, using these reference sequences, in nucleotide mode; preliminary runs in amino acid mode were too inclusive in the read-sorting step: even distantly related paralogs were captured and grouped together. HybPiper annotates intron sequences using Exonerate (Slater and Birney, 2005) via a separate script (intronerate. py). Because our paired-end reads were approximately 500 bases, we reasoned that we should be able to capture up to that length of intron from the end of each exon (more if unbaited reads are included).

Phylogenetic analyses
The exon-only "paralog" sequences assembled by HybPiper were aligned for each locus in MUSCLE version 3.8.31 (Edgar, 2004). These alignments were manually refined while viewed as their translated amino acids in AliView version 118 (Larsson, 2014a), with ambiguous regions excluded from subsequent analysis. In several cases, the MUSCLE alignments resulted in long contiguous segments of sequence from individual accessions that were out of alignment compared to the rest of the sequences. BLAST results of these sequences demonstrated that they consisted of exons that were in the incorrect order, due to an Exonerate error that can occur when read coverage is low. The exons of these sequences were manually reassembled in the correct order (based on the BLAST coordinates), and the resulting sequences realigned. Preliminary maximum parsimony trees were inferred from these alignments in PAUP* v4.0 build 159 (Swofford, 2002) and, based on these trees, a small number of mis-sorted sequences were transferred to their correct locus (e.g., some CRY3 sequences were incorrectly sorted into the CRY4 pool). In the few cases where multiple sequences were recovered for single accessions within a locus, the longest sequence that was resolved in a phylogenetically appropriate location in the maximum parsimony trees was retained and any others erased; these "extra" sequences presumably reflect mis-assemblies, recent gene duplications, or homeologous sequences. These 25 individual-locus alignments were concatenated (while retaining exset and codon position information) into a single alignment with the "seqconcat" command in abioscripts version 0.9.3 (Larsson, 2010). This alignment was subject to maximum likelihood analysis in GARLI version 2.0 (Zwickl, 2006), with the data partitioned by codon position (see Rothfels et al., 2015b) and each data subset given an independent GTR+I+G substitution model. The best-tree searches were performed from 10 random-additionsequence starting trees, and support was assessed with 720 bootstrap pseudoreplicates (each searched from a single random-additionsequence starting tree). All GARLI analyses were performed on the CIPRES Science Gateway version 3.3 (Miller et al., 2010). The maximum likelihood tree was annotated with support values summarized from the 720 bootstrap trees using the SumTrees version 3.3.1 command in the DendroPy phylogenetic computing library version 3.12.0 (Sukumaran and Holder, 2010).
To further explore the phylogenetic signal in these data and their utility in combination with previously published data, we combined the 25 exon-only alignments with the corresponding transcript alignment from Rothfels et al. (2015b). We used a custom Python script to profile-profile align the alignment pairs (each Hyb-Seq exon-only sequence paired with its transcriptome-based partner) for each locus, using MUSCLE version 3.8.31 (Edgar, 2004). These alignments were then manually refined, concatenated, and analyzed as above for the Hyb-Seq-only alignments (Swofford, 2002;Zwickl, 2006;Larsson, 2010Larsson, , 2014aMiller et al., 2010;Sukumaran and Holder, 2010).

RESULTS
Raw reads, target sequence files, sequence alignments, taxon information, and analysis details are available on Digital Commons (https:// doi.org/10.15142/t3mg95). In addition, the intron-containing "supercontigs" for each sample are deposited in GenBank (MG817735-MG818130), and the exon-only alignments are deposited in TreeBASE (study S22073 [https://treebase.org/treebase-web/search/ study/anyObjectAsRDF.rdf?namespacedGUID=TB2:S22073]). We averaged 1,071,998 reads per sample for the MiSeq samples; the one HiSeq sample returned 15,570,735 reads (Table 1). We measured success of the baits by the proportion of the average gene length in the target file represented in the recovered exon-only sequences. By this measure, we recovered an average of 51% coverage (across all genes and all accessions; Fig. 1). However, there was considerable variation across genes and taxa. Recovery efficiency was poorest for Phlebodium (R. Br.) J. Sm., Hymenophyllum Sm., Saccoloma Kaulf., and Metaxya C. Presl and the genes NDUFS6 and COP9. We did not see a clear pattern of taxon success relative to genera used for baits, although the results suggest a combination of bait-related and taxon-specific factors. For example, Metaxya and Saccoloma are both phylogenetically distant from any sample included in the bait set, whereas Hymenophyllaceae species, while included in the bait set, have sequenced poorly in other broad analyses (e.g., Rothfels et al., 2015b). The mean per-base read depth for exon regions was 31.2. However, this includes regions that failed to get captured or were absent in a particular taxon; mean depth excluding samples with zero depth was 50.1.
Average supercontig lengths varied from 700 (NDUFS6) to 6356 (IBR3) bp (mean = 3426; Table 2). The mean length of (partial) intron sequences captured was 267 bp. For all introns detected, including those with no coverage, the mean was 135 bp. On average, sequences of introns and intron fragments made up 41% of the supercontig data; these intron sequences could be useful in future studies with more focused phylogenetic samples, but were effectively unalignable at the phylogenetic depth of this study.
Phylogeny estimates from our Hyb-Seq data were generally well supported (Fig. 2) and are broadly consistent with previous studies based on both nuclear and chloroplast genes (e.g., Schuettpelz and Pryer, 2007;Lehtonen, 2011;Rothfels et al., 2015b;Testo and Sundue, 2016). The analyses of our data combined with those from Rothfels et al. (2015b) find the same major relationships as in the earlier study, with some refinements, and with novel inferences for taxa not previously included in broad nuclear phylogenies (Fig. 3 and Discussion).

DISCUSSION
Our results demonstrate that our target bait collection is effective on a broad selection of fern taxa across the tree of life. We did not see a clear pattern of capture success based on relationships between taxa used in bait design and taxa sampled for capture. Nor is there a clear pattern with respect to genome size, for those few taxa for which we have estimates. We suspect that the main factors that affect success are DNA quantity and quality. Because the baits are designed from exon regions, they are more effective at capturing exon regions. However, the long read lengths enabled us to capture considerable (41% of the supercontig sequence) intron data. We were unable to align the intron regions because of extensive divergence across the taxa sampled. However, these regions can potentially be useful for more focused phylogenetic studies, and it should be possible to use information on exon/intron boundaries in our supercontig files to design more intron-aware baits for more phylogenetically targeted studies. We were able to get approximately 1 million reads per sample for the 25 genes by loading 24 samples on a single Illumina MiSeq lane. This appears In comparison with the phylogenetic results of Rothfels et al. (2015b), we find the same broad relationships, including Equisetum L. sister to the remainder of ferns, Marattiales sister to leptosporangiate ferns, Cyatheales sister to Polypodiales, and Dennstaedtiaceae sister to the eupolypods (for the latter three relationships, we find increased maximum likelihood bootstrap support: 100% vs. 79%, 100% vs. 57%, and 100% vs. 89%, respectively). In addition, we resolve (but do not support) a monophyletic Gleicheniales and we find marginal support (72%) for a Gleicheniales + Hymenophyllales clade; neither relationship was resolved by Rothfels et al. (2015b) but has been reported in other previous studies (Pryer et al., 2004;Shen et al., 2017). The novel taxa sampled add further insights: Dicksonia L'Hér. and Metaxya, representing two previously unsampled families, are well supported in a clade with Cyatheaceae, although the relationships among the three families are uncertain; Cystodiaceae (Cystodium J. Sm.) is resolved as sister to Lindsaeaceae with 86%  FIGURE 3. Maximum likelihood phylogeny from our sequence-capture data set combined with the transcriptome-derived data from Rothfels et al. (2015b). Thickened branches indicate bootstrap support greater than 80%, and bold taxon names followed by an asterisk indicate accessions from the sequence-capture data set (see Fig. 1 (Korall et al., 2006a(Korall et al., , 2006bLehtonen et al., 2012;PPG, 2016). Taxon sampling from Dennstaedtiaceae and Lindsaeaceae remains sparse and more work is required in this area to robustly resolve the early evolutionary history of the Polypodiales. Nevertheless, results from phylogenetic analyses suggest that our approach is capable of capturing DNA sequences with phylogenetic signal. The genes we chose were specifically selected because they have phylogenetic utility across ferns. Not all of these genes may necessarily be useful at shallow levels, e.g., comparing closely related species. Furthermore, we have not specifically tested how effective these genes are at detecting hybrids, because we did not include any suspected hybrids in the study. However, by choosing genes for which nucleotide sequences are already available, researchers will more likely be able to detect hybrids in future studies using a growing collection of available data. The main challenge in analyzing nuclear gene data for analysis of hybrids is phasing of alleles and homeologs, especially in young hybrids. This is an issue with all such data and can be minimized by favoring long-read sequencing of MiSeq PE250 over the shorter reads from HiSeq. It should also be possible to combine the target sequence capture approach with sequencing on the PacBio platform for even longer reads, which has been shown to be effective for resolving polyploid and hybrid origins (Rothfels et al., 2017). As evolutionary studies increasingly switch to nuclear gene data, our bait collection should provide a useful tool for a range of projects across ferns. Design of baits for additional genomic regions should soon be possible using accumulating transcriptome data  as well as sequences from two fern genome projects (Azolla and Salvinia Genomes Consortium, unpublished data). Our approach to gathering data for ferns is complementary to other projects in progress at a broader scale, including the Genealogy of Flagellate Plants (GoFlag) project (http://flag ellateplants.group.ufl.edu/) and the Plant and Fungal Trees of Life (PAFTOL) project (https://www.kew.org/science/who-we-areand-what-we-do/strategic-outputs-2020/plant-and-fungal-treeslife). Our hope is that the accumulation of nuclear gene data will enhance our understanding of fern phylogeny and contribute to an increased understanding of the role of hybridization to fern evolution, across both deep and shallow phylogenetic depths. Although we are not yet able to evaluate the level of overlap between the regions captured in the method described here and those captured in the broader-scale studies mentioned previously, we believe there is value in data sets that differ, as well as in those that have sufficient overlap to be combined. Just as phylogenies generated with sequence data from different regions of chloroplast DNA in ferns provided us with a better understanding of the general pattern of fern evolution (Hasebe et al., 1995;Schuettpelz and Pryer, 2007;Kuo et al., 2011) and identified challenging taxa (Pryer et al., 2004;Rothfels et al., 2012;Testo and Sundue, 2016;Wei et al., 2017), those generated from differing subsets of the nuclear genome should further elucidate the evolutionary history of the group. As a diversity of methods become available over the coming years, researchers will be able to choose an approach that is best suited for their needs.