Conservation genomics of Wyethia reticulata (Asteraceae)

Wyethia reticulata is a rare perennial herb found only on gabbro-derived soils of the Pine Hill formation in the Sierra Nevada foothills of El Dorado County, California. Wyethia reticulata is capable of both vegetative (clonal) and sexual reproduction, although the former is thought to be more common than the latter, with sexual reproduction tied to fire. The potential dominance of vegetative reproduction has conservation implications for W. reticulata, as populations could appear healthy in terms of stand density, while losing genetic diversity due to a lack of sexual recruitment. Past genetic assessments of W. reticulata relied on low-resolution markers from a small number of populations and did not directly address conservation issues. Here, double-digest random amplified DNA sequencing was used to develop a high-resolution genomic dataset for a sample of 168 W. reticulata samples from seven locations that represent every known population center of the species. The data was used to 1) calculate the frequency of asexual reproduction, 2) determine which populations contain unique genetic variation, and 3) test whether there are any ecological correlates of genetic patterns. In contrast to past genetic research, the results presented here indicate that clonality is not common; less than 2% of plants were identified as being derived from asexual reproduction. Nevertheless, some locations seem to be more affected by asexuality than others; the index of association varied from 4.4 to 19.4 (average = 10.9). Genetic differentiation among locations was low (Fst 0.03-0.04). However, clustering identified three distinct population segments that were not previously recognized. Elevation was the only environmental variable significantly associated with genetic differentiation. Based on these results, it is recommended that conservation effort be focused on genetically isolated population segments.Funding provided by: U.S. Bureau of Land ManagementCrossref Funder Registry ID: http://dx.doi.org/10.13039/100007149Award Number: L16PX01835Sampling For sampling of W. reticulata, seven locations in the Pine Hill area were selected (Table 1; Fig. 1). Locations were selected to represent every major population center of the species and the range of environmental conditions under which the species is found. Locations were selected based on information obtained during preliminary visits to known population centers carried out between 15 and 19 June, 2018. Preliminary visits allowed the author to select large populations within which to implement the sampling scheme outlined below. To prevent confusion between populations and subpopulations, collecting locations are referred to as 'locales.' Locales are spatially discontinuous groups of plants with centers at least 1000 m apart and edges at least 100 m apart. Sampling was carried out using a spatially explicit, repeatable approach based on past research into the population biology of W. reticulata (Ayres and Ryan 1997) and a recent review of best practices for population genetic sampling of plants with mixed clonal and sexual reproduction(Arnaud-Haond et al. 2007b), as in W. reticulata Ayres and Ryan (1997, p. 765) sampled three contiguous patches of W. reticulata with a total area of 900 m2 and identified a total of 9 genets in this area; this implies that on average, there is one genet present in a 100 m2 patch, meaning t­­­hat in an area of continuous cover, genets are separated, on average, by 5.6 m (the radius of a circle with an area of 100 m2). Ayres and Ryan (1997, p. 770) also found significant genetic differentiation among populations separated by as little as 500 m. Based on this existing information for the focal species, a modification of the stratified random sampling design recommended by Arnaud-Haond et al. (2007b) was used, as this method was developed for species with patchy cover. A custom R script (R Development Core Team 2019) written in by the author was used to obtain 48 random GPS points for each locale; the points were at least 5.6 m apart inside a 300 m square. During a locale visit, each of the GPS points was visited. At some locales, visits to the points required the use of protective clothing and knee pads in order to access plants and avoid injury from the dense chaparral. If a colony was present at a point, a sample was obtained there. This continued until all 48 points were visited. If the 48 points yielded less than 24 samples, the points were re-visited in random order, and the nearest colony sampled, selecting a stem at least 5.6 m distant from the nearest sampled stem. This was continued until 24 samples were obtained. At locales other than Square Nail (Table 1; Fig. 1), this strategy yielded sufficient samples without difficulties. Due to the small number of plants at the Square Nail Locale, most of which were associated with a disused fire road (Fig. 1), an additional eight stems were sampled to the south of the originally targeted area (Fig. 1); samples at this location were haphazardly selected, no less than 5.6 m apart. At each locale, up to four duplicate samples were obtained to be used as biological replicates. Duplicate samples were obtained from randomly selected stems, and were taken from the same flower on that stem. For all samples, flower petals of a single immature flower (anthers not yet dehisced) were taken. Tissue was preserved in two mL of absolute ethanol. Flower petals were selected as the target tissue because preliminary research on DNA extraction comparing immature leaf to petal tissue (unpublished data) suggested that only flower petals yielded DNA in sufficient quantities for the successful implementation of the ddRAD method employed here. Due to the highly patchy nature of W. reticulata populations (Ayres and Ryan 1997), in which demography is known to be strongly influenced by past fires and mechanical disturbance (Ayres and Ryan 1997, Ayres 1997), the resulting samples were not evenly distributed across the targeted areas of habitat (Fig. 1). What is more, the requirement for flowering stems occasionally meant that non-flowering stems were passed over for sampling. Both the location and flowering status of stems can be affected by both fire and mechanical disturbance history (Ayres and Ryan 1997, Ayres 1997). However, other than fire history at a locale (see below), the influence of non-random historical factors such as mechanical disturbance was not considered in the sampling due to the difficulty of measuring and quantifying such disturbances. For example, at the Square Nail locale few plants were found, even under chaparral that had re-grown since a prescribed fire that burned in 1990 (Ayres 1997); nearly all plants (flowering or not) located during extensive visits to the area prior to sampling were located within 10 m of the fire road (unpublished data; Fig. 1). Environmental Data Climate and soil data for each sampled plant were inferred for use in tests of association between genetics and environment (see Population Genetics, below). The following six climatic variables were selected to encompass important drivers of vegetation patterns (Flint et al. 2013): 1) climatic water deficit, 2) actual evapotranspiration, 3) average minimum winter temperature, 4) average maximum summer temperature, 5) annual precipitation, and 6) elevation. Using the R package raster (Hijmans 2017), GIS data from Flint et al. (2013) were used to obtain the first five climatic variables for each sampled plant. Average minimum winter temperature was calculated by taking the average of the monthly minima for December, January, and February. The average maximum summer temperature was calculated by taking the average of the monthly maxima for June, July, and August. Climatic water deficit and actual evapotranspiration are described by Hijmans (2017). Elevation was taken using a Garmin 12 GPS unit (Olathe, Kansas). Soil data was taken from GIS data layers in the Soil Survey Geographic (SSURGO) database for El Dorado Area, California (U.S. Department of Agriculture, Natural Resources Conservation Service; https://www.nrcs.usda.gov/wps/portal/nrcs/main/soils/survey/). For each locale, the map unit name (MUNAME) for the soil type underlying the area was recorded. For all raw environmental data, see Appendix S1 (Dryad Digital Repository: doi pending). Information on past fires was obtained from the CAL FIRE Resource Assessment Program Fire Perimeters (http://frap.fire.ca.gov), which contains data on all fires of greater than ten acres in the state of California. The data covers all such fires since 1898. If a fire burned all or part of the locale, the date of the fire was recorded. For fires smaller than ten acres, records maintained by the Motherlode Office of the Bureau of Land Management were examined; no such fires recorded by the office had affected the targeted sampling areas. DNA Extraction and Sequencing Genomic DNA was extracted from the flower petals by Eurofins BioDiagnostics (River Falls, Wisconsin) using the LGC sbeadex 96 well plant kit. Two kits were used. In order to avoid systematic DNA-quality biases among locales based on extraction, a blocking and randomization scheme was employed. Samples were blocked by locale, assigning equal numbers of individuals from each locale to each of the kits, and then randomly assigning each individual to a well within a kit. Extraction was done according to the manufacturer's instructions, with several proprietary modifications developed by Eurofins BioDiagnostics. DNA was checked for degradation using a 1% agarose gel, and quantified using a spectrophotometer. Total genomic DNA was standardized to 30 ng/mL by Eurofins BioDiagnostics and sent directly to Floragenex, Inc. (Portland, Oregon) for further quality control, library construction, and DNA sequencing. The library preparation and sequencing methods employed by Floragenex followed those described by Lozier (2014), with the exception that the present work used the restriction enzyme PstI (Hipp et al. 2014). In brief, sequence identifier barcodes and sequence adapters were added to the genomic DNA after digestion by the endonuclease PstI. Barcoded samples were combined and the fragments sequenced outwards from restriction-sites using single-end DNA sequencing reads 100 bp long. All DNA sequencing was done using a run of a HiSeq 2500 DNA sequencing instrument (Illumina, San Diego, CA). Following sequencing, DNA sequences from individual plants were separated based on their barcodes using the program fastq-multx (Aronesty 2013); the barcodes were then removed from the sequences. Variant Detection Due to the lack of an existing reference genome for Wyethia or a closely related member of Asteraceae, detection of variation among individual plants was done using a 'RAD reference' approach (Lozier 2014), in which a kind of reference genome was 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 for alignment of sequences from other individuals. The plant selected for use as the reference was 2344-22 (Table 1; Appendix S1, Dryad Digital Repository: doi pending). In developing the RAD reference, custom methods developed by Floragenex were used to cluster identical sequences that had 5-500× sequencing coverage (Lozier 2014), 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), 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 mis-matched 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. Criteria for SNP calling required a minimum phred score of 10, a minimum of 6× sequence coverage and a maximum of 25% missing data across samples. Following SNP calling, the data was filtered based on the repeatability of genotype calls, as revealed by the 22 biological replicates included in the sequencing and SNP calling. Each of the loci was checked for genotype concordance between replicates. A locus was removed when the biological replicates had different genotype calls. Such discordant loci were removed in order to ensure that only repeatable loci were included in the final dataset. The process of identifying and removing discordant loci was carried out using a proprietary bash script (Floragenex, Inc.) that employed VCFtools (Danecek et al. 2011). Quantification of Clonality Previous research indicated that W. reticulata is strongly clonal (Ayres and Ryan 1999), reproducing asexually via underground rhizomes. As is the case with all organisms capable of clonal reproduction, this can lead to one genetic individual (genet) being represented by multiple independent vegetative parts (ramets) derived from that genet (Harper 1977). Thus, an effort was made to quantify clonality in the dataset. Clonality was quantified by identifying the number of unique multilocus genotypes (MLGs) and assessing whether these are the product of asexual reproduction. When an MLG is shared among multiple individuals, this does not mean that such individuals belong to one genet and are the product of asexual reproduction; the dataset may simply lack the power to detect the genetic differences between closely related plants derived from sexual reproduction. Conversely, it is possible that multiple MLGs belong to one genet, with subsequent somatic mutation. In the latter case, the group is better referred to as a multilocus lineage, or MLL. Overall, the definition of MLLs relies on selection of a genetic distance threshold to group MLGs into MLLs (reviewed by Bailleul et al. 2016). Some researchers have solved this problem by collecting genetic data on known ramets in order to directly quantify the amount of genetic differentiation expected among ramets (Douhovnikoff and Dodd 2003). However, such studies are very rarely done due to the logistical hurdles involved. More frequently, workers have relied on the following categories of methodology to group MLG: 1) clustering algorithms (Kamvar et al. 2015), 2) modeling asexual versus sexual reproduction (Bailleul et al. 2016), 3) visual examination of frequency distributions of genetic distance (Meirmans and Van Tienderen 2004; Arnaud‐Haond et al. 2007a; Clark and Jasieniuk 2011), and 4) various combinations of these methods. In all such methods, however, the definition of MLLs depends, either directly or indirectly, on assumptions made by the user. Methods to quantify clonality in the present study were selected to minimize the number of assumptions, especially with respect to population genetic or demographic parameters that were not directly measured. First, the R package poppr v. 2.8.3 (Kamvar et al. 2014) was used to group MLG into MLLs based on a-priori selection of a genetic distance threshold. A Euclidean distance matrix was created using the bitwise.dist function of poppr. The filter_stats algorithm was then applied to this matrix under default parameters, identifying three alternative genetic distances thresholds that could be used to identify MLLs. Requiring a minimum genetic distance to distinguish one MLL from another allows for small differences among MLGs, such as those derived from DNA sequencing error or somatic mutations. The mlg.filter option of poppr was used, under default settings, to collapse MLGs into MLLs based on each of the three genetic distance thresholds, farthest_neighbor (merges MLGs based on maximum distance between individuals in either of the clusters; most strict), average_neighbor (merges MLGs based on average distance between every pair of individuals between clusters; intermediate strictness), and nearest_neighbor (merges MLGs based on minimum distance between individuals in either cluster; least strict). Results of the three methods were similar: farthest_neighbor identified one MLL represented by six individuals from four locales (2340, 2343, 2344 & 2345); average_neighbor identified one MLL represented by five individuals from three locales (2343, 2344, & 2345); nearest_neighbor identified two MLLs represented by four individuals from three locales (2343, 2344, & 2345). To test sensitivity to the genetic distance methods used, the analysis was run again using dissimilarity distance (diss.dist in poppr) and Nei's genetic distance (nei.dist in poppr). The results were identical. The R package RClone v. 1.0.2 (Bailleul et al. 2016) was used as a graphical means to identify MLLs, and thereby separate ramets from genets among the sampled stems. The pairwise allele distances among the sampled stems were calculated (genet_dist function of RClone) and plotted as histograms for each local. As in Meirmans and Van Tienderen (2004) and Clark and Jasieniuk (2011), it was assumed that the frequency distributions of pairwise distances would each contain a wide valley separating stems derived from asexual reproduction (low end of the frequency distribution) from those derived from sexual reproduction (high end). However, no such valley was identified for any local; in all cases, the major cluster began at an allele distance of more than 1000 (Appendix S2, Dryad Digital Repository: doi pending). As in Meirmans and Van Tienderen (2004), this cluster of greater distances was inferred to belong to the group of samples derived from sexual reproduction; additional clusters were present in some histograms, but these were at greater distances. Overall, the graphical results were taken to indicate that no stems were derived from clonal reproduction; such asexual reproduction would be expected to result in a separate cluster of allele distances at the low end of the distribution (Meirmans and Van Tienderen 2004; Clark and Jasieniuk 2011), which was not observed. As a final alternative method of inferring the degree to which plants are reproducing asexually, the index of association (IA, Smith et al. 1993) and the index of multi-locus linkage disequilibrium (r d, Agapow and Burt 2001) were inferred for each population. Calculations were done in poppr (Kamvar et al. 2014). Both these indices express linkage disequilibrium. As such, both tend to be higher in populations where there is a higher level of asexual reproduction. The IA index is similar to r d, but the latter accounts for the number of loci sampled, while the former does not. Both of the indices suggest that asexual reproduction is taking place mainly in locales 2340, 2344 and 2345, (Table 2). This supports the results of the MLG-based analyses above, which identified a small number of potential clones from some of the same locales. Although the indices of association indicate that asexual reproduction may be taking place in some locales (Table 2), the MLL results suggest that nearly every sampled plant is a separate genetic individual. Thus, it appears that asexual reproduction (clonality) is not common in W. reticulata, at least at the scale (above 5.6 m horizontal distance) that sampling was done. Given the low rate of clonality (only two MLLs identified using the least conservative filtering algorithm of poppr; less than 1.2% of samples), asexual reproduction was ignored in subsequent population genetic analyses; all of the work described below was carried out on the complete dataset (168 samples; Table 1), excluding the biological replicates Genetic Differentiation Among Locales Population genetic analysis was used to obtain genetic diversity statistics for each locale and estimate the amount of differentiation among locales. The R packages adegenet, v. 1.3-1 (Jombart 2008; Jombart and Ahmed 2011) and poppr, v. 2.8.0 (Kamvar et al. 2014, 2015) were used to calculate population genetic statistics. The package adegenet was used to calculate allelic diversity and heterozygosity for each locale, and to estimate differentiation (Fst) among locales. The significance of differentiation was tested using Goudet's G-statistic Monte Carlo test, which was implemented in adegenet with 500 replicates (Jombart 2008). Isolation by distance at the level of the locale was assessed using partial Mantel tests carried out in the R package vegan v. 2.5-1 (Oksanen et al. 2018), controlling for the latitude and longitude of the locale centroid (Table 1); Pearson and Spearman correlation coefficients were used with 1000 replicates of the permutation test. As an additional means to quantify the relationship between genetics and spatial distance, the R package EcoGenetics v 1.2.1-5 (https://cran.r-project.org/package=EcoGenetics) was used to test for spatial autocorrelation; the eco.malecot algorithm was implemented, testing for global autocorrelation between great circle distance and Kinship (Nason's Fij; Loiselle et al. 1995) Clustering Analysis The program STRUCTURE, v. 2.3.4 (Pritchard et al. 2000; Falush et al. 2003, 2007; Hubisz et al. 2009) was used to cluster the sampled individuals. 100,000 Markov chain Monte Carlo (MCMC) replicates were run, discarding 20,000 replicates as burnin (Hubisz et al. 2009). The analysis was repeated eight times, varying K (the number of theoretical clusters) from one to eight, with 10 MCMC replicates at each K. The optimal K was inferred using STRUCTURE HARVESTER (Earl and vonHoldt 2012), according to the method of