The genus Yucca (L) is a group of ~50 species of woody monocots endemic to the North American arid regions. Their obligate pollination mutualism with yucca moths is considered a 'textbook example' of coevolution, and is hypothesized to have promoted rapid diversification. However, testing this hypothesis has been difficult due to uncertainty about the placement of a rogue taxon, Yucca queretaroensis, a rare endemic of the Sierra Gorda region of central Mexico. Past work placed this species in different positions within the Agavoideae, producing starkly different age estimates for Yucca, from 25 MY to 4 MY.

We generated new sequence capture data for 353 nuclear genes and all coding regions of the plastid genome from wild-collected plants and samples included in previous studies, producing a new phylogeny and a new age estimate for Yucca.

Here, we have archived sequence data and morphological measurements associated with a phylogenomic study of Yucca. Includes gene sequence assemblies from the Angiosperm 353 probe set, gene sequence assemblies from the coding regions of the plastid genome, Sanger sequence data from the chloroplast trnL-F region, Sanger sequence data from the insect mitochondrial cytochrome oxidase one gene, partial outpus from the Hybphaser program used to estimate paralogy and hybridization from Angiosperm 353 data, and a table of morphological measurements. 

Collections: To generate new genetic data for Y. queretaroensis, we obtained leaf tissue from 33 individuals of that species. This number includes two samples obtained from the Greg Starr Nursery (a commercial horticultural collection in Tucson, Arizona), two samples from the Pellmyr tissue collection used in past phylogenetic studies (Pellmyr 146q and Pellmyr 311), four from the Botanical Garden at the National Autonomous University in Mexico City (hereafter Jardin Botánico, UNAM), and 25 samples collected in the wild. In addition, we obtained leaf tissue from five known or suspected hybrids between Y. filifera and Y. queretaroensis - two from the Gregg Starr Nursery, one from the Cadereyta Regional Botanical Garden in Querétaro, Mexico, and two identified in the wild.

We identified wild populations of Y. queretaroensis using existing collection records provided by the Cadereyta Regional Botanical Garden and crowd-sourced, georeferenced observations from iNaturalist, and using species distribution models contained in the CITES petition (Magallán-Hernández, 2013) to predict where potential new localities might exist. In addition, we met with local municipal officials, interviewed local landowners, and worked with private guides to identify previously undocumented populations. Finally, three new localities were identified by chance while driving between existing sites.

From each plant sampled in the wild we measured the length and width of one leaf, with width measured at the widest point. We recorded the location of each plant using a Garmin eTrex GPS unit (Garmin International, Olathe, KS) and photographed each plant including a standard measuring pole for scale. At each sampled site, we collected leaves from a single plant to serve as vouchers. We preserved wild-collected leaf tissue by drying in a standard plant press, and by placing sections of leaf cut into ~5 cm by ~1 cm strips in silica gel.

We supplemented these collections with 55 additional accessions from other species within the Agavoideae. These include leaf tissue collected in the wild and in botanical gardens, tissue from existing collections provided by other investigators, and by Gregg Starr Nursery and Cistus Gardens (Portland, Oregon). These materials were preserved in several different ways. Leaf tissues from the Pellmyr collection were initially maintained in a cooler in the field before being transferred to a -80 ℃ freezer upon return to the lab and were maintained at this temperature for decades. Materials from the Clara-Arteaga collection and those sampled from horticultural gardens were dried on silica. The full set of tissue samples included 84 individuals representing 25 species of Yucca, two individuals from different species of Agave, and two samples of D. longissiumum.

Moth collections: There are no existing collections of adult yucca moth pollinators of Y. queretaroensis. Instead, we obtained three larvae extracted from fruits of an individual Y. queretaroensis growing in the Cadereyta Regional Botanical Garden. The botanical garden at Cadereyta is outside of the natural range of Y. queretaroensis, but one individual had bloomed coincident with local Y. filifera plants and subsequently set fruit. We presume that this plant was cross-fertilized by Y. filifera, which is pollinated by Tegeticula tambasi (Pellmyr et al., 2008). Seeds and larvae from these fruits were extracted. The larvae were preserved in ethanol and then dried.

Sanger sequencing: To evaluate whether misidentification or laboratory error may have resulted in an incorrect placement of Y. queretaroensis in previous studies, we generated new Sanger sequence data for the plastid tRNA leucine (trnL) and the trnL-trnF intergenic spacer.

