Agave utahensis Engelm. (Utah agave; Asparagaceae) and its putative subspecies, A. utahensis subsp. kaibabensis (McKelvey) Gentry and A. utahensis subsp. utahensis, play critical roles as keystone taxa in their native habitats. Utah agave contributes significantly to soil formation and the creation of oases of organic matter in the Mojave Desert and Colorado Plateau in the southwestern United States (Gentry, 1982). Their large inflorescences produce copious amounts of nectar for a wide array of pollinators, including bats, birds, bees, and hawkmoths (Slauson, 2001). The rise in frequent and intense fires across the distribution of Utah agave is threatening native populations and consequently the genetic diversity of the species. As with other agaves, Utah agave reproduces clonally via vegetative rhizomes (Gentry, 1982; Arizaga and Ezcurra, 2002). As such, the risk exists that its asexual reproductive strategy could hasten its ecological decline if fire frequencies continue to increase.

Moreover, taxonomic questions exist concerning the subspecies ranking within this highly variable species complex. The subspecies ranking of the two Utah agave (subsp. kaibabensis and subsp. utahensis) is based solely on morphological traits, with subsp. kaibabensis having wider and greener leaves (3.0–5.5 cm) and larger rosettes (40–100 cm) than subsp. utahensis (1.0–3.0 cm, 15–40 cm, respectively) (Gentry, 1982; Reveal and Hodgson, 2003; Janeba, 2008). DNA marker information will help characterize the genetic diversity of the Utah agave species complex and confirm whether subspecies ranking is justified.

Despite its importance as a keystone species, no molecular tools have been developed to manage Utah agave or to investigate population structure, subspecies ranking, hybridization, or polyploidy within the species. We constructed primers to amplify microsatellite markers from both subspecies along a west-east transect across southwestern Utah and north-central Arizona.

METHODS AND RESULTS

The genomic reduction protocol used here is described in detail by Maughan et al. (2009) for single-nucleotide polymorphism discovery. Here we use the same methodology, but analyze the results for putatively polymorphic microsatellite sequences. The method is based on the conservation of restriction sites across individuals followed by genomic selection and en masse fragment sequencing via next-generation sequencing technology. Specific DNA barcodes are incorporated into the DNA fragment of each member of the screening panel, which are subsequently used to deconvolute the sequence read pool and to search for polymorphic microsatellite loci across the assembled DNA fragments. Leaf tissue was lyophilized and DNA extracted using a modified mini-salts protocol as reported by Todd and Vodkin (1996). In short, 450 ng of DNA from two A. utahensis subsp. kaibabensis and two A. utahensis subsp. utahensis agave accessions were separately double-digested for 1 h at 37°C with 3 units of the restriction enzymes EcoRI and BfaI (New England Biolabs, Beverly, Massachusetts, USA) in 1× NEB4 restriction buffer. The resultant DNA fragments were immediately ligated with 1.5 µM 5′-TEG biotinylated/3′-phosphorylated EcoRI adapters and 15 µM 3′-phosphorylated BfaI adapters (see Maughan et al., 2010 for adapter sequences) using 3 units of T₄ ligase (New England Biolabs) at 16°C for 3 h. Non-biotin-labeled DNA fragments (BfaI to BfaI restriction fragments) were removed from the reaction using a biotin-streptavidin paramagnetic bead separation using M-280 streptavidin beads (Invitrogen, Carlsbad, California, USA) according to the manufacturer's specifications. The remaining EcoRI-BfaI and EcoRI-EcoRI fragments, with the biotin label still attached, were resuspended in 100 µL of 10 : 1 TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0).

Four primer pairs, designed to be complementary to the EcoRI and BfaI adapter sequences and to carry unique 5′ 10-base barcode sequences (Maughan et al., 2010), were synthesized by Integrated DNA Technologies (Iowa City, Iowa, USA). A single barcode primer pair was used to amplify 1 µL of the streptavidin-cleaned DNA fragments for each DNA sample. Amplification of each sample was performed in 50-µL PCR reactions using 1× Advantage HF 2 PCR Master Mix (ClonTech, Mountain View, California, USA) and 0.2 µM of the MIDX-EcoRI and MIDX-BfaI primer pairs, with the following thermocycling profile: 95°C for 1 min followed by 22 cycles of 95°C for 15 s, 65°C for 30 s, and 68°C for 2 min. Equimolar amounts of DNA for each of the PCR reactions were pooled and sequenced using standard protocols for 454-pyrosequencing as a service at the Brigham Young University DNA Sequencing Center (DNASC, Provo, Utah, USA) using a Roche-454 GS FLX instrument and Titanium reagents (454 Life Sciences, a Roche Company, Branford, Connecticut, USA) without DNA fragmentation. DNA reads were bioinformatically trimmed and separated into barcode pools representing the four Utah agave accessions and then de novo assembled together using the default parameters of Roche Newbler assembler (version 2.3; 454 Life Sciences, a Roche Company). All 786,104 sequence reads were assembled into 22,386 large reference contigs (>400 bp). Using the computer program MISA (Thiel et al., 2003), the reference contigs were scanned for perfect di-, tri-, and tetranucleotide motifs with minimum repeat units of eight, six, and five, respectively .

