Taxon Sampling

For sampling of the Del Mar manzanita, we visited 14 locations in California and Baja California, Mexico (Table 1, Fig. 1, Appendix 1). In some large populations, we sampled from several subpopulations. To prevent confusion between populations and subpopulations, we refer to plant collecting locations as locales, which we define as groups of plants that are discontinuous from one another, with centers at least one km apart. In every case, the edges of the locales were at least 0.5 km apart. We were not able to obtain estimates of population size. However, the spatial density and size of populations/subpopulations vary dramatically across the range of the Del Mar manzanita (AECOM 2015), with consequences for population genetics.

Table 1.

Collecting Locales and Sampling Information. Taxon represents A. glandulosa individuals identified according to dichotomous keys in the Jepson Manual (Parker et al. 2012). Code indicates a field collection code for each locale used as a short-hand to refer to individual locales; each code corresponds to an herbarium voucher specimen housed at the UC Davis Center for Plant Diversity or the San Diego Natural History Museum (only one herbarium voucher was collected per site). M. Mulligan sampled locales 3318 and 3321; D. Burge sampled all other locales. Samples correspond to the number of individuals sequenced. Latitude and longitude are given in the WGS 84 datum.

Fig. 1. 

Sampling map. Grid dimensions are in degrees of latitude and longitude (WGS84 datum). See Table 1 for more information on locales.

We sampled intensively in central San Diego County, especially at the Marine Corps Air Station (MCAS) Miramar, a federally owned property managed by the United States Department of Defense (Fig. 1). MCAS Miramar is thought to contain large populations of Del Mar manzanita, and is also a zone of potential contact between this taxon and Eastwood manzanita (Rebman and Dossey 2002).

To aid in determining the correspondence of genetics with morphometric taxon concepts, we sampled from the taxon thought to be most closely related to Del Mar manzanita, the Eastwood manzanita (six locales, Table 1, Fig. 1), and from Arctostaphylos glandulosa subsp. atumescens J.E.Keeley, M.C.Vasey & V.T.Parker, a Mexican endemic that is the only subspecies of A. glandulosa other than the Del Mar manzanita and the Eastwood manzanita known to occur close to the coast in southern California and Baja California, Mexico (Table 1, Fig. 1). This latter species is a narrow endemic found only on Punta Banda Peninsula in Baja California, Mexico (Fig. 1), where it is the only manzanita. Thus, while this species does occur in the same region as Del Mar manzanita, they are not known or expected to come into close contact.

Here we use scientific names that are based on the morphology of the plants prior to DNA sequencing. Because the purpose of our study was to find out whether plants identified as Del Mar manzanita using morphological criteria could be reliably grouped together based on genetic similarity, we began by using names based on the morphology of the plants that we sampled (Table 1). Plants were identified according to the taxonomic keys and descriptions in the Jepson Manual (Parker et al. 2012), supplemented by relevant taxonomic revisions (Keeley et al. 2007). In both Parker et al. (2012) and Keeley et al. (2007), the presence of glandular hairs on leaves, twigs, and infloresences, as well as color of the leaves and shape of the fruit are used to discern Del Mar manzanita from other taxa. Using this means of identification allowed us to begin with a hypothesis as to the correct identification of the plants, which was then tested using DNA. In the results presented below, we call attention to groups of plants in which the DNA results do not agree with the morphology-based names.

We took special care to sample multiple plants from the type localities of Del Mar manzanita (Del Mar, CA; Jepson 1922) and Eastwood manzanita (Mount Tamalpais, CA; Eastwood 1897). We included the type localities in the analysis so that we could make progress toward establishing the genetic limits of these taxa. We recognize that this is not ideal in the eyes of some, and that it might be preferable to sample DNA directly from the type specimens for such work. However, preliminary tests showed that DNA is not recoverable from older herbarium specimens in sufficient quantities to allow RAD DNA sequencing, the genomic method employed here. Therefore, sampling living plants at the type localities, rather than from the type specimens themselves, represents a next-best approach. This may also have the advantage of identifying relevant variation in morphology and genetics currently present in the population of the type locality.

For all sampled plants, we collected young leaves and flower buds, which were frozen on dry ice within 48 hours of collection. Up to 10 individuals were sampled from each locale, but not all collected individuals were sequenced (Table 1). We collected more individuals than required due to the low rate of successful DNA isolation; only high-quality samples were used.