We amplified the trnL-trnF region by PCR using the trnL forward primer (CGAAATCGGTAGACGCTACG), and the trnL-F reverse primer (ATTTGAACTGGTGACACGAG) (Taberlet et al., 1991). We performed PCR reactions using the Apex Bioresearch 2.0X Taq RED Master Mix Kit (APExBIO Technology LLC , Houston, Texas) in 25 uL volumes including 1uL of each primer, 9.5 ul of reagent grade water, and 2ul of genomic DNA per reaction. The thermocycler program began with an initial denature at 95 ℃ for 60s, followed by 40 cycles of the following: a 30 s incubation at 95 ℃, then a 60 s incubation at 47 ℃, followed by 90 s incubation at 72 ℃. The program ended with a final long extension step, incubating at 72 ℃ for 3 mins. We prepared samples for sequencing using an ExoCIP Rapid PCR Cleanup Kit (New England Biolabs, Ipswich, Massachusetts), and they were sequenced in forward (5’ to 3’) and reverse (3’ to 5’) directions by Eurofins Genomics (Louisville, Kentucky) using an ABI 3700 capillary sequencer.

To determine the identity of insect larvae extracted from fruits of Y. queretaroensis, we amplified the mitochondrial Cytochrome Oxidase One (COI) gene using the primers S2183 (5′-CAACATTTATTTTGATTTTTTGG-3′) and A3020 (5′-TCCAATGCACTAATCTGCCATATTA-3′) (Hebert et al., 2016). Each 25 μL reaction contained approximately 10–100 ng of genomic DNA, 1X PCR buffer, 1.5 mM MgCl₂, 0.2 mM of each dNTP, 0.4 μM of each primer, and 0.5 U of Taq DNA polymerase (Promega Corporation, Madison, WI). Amplifications were carried out under the following thermal cycling conditions: an initial denaturation at 94 °C for 5 min; followed by 35 cycles of 94 °C for 30 s, 48 °C for 30 s, and 72 °C for 1 min; with a final extension at 72 °C for 10 min. The PCR reactions were cleaned and then sequenced by Eurofins Genomics (Louisville, Kentucky) using an ABI 3700 capillary sequencer.

Seq capture / Shotgun seq: We generated Illumina sequence data in two batches: an initial pilot data collection batch including only 32 samples, and a second batch containing all 88 samples. In the first batch, we fragmented the whole genomic DNAs prior to library prep using a Covaris Focused-ultrasonicator M220 (Covaris, LLC., Woburn, Massachusetts) and then prepared them for sequencing using the NEBNext Ultra II DNA Library Prep kit for Illumina (New England BioLabs, Ipswich, Massachussetts). In the second batch, we prepared libraries using the NEBNext Ultra II FS DNA Library Prep kit for Illumina (New England BioLabs, Ipswich, Massachussetts), which incorporates an enzymatic fragmentation step as part of the protocol. In both cases, we added sample-specific indices using the NEBNext® Multiplex Oligos for Illumina Dual Index Primers Sets 1 and 3 (New England BioLabs, Ipswich, Massachussetts). We enriched our libraries for target genes using the Arbor Biosciences myBaits® Angiosperms 353 v1 Target Sequence Capture Kit (Arbor Biosciences LLC, Ann Arbor, Michigan), following the manufacturer’s protocol.

In the initial sequencing effort, we sequenced only enriched libraries. The second batch contained enriched libraries for 56 new samples, as well as unenriched whole-genomic libraries for all 88 samples. We included unenriched libraries to increase the number of reads from the plastid genome, which predominates in whole genomic DNA from leaf tissue. All libraries were sequenced on an Illumina NOVA Seq 6000 platform at the University of Oregon Genomics and Cell Characterization Core Facility, each using paired-end 150 base pair reads, but with differing data output. The initial pilot sequencing effort used a NovaSeq 6000: SP PE 150 nt run; the second batch used a NovaSeq 6000: S1 PE 150 nt run. FASTQ data were demultiplexed by the core facility based on barcode sequences we provided.

Analysis Illumina sequence data: We analyzed the initial FASTQ files using FastQC version 0.12.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). FastQC identified two over-represented sequences, which matched the NEBNext adapter sequences. We then cleaned our sequences with Trimmomatic version 0.39, (Bolger et al., 2014) using a custom file containing the NEBNext adapter sequences. In the process, we renamed our cleaned output files to include individual accessions and sample IDs in the filename. After cleaning, the data passed all FastQC tests, but with warnings (per tile sequence quality, sequence length distribution, and sequence duplication levels). We assembled our cleaned FASTQ data into gene sequence alignments using the Hybpiper v 2.17 pipeline (Johnson et al., 2016).

Assembly of nuclear (Angiosperm 353) loci and phylogenetic inference: For the Angiosperm 353 genes, we constructed a custom target file using sequences from species within the Agavoideae obtained from Kew Tree of Life Explorer (Baker et al., 2022), including Agave fourcroydes Lem., Agave tequilana F.A.C.Weber, Beschorneria yuccoides K. Koch, Camassia quamash (Pursh) Greene, Hesperaloe nocturna Gentry, Hesperaloe parviflora (Torr.) J.M.Coult, Yucca brevifolia Engelm., and Yucca filamentosa Riddell. We aligned the cleaned FASTQ data against this target file using the Hybpiper ‘assemble’ pipeline. We extracted exon sequences, introns, and ‘supercontigs’ (containing both exons and introns) using the ‘retrieve_sequences’ pipeline.