A total of 1013 putative microsatellite markers were identified, including 762 di-, 205 tri-, and 46 tetranucleotide motifs. The most common motifs were AG/CT, AAG/CTT, and AAAC/GTTT. The MISA output files and the reference contig.ace file were then processed using custom Perl scripts (Stajich et al., 2002) to preselect microsatellite loci that varied in repeat length as determined from consensus reads for each of the agave subspecies (utahensis vs. kaibabensis). A total of 298 putatively polymorphic microsatellites were identified, of which 15 trinucleotides representing the range of motifs (AAT, ACA, CAC, CCA, GAG, GCT, GTT, TAT, TTC) and repeat lengths (5–10 repeats) were selected for primer design using Primer3 version 2.0 software (Rozen and Skaletsky, 2000) and validation on a diversity panel of 104 Utah agave accessions collected from eight populations in southern Utah and northern Arizona (Appendix 1). Amplification of each sample was performed using multiplex PCR reactions, which included three primers, each with different fluorescent tags (VIC, 6-FAM, NED; Table 1). Multiplex PCR reactions were performed using a Type-it Microsatellite PCR Kit (QIAGEN, Germantown, Maryland, USA) in 10-µL reactions with 30 ng of genomic DNA and 0.2 µM of all primers according to the manufacturer guidelines with the following amplification parameters: 95°C for 5 min, followed by 28 cycles of 95°C for 30 s, 57°C for 90 s, 72°C for 30 s, followed by a final extension cycle of 60°C for 30 min. The amplified products were diluted 1 : 60, and 1 µL of the dried products were analyzed using fragment analysis on an ABI 3730xl using GeneScan 500 ROX (Applied Biosystems, Carlsbad, California, USA) as the size standard at the Brigham Young University DNASC. The data generated from the fragment analysis were scored using the microsatellite plugin in Geneious R6 (Drummond et al., 2011) and analyzed using default parameters of the computer program Arlequin version 3.5 (Excoffier and Lischer, 2010).

Four of the 15 preselected markers failed to amplify or produced poor, weak, or ambiguous amplification products. Eleven markers were polymorphic, producing a total of 61 alleles across both subspecies. The number of alleles per locus varied from a low of three alleles (BYU5866) to a high of eight alleles (BYU5164 and BYU8677), with an average of 5.5 alleles per marker (Table 2). Interestingly, two of the microsatellite markers (BYU4463 and BYU4988) exhibited individuals (∼25%) that amplified between two and four alleles, presumably due to the amplification of two paralogous DNA sequences. The mean observed and expected heterozygosity values for A. utahensis subsp. utahensis and A. utahensis subsp. kaibabensis were 0.357 and 0.491 and 0.357 and 0.408, respectively (Table 2). Division of the diversity panel by putative subspecies ranking exhibited a between-subspecies genic differentiation (F_ST) of 0.24. When all individuals were analyzed for Hardy–Weinberg equilibrium (HWE) on a subspecies basis, irrespective of subpopulation, several markers violated HWE (Table 2). This was likely due to population substructure (i.e., Wahlund effect). Indeed, when the data were analyzed for each of the eight subpopulations (Appendix 1), all microsatellites were in HWE for five of the populations, with only BYU4012 or BYU8677 violating HWE in the remaining three populations (note that the sample sizes for these populations are too small for definitive conclusions).

Table 1.

Characteristics of 11 polymorphic microsatellite primers developed for Agave utahensis.

Table 2.

Genetic properties, including number of alleles detected at each microsatellite marker and observed and expected heterozygosity values, for each subspecies.^a

CONCLUSIONS

We used a novel genomic reduction approach tied to 454-pyrosequencing and DNA sample barcoding to identify polymorphic marker loci in A. utahensis. Eleven microsatellite markers were validated on 104 accessions of A. utahensis subsp. utahensis and subsp. kaibabensis. These markers should prove invaluable for performing larger investigations of population structure, genetic diversity, and phylogenetic relationships in agaves.