Morphometrics

Morphologically, Del Mar manzanita is distinguished from other subspecies of A. glandulosa by qualitative characteristics of the reproductive and vegetative parts (Parker et al. 2012; Keeley et al. 2007; Table 2), many of which are subjective, particularly in the case of leaf color. Because most of the sites were visited in the spring, to obtain fresh flower buds for DNA extraction (see below), we were unable to collect material that included fruit and were therefore unable to use fruit characters in our morphometric analysis. We also did not collect data on leaf color. Though leaf color is often used to help distinguish subspecies of A. glandulosa from one another (Table 2), we found that it was difficult to assess this character in our dried herbarium specimens, due to differences in drying conditions. Instead, we focused on the hairs of the leaves and youngest stems (Table 2). If present at all, these hairs are typically in two layers, one of longer, stiff hairs that sometimes are tipped with glands, and another of shorter hairs that are typically appressed to the stem below the level of the long hairs.

Table 2.

Key Morphological Traits Traditionally Used to Distinguish Subspecies of A. glandulosa. Traits are taken from Keeley et al. (2007) and Parker et al. (2012).

For each of the specimens from which we obtained genetic data, we also collected data on the following traits: (A) long hair density, the density of long hairs on the young stem (0, absent; 1, sparse; 2, intermediate; 3, very dense), (B) long hair type, the texture of the long hairs (1, stiff; 2, soft), (C) long hair glands, the presence of glands at the tips of the long hairs (0, absent; 1, present), (D) short hair presence, the presence of short hairs below the layer of the long hairs (0, absent; 1, present), (E) short hair type, the type of the short hairs (1, stiff and erect; 2, soft and appressed), (F) margin hairs, the presence or absence of hairs on the margin of the leaf, (0, absent; 1, present), (G) margin hair type, type of the margin hairs, if present (1, long, glandular; 2, short, non-glandular; 3, stiff, long, nonglandular), (H) quantitative leaf shape (leaf length (mm) divided by width (mm); average for three mature leaves). Raw morphological data were deposited at the Dryad Digital Repository (Appendix S1; doi:  10.5061/dryad.kv573c7).

The vegetative morphometric data were treated in a multivariate framework. We visualized data using principal components analysis (PCA) and tree reconstruction, carried out in R, version 3.1.2 (R Foundation for Statistical Computing, Vienna, Austria). For PCA, we used the prcomp function in R to model the seven categorical variables and the single quantitative variable (H, leaf shape). Quantitative leaf shape was transformed into a Z-score before analysis. Also, due to the needs of the model, all individual plants with “missing” data (cases in which characters were not applicable to the individual in question) were excluded from analysis. The first two principal components were visualized in bivariate space to examine relationships. The contribution of each morphological character to the principal components was determined based on vector loadings. Based on preliminary results of the PCA analysis, it was determined that only A, C, D, and G contributed significantly to the model (absolute value of maximum vector loading on first and second principal component < 0.09).