Of the 353 target gene regions, 341 were recovered in at least 90% of samples (see Appendix S7 for details of gene recovery and gene lengths). To test for the presence of paralogous gene copies and to evaluate their impact on phylogenetic inference, we used the ‘paralog retriever’ pipeline in Hybpiper2 to identify putative paralogous gene copies in each sample. Paralog retriever identified 91 genes in which at least one sample showed evidence of multiple contigs or higher than expected coverage. To reduce their potential impact on phylogenetic estimation, we constructed ASTRAL trees including only the 262 genes that showed no evidence of paralogs. We aligned the resulting FASTA files for these genes using MAFFT (Katoh et al., 2002), and converted the alignments to phylip format using PGDSpider version 2.1.1.5 (Lischer and Excoffier, 2011) and a custom shell script (provided in Appendix S8) to insert a space between the end of each taxon name and the start of sequence data. We estimated gene trees with RAXML version 8.2.12 (Rokas, 2011) from exon sequences, introns, and supercontigs using a general time reversible model with gamma distributed rates. We concatenated the best trees generated for each alignment into one of three tree files, representing gene trees based on exons only, introns only, or supercontigs. We constructed trees using Astral v. 5.7.7 (Zhang et al., 2018) using the default settings.

Reconstruction of plastid genes and phylogenetic analysis: To reconstruct coding regions from the plastid genome we again used Hybpiper2 to align Illumina sequence reads against sequence targets. Unlike the nuclear Angiosperm 353 data, plastid genes were not enriched in the sequencing libraries using sequence capture probes; because the plastid genome is much more abundant than the nuclear genome, unenriched whole genomic shotgun sequence data typically provides adequate depth of coverage to assemble plastid genes. We created a custom target file containing sequences of 102 plastid genes from five species of Yucca, originally described in McKain (2016) (Genbank accessions KX931469.1, KX931468, MW281827, KX931467.1, and KX931466). As with the Angiosperm 353 data, we aligned the cleaned FASTQ data against this target file using the ‘assemble’ pipeline, and gene sequences were each extracted using the ‘retrieve_sequences’ pipeline. Introns and supercontigs were not retrieved. We edited the resulting FASTA files to remove information describing the stitching of contigs using a custom shell script (see Appendix S8) and aligned sequences from each gene against outgroup sequences region using MAFFT. We selected outgroups following Smith et al. (2021), with a goal of providing multiple age calibration points within Asparagales. For a complete list of outgroups and GenBank Accessions see Appendix S9. We concatenated the plastid gene alignments to create a single nexus file using phyutility v. 2.7.1 (Smith and Dunn, 2008).

We estimated phylogenetic relationships and divergence times between plastid gene sequences using BEAST v. 10.5.0-beta2 (Drummond and Bouckaert, 2015). We estimated rates of evolution assuming an uncorrelated relaxed clock model with a log-normal distribution using a mean of 1.0 and standard deviation of 1.0. The ages of key nodes were constrained based on fossil calibrations described in Iles et al. (2015), as implemented in Smith et al. (2021). All age constraints used a lognormal prior distribution, with shape parameters and offsets as follows: we set the age of the common ancestor of Yucca and Agave to 14 MY with a mean of 1.0 and a standard deviation of 1.5; we set the age of  Arecacace to 83 MY with a mean of 8.0 and a standard deviation of 1.0; we set the age of  Goodyearinae to 15 MY with a mean of 1.5 and a standard deviation of 1.0; we set the age of  Zingiberales to 72 MY with a mean of 7.0 and a standard deviation of 1.0; we set the age of the root to 113 my with mean of 5.0 and a standard deviation of 4.0. The model of sequence evolution assumed a general time reversible model with gamma-distributed rates using default priors. The analysis used a 30 million generation Markov chain, discarding the first 15 million generations as burn-in. We estimated the maximum clade credibility tree using TreeAnnotator v. 10,5.0-beta2 (a software distributed as part of the BEAST package), calculating the median node height. The posterior distribution and maximum a posteriori estimates of the ages of the common ancestors of key groups were calculated using Tracer v. 1.7.2 (Rambaut et al., 2018)

Tests of hybridity: To evaluate evidence of hybridization, we used Hybphaser v 2.0 (Nauheimer et al., 2021) to calculate the average heterozygosity across all loci, and mean sequence divergence between alleles at each locus, for each individual. We then compared these values between four groups of samples: putative hybrid individuals of Y. queretaroensis (identified based on vegetative morphology), all individuals of Y. queretaroensis including putative hybrids, all Yucca samples, and the complete dataset.