Based on the seven categorical variables, we constructed a classification tree using the rpart R package, version 4.1-10. The rpart package implements recursive partitioning (Therneau and Atkinson 1997; De'ath 2002). In recursive partitioning, a tree linking all of the observations (individual plants in this case) is built by a simple process: a categorical variable is selected that “best” splits the individuals into two groups. This process is applied until the subgroups either reach a minimum size, or all individuals are grouped. We set the “minsplit” parameter of rpart to 2, the “minbucket” parameter of rpart to 1, and the “cp” parameter to 0.001. We did not prune the resulting tree, because cross-validation suggested that no improvement could be made.

In addition to recursive partitioning, we constructed a distance-based tree to visualize relationships among individual plants based on their overall dissimilarity across the categorical variables that we scored. We constructed a UPGMA (Unweighted Pair Group Method with Arithmetic Mean) tree using the vegan package of R, version 2.3-4. We used the hclust function of vegan, with ‘method' set to average. Based on the preliminary results of the PCA, we excluded traits that did not contribute strongly to the preliminary PCA analysis, as these would likely add noise to the UPGMA tree; the UPGMA tree used traits A, C, D and G (see above).

DNA Sequencing & Variant Detection

Total genomic DNA was extracted from flower buds using the DNeasy Plant Mini Kit (Qiagen, Germantown, MD) according to the manufacturer's instructions. Total genomic DNA was checked for degradation on a 1% agarose electrophoresis gel; samples with high molecular-weight DNA (∼20-50 kb) were standardized to 30 ng/μL (diluted to a volume of 150 μL in TE) and sent for quality control, library construction, and sequencing at Floragenex, Inc. (Portland, OR). Methods for this process generally followed those described by Lozier (2014), with the exception that our work used the restriction enzyme PstI (Hipp et al. 2014). In brief, sequence identifier barcodes (a unique one for each plant) and sequence adapters were added to genomic DNA after digestion by the endonuclease PstI. The barcoded samples were then combined and the fragments sequenced outwards from restriction-sites using single-end DNA sequencing reads that were 100 base pairs long. All DNA sequencing was done on a HiSeq 2500 DNA sequencing instrument (Illumina, San Diego, CA). Following sequencing, DNA sequences from individual plants were separated based on their unique barcodes using the program fastq-multx (Aronesty 2013); the barcodes were then removed from the sequences.

Due to the lack of an existing reference genome for Arctostaphylos or a closely related member of the Ericaceae, detection of variation among individual plants was done using a “RAD reference” approach (Lozier 2014), in which a kind of reference genome is constructed using sequence data from the plant with the greatest number of unique RAD clusters. We employed this strategy to ensure the best reference genome “target” for subsequent alignment of sequences from other plants. In developing the RAD reference, custom methods developed by Floragenex (Lozier 2014) were used to cluster identical sequences that had 5–500× sequencing coverage, which produced a preliminary assembly. The assembly was then collapsed back to separate sequences and these sequences were realigned against the preliminary RAD reference genome using the program BWA [Li and Durbin 2009; aln function, edit distance (–n) 3, –N (disable iterative search)], allowing at most four mismatches among reads within a cluster. The purpose of self-alignment was to identify and remove repetitive DNA regions.

Reads for each individual plant were aligned to the RAD reference using BOWTIE (Langmead et al. 2009), relying on sequence quality information to aid in the process of match-making, allowing a maximum of three mismatched bases per read, and permitting alignment of each read to no more than one region of the reference. SAMTOOLS (Li et al. 2009) was used to detect SNPs and call genotypes. Filters for SNP calling required a minimum phred score of 20, a minimum of 15× sequence coverage and a maximum of 10% missing data across samples.

Following SNP calling, we also excluded: 1) all positions in which data were missing for any individual plant (Arnold et al. 2013; Davey et al. 2013), 2) all invariant sites (where every individual had the same allele call), and 3) all sites in which more than two allelic variants were detected. The final filtering step restricted the dataset to only biallelic SNPs, thus rendering the data diploid. We used this filter because most population genetic software are not able to deal with polyploid data (where more than two allelic variants are allowed per locus). By rendering our data as biallelic SNPs, we are able to apply population genetic methods that assume diploid loci. We recognize that this is not ideal, as the underlying genetics of the organisms is tetraploid. We think that this filtering method is the best option given the limitations of available SNP calling methods. We also note that this method of dealing with polyploids in population genetic analyses has been employed in other studies (Qi et al. 2015; reviewed by Dufresne et al. 2014). The rendering of loci as diploid when the underlying genetics is polyploid can lead to conflicts with assumptions of software, mainly due to the issue of calculating allele dosage in polyploids (Dufresne et al. 2014). Methods are available to deal with polyploids in genomic data by calling alleles using a draft reference genome (Garrison and Marth 2012). Unfortunately, the lack of a manzanita reference genome makes it impossible to apply such methods at this time.

The data resulting from the above filtering processes was used to create two datasets, one using all of the A. glandulosa samples (54 individuals; Table 1), and the second focused on just the locales classified as Del Mar manzanita, and for which we had more than one sample (39 individuals; Table 1). Hereafter, we refer to the former as the “Complete” and to the latter as the “Reduced” dataset. In both datasets, we excluded invariant sites, as well as sites with a minor allele frequency of less than 0.02 (Bradbury et al. 2007).

Population Genetics

To objectively identify statistically meaningful groups of individuals, we used a combination of parametric and non-parametric approaches. These approaches were applied in exactly the same way to all datasets. For the non-parametric approach, we carried out principal coordinates analysis (PCoA) using an identity by site genetic distance matrix. For the parametric approach, we used the program STRUCTURE version 2.3.4 (Pritchard et al. 2000; Falush et al. 2003, 2007). For each STRUCTURE run, we used 50,000 Markov chain Monte Carlo replicates, applying an admixture ancestry model and assuming diploid loci. In each of the runs, the first 10,000 replicates were discarded as “burnin”, a method that reduces noise in the final sample by excluding samples that were taken while the model was still unstable (Hubisz et al. 2009). For each dataset, we did 10 replicate runs at each level of K from one to thirteen. The appropriate K for each dataset was inferred using the software program STRUCTURE HARVESTER (Earl and vonHoldt 2012), according to the method of Evanno et al. (2005).

Morphology Versus Genetics

To test whether genetic patterns match groupings of plants based on morphology, we used a Mantel test, carried out in the R package vegan (version 2.3-4). For the genetic data, we used the same identity by state distance matrix used for PCoA (see above). For the morphological data, we created a Euclidean distance matrix using the eight traits used for the morphometric analyses described above. We ran a partial Mantel test using both Pearson and Spearman correlation coefficients, and controlling for geographic distance.

Vegetative Morphometric Patterns

Patterns of morphological variation across subspecies, locales, and individuals were complex, as suggested by visualizations of multivariate analysis (Fig. 2A). In PCA, preliminary analysis suggested that only four variables contributed strongly to the model (margin hair type, short hair type, long hair density, long hair glands), and so only these variables were retained for the final analysis. In the final model (Fig. 2A), the first two principal components (those plotted) account for more than 84% of the variance in the data. Long hair density is strongly negatively correlated with the first principal component (vector loading: -0.93), followed closely by margin hair type, which is positively correlated (vector loading: 0.32). Long hair glands and margin hair type both contribute very strongly to the second principal component, but in opposite directions (vector loadings of 0.62 and -0.68, respectively).

Fig. 2. 

Results of morphometric analysis. A, results of principal components analysis (PCA); size of circle indicates the number of individuals. B, classification tree for manzanita individuals; character states at left define all individuals below that branch; the character ‘long hair density' appears in this tree twice because some sub-groups of plants are distinguished by very dense hairs (3) and some by sparse hairs (1). C, UPGMA tree for manzanita individuals; branch length represents morphological distance.

The classification tree (Fig. 2B) suggests some of the same patterns as the PCA, indicating that a large number of putative Del Mar manzanita individuals, including some from the type locality, have a lower density of long hairs than other sampled plants. As in the PCA, other sampled individuals are also differentiated by margin hair type and short hair presence; with the exception of long hair glands, the most informative characters are the same in both the PCA and the classification tree.

The UPGMA tree (Fig. 2C) allows for examination of relationships among individual plants. This tree shows that most groupings of individual plants do not conform to expectations based on the morphology of the plants, or even the locality where they were collected (Table 1). Subspecies do not form cohesive groups at any level of the tree, and in many cases individuals from the same locale are not cohesive. Finally, collections from the type locality of Del Mar manzanita are not strongly cohesive and are found in multiple portions of the tree. Overall, these results demonstrate the morphological heterogeneity of the sample, with morphologically distinct plants co-occurring at the same locale.

Genetic Variants

A total of 54 samples were successfully sequenced (Table 1). Quality was high for all samples; none were excluded due to low DNA sequencing coverage or excessive missing data. Sequencing rates averaged 6,903,818 reads per individual (SD = 3,184,037). For information on sequencing coverage for each plant, see Appendix S2 (Dryad Digital Repository). Raw DNA sequence data is on the NCBI BioProject database under accession PRJNA396085.

A total of 16,083,660 reads were available for the sample used to construct the RAD reference genome (D. Burge 2071_4, A. glandulosa subsp. crassifolia, Table 1). 15,606,332 reads passed quality filters and were assembled into 25,321 contigs 92 of bp in length, for an assembled RAD reference genome length of 2,329,532 bp.

Using the RAD reference genome as a guide, we called a total of 163,793 loci. The complete variant set is available at the Dryad Digital Repository (Appendix S3; doi:  10.5061/dryad.kv573c7). After the application of post-variant calling filters, the Complete dataset contained 60,865 loci for 54 individuals; the Reduced dataset contained 53,850 loci for 39 individuals. These datasets are available at the Dryad Digital Repository (Appendices S4 & S5; doi: 10.5061/dryad.kv573c7).

Population Genetics

Principal coordinates analysis of the Complete dataset (Fig. 3) revealed several distinct groups of plants, including: 1) a group comprising all the putative Del Mar manzanitas collected at the type locality, plus several nearby locales (Fig. 3A), 2) a group comprising all remaining samples of putative Del Mar manzanita, plus two individuals of Eastwood manzanita and one of A. glandulosa subsp. atumescens (Fig. 3B), and 3) a group containing most of the remaining Eastwood manzanitas collected in northern California (Fig. 3C). STRUCTURE analysis of the Complete dataset revealed that the most optimal number of genetic groups was two (K = 2). The STRUCTURE HARVESTER output demonstrating this result is available at the Dryad Digital Repository (doi:  10.5061/dryad.kv573c7; Appendix S6). The ancestry proportions inferred by STRUCTURE (Fig. 4) reveal that most of the plants identified as Del Mar manzanita are dominated by one genetic group (Fig. 4, Ancestral Group AG1), while most of the plants identified as Eastwood manzanita are dominated by the second genetic group (Fig. 4, Ancestral Group AG2). Results of the STRUCTURE analysis also suggests some genetic admixture between these groups in the case of ten plants (Fig. 4, 2006-1, 1709-6, 1717-4, 1717-5, 1717-9, 1688-1, 1688-3, 1688-4, 2087-6, & 2090-7). All of these plants are from MCAS Miramar.

Fig. 3. 

Genetic relationships based on principal coordinates analysis (PCoA) of all SNP data. Individuals from the type locality of Del Mar manzanita (1729; Del Mar) are indicated with a dashed halo. For clarity of presentation, some of the plants are not labeled individually.

Fig. 4. 

Proportional genetic ancestry based on analysis of all SNP data using STRUCTURE. Results of analysis with the program STRUCTURE, assuming two ancestral groups (see Methods). Under the assumption that two ancestral groups are present, STRUCTURE estimated the proportion of genetic variation in each sampled individual that was assignable to each of these groups. Proportions are indicated by shaded columns.

PCoA of the Reduced dataset (Fig. 5) revealed groups that match those revealed by the Complete dataset, including: 1) a group comprising all the putative Del Mar manzanitas collected at the type locality, plus several nearby locales (Fig. 5A), and 2) a set of three weakly supported groups comprising all remaining samples of putative Del Mar manzanita (Fig. 5B). STRUCTURE analysis of the Reduced dataset revealed that the most optimal number of genetic groups was two (K = 2). The STRUCTURE HARVESTER output demonstrating this result is available at the Dryad Digital Repository (doi:  10.5061/dryad.kv573c7; Appendix S7). The ancestry proportions inferred by STRUCTURE (Fig. 6) reveal that putative Del Mar manzanita individuals come from two fairly distinct genetic groups (Fig. 6, Ancestral Groups AGC1 and 2). STRUCTURE results also suggest some genetic admixture between these ancestral groups, represented by six plants that were inferred to be of hybrid origin (P << 0.01; Fig. 4, 1719-25, 1719-29, 1722-3, 1722-4, 1722-7, and 1725-3). All of these plants are from MCAS Miramar.

Fig. 5. 

Genetic relationships based on PCoA of SNP data for only putative Del Mar manzanita. Individuals from the type locality of Del Mar manzanita (1729; Del Mar) are indicated with a dashed halo. For clarity of presentation, some of the plants are not labeled individually.

Fig. 6. 

Proportional genetic ancestry for putative Del Mar manzanitas based on STRUCTURE. Results of analysis with the program STRUCTURE, assuming two ancestral groups (see Methods). Under the assumption that two ancestral groups are present, STRUCTURE estimated the proportion of genetic variation in each sampled individual that was assignable to each of these groups. Proportions are indicated by shaded columns. Individuals marked with an asterisk are considered to be meaningfully admixed at the P < 0.05 threshold.

Genetics Versus Morphology

The Mantel test using the Pearson correlation coefficient indicated a significant relationship between morphology and genetics (Mantel statistic r = 0.094; P = 0.005), as did the test using the Spearman correlation coefficient (Mantel statistic r = 0.069; P = 0.022).

Does Morphology Predict Genetic Relationships?

Overall, our Mantel test results suggest that morphology is a reliable predictor of the underlying genetic groups. However, there is a very weak association between genetic groups and the names that were assigned to the sampled plants based on current morphological circumscriptions of A. glandulosa subspecies (Figs. 3 and 4). This lack of correspondence suggests that the current circumscriptions should be modified. Unfortunately, our sampling of A. glandulosa subspecies and populations is not broad enough to provide the evidence necessary to support such changes. Future studies should aim to expand on the sampling employed in this study, ideally including a geographically broad sample of A. glandulosa populations from all of the known subspecies.

Our morphological analyses did not use fruit characters, due to the time of year in which we sampled. Fruit and seed characters are frequently used in manzanita systematics, and some putative Del Mar manzanita populations are known to have a distinctive (within A. glandulosa), flattened fruit (Parker et al. 2012), especially near the type locality, close to the sea (V. T. Parker, San Francisco State University, personal observation). Future research should aim to measure a broader suite of characters than those analyzed in this study, especially those relating to fruit and inflorescence morphology.

How Widespread is Del Mar Manzanita?

At the broadest scale, our genetic data suggest that the Del Mar manzanita forms part of a very widespread genetic lineage (Fig. 4, AG1) found mainly in southern California and northern Baja California, Mexico. This result argues in favor of a greatly expanded circumscription of Del Mar manzanita. In this case, Del Mar manzanita would be the dominant burl-forming manzanita in coastal portions of southern California and northern Baja California, Mexico (Fig. 4, AG1). On the other hand, the more focused analysis of only the samples from coastal San Diego County suggests that plants from the type locality and nearby locations close to the coast fall into one genetic group (Fig. 6, AGC1), while plants from more distant locations, especially more inland locations like MCAS Miramar, do not form part of this group. This result suggests that a more restrictive view of Del Mar manzanita might be possible; the circumscription could be modified to include only those plants from the second analysis that fall into the same genetic group as the plants from the type locality (Fig. 6, AGC1), which are from sandstone derived soils very close to the sea. Unfortunately, we did not obtain a broad enough sample of A. glandulosa subspecies and populations to reliably test taxonomic concepts in the group based on DNA, and so we are not able to modify circumscriptions according to genetic patterns. More research is needed to determine taxonomic limits in the subspecies of A. glandulosa. Ideally, future research will expand on the present work by adding samples of more subspecies and populations from especially from coastal San Diego County and northern Baja California, Mexico.

Conservation and Recovery Implications

Management of rare plants demands consistent, objective tools to identify these plants to the exclusion of other taxa. Identification is usually done using morphological features. However, many rare plants are difficult to identify morphologically, or only display their diagnostic features during a particular phenological period (e.g., during flowering or fruiting). There are a variety of methods available that seek to overcome this limitation using genetic markers (reviewed by Hollingsworth et al. 2016). However, many of these techniques, particularly DNA barcoding methods (Kress et al. 2005), fall short in groups with low levels of genetic variation. This situation limits objective identification of rare plants, which in turn limits management, where precise identification of rare plants is critical for effective conservation. As demonstrated here and in other recent work (Hollingsworth et al. 2016; Andrews et al. 2016), high-throughput DNA sequencing, particularly RADseq, can be used to target large numbers of variable SNPs without any prior knowledge of the SNPs or the genome of the target organism, allowing for the rapid development of useful panels of genetic markers for precise and repeatable conservation work.

Our results show that in the case of Del Mar manzanita and its close relatives, there is a mismatch between genetic groups and the groups based on current taxonomic concepts. As we explained above, our genetic results could be used equally effectively to argue in favor of either a greatly reduced or a greatly expanded circumscription of Del Mar manzanita. Clearly, these two alternatives have dramatically different implications for conservation and recovery; if we were to choose the first option, this already rare taxon would become even more rare; if we were to choose the second option, the taxon would become very widespread, and would probably not require conservation measures. As outlined above, future research should expand upon our work by sampling more populations and taxa, in order to thoroughly test the circumscriptions of the A. glandulosa subspecies. Such work will lead to conservation outcomes that are consistent with the underlying genetics of the plants being conserved, rather than the names applied to them.