A Comparison of Non-Destructive Visceral Swab and Tissue Biopsy Sampling Methods for Genotyping-by-Sequencing in the Freshwater Mussel Fusconaia askewi

1. Introduction

Studies examining the population genetic structure - or lack thereof - of rare and threatened species are critical tools for enabling biologists to recommend management strategies of these organisms [1,2,3,4,5]. Large sample sizes are often needed to obtain robust estimates of population genetic parameters, and it is important to minimize the degree to which genetic sampling results in the disruption or death of individual organisms, especially for rare species. A number of nondestructive sampling techniques (e.g., saliva, hair, swabs, etc.) have thus been developed to reduce the potentially harmful impact of genetic sampling [6,7,8,9,10].

In North America, Family Unionidae is a speciose group of freshwater mussels consisting of both widely distributed and locally endemic species, many of which exhibit varying conservation statuses [11,12,13]. A diverse array of DNA collection methods has been proposed for this group that reduce the mortality of sampled and re-released individuals. Collection methods include clipping of small pieces of mantle or foot tissue (i.e., “tissue biopsies”), swabbing of the viscera, foot scraping, and haemolymph extraction via hypodermic needles [14,15,16]. These four methods were specifically evaluated for the Freshwater Pearl Mussel (Margaritifera margaritifera) with respect to survivorship, growth characteristics, and DNA quantity and quality of the sampled individuals [16]. All four methods reliably produced DNA, though yields from haemolymph and foot scraping were significantly lower than those produced from tissue biopsies and swabbing. Furthermore, survivorship (after 128 days) was 100% for all methods, leading the authors to conclude that visceral swabbing was the most appropriate and least invasive method for tissue collection in M. margaritifera. A separate investigation in the Pimpleback mussel Quadrula pustulosa further determined that tissue swabbing was preferred for DNA collection over mantle clips because it produced quality samples with less disruption to the organisms [15].

Visceral swabbing methods have produced quality DNA appropriate for use in traditional Sanger sequencing methods that include generating ND-1 mitochondrial sequences [15] and microsatellite loci [16]. However, such methods do not require large amounts of targeted DNA be extracted, and contamination with other non-target organisms is largely not problematic because primers are oftentimes designed to be species-specific. The advent of genotyping-by-sequencing (GBS) methods, on the other hand, has allowed for the production of single nucleotide polymorphism (SNP) data that can routinely produce orders of magnitude more genetic markers at much lower cost compared to traditional Sanger sequencing methods [17]. Massault et al. [18] evaluated the effectiveness of swabbing methods for GBS sequencing in the Silver-Lipped Pearl Oyster (Pinctada maxima, Family: Margaritaferidae) and found such methods to be appropriate for parentage analysis and pedigree reconstruction in these marine bivalves. Swab collections produced DNA of sufficient quantity and quality for genotyping, and SNP counts were highly correlated among tissue and swab collection methods [18]. However, it remains unclear whether swab methods are similarly appropriate for GBS sequencing of freshwater Unionid mussels, and further investigation is needed to address this question.

In the current study, we evaluated the degree to which visceral swab and tissue biopsies generate high quality GBS/SNP genotype data for the Unionoid freshwater mussel Texas Pigtoe (Fusconaia askewi). This was achieved by collecting and processing swab and tissue biopsy samples from the same individuals and generating separate SNP datasets for both types of sampling methods, with swab samples being collected prior to performing tissue biopsies. Our primary goal was to determine whether the less-disruptive swab sampling method was sufficient to undertake large-scale population genetic studies using GBS methods for freshwater mussels. Therefore, we sought to 1) compare DNA concentrations and number of raw reads produced from visceral swab samples and tissue biopsies collected across 14 different individuals, 2) determine whether and how SNP coverage and numbers differed across the two sampling methods and 3) determine whether the patterns of genetic variation among individuals are consistent across both datasets when assembled independently from one another, as well as when assembled to a larger tissue-based de novo reference assembly.

2. Materials and Methods

2.1. Sample Collection

Fourteen mussels were collected via snorkeling and tactile searches at two separate locations along the Sabine River in Texas, USA at 32.62986° North 32.62986° West (n = 6), and at 32.462220° North and 94.845864° West (n = 8). The mussels were identified in the field as F. askewi on the basis of morphology and were subsequently sampled using Isohelix SK-1S buccal swabs by carefully rubbing both sides of the swab along the mussel’s foot tissue for 30 seconds. Swab heads were immediately placed in individual containment tubes and kept on dry ice until storage at −20 °C. After swabbing, a ~0.5 cm³ tissue plug was collected from the same individual using a nasal biopsy tool, and the mussel was then returned to its original location. Tissue samples were stored separately in 95% ethanol on dry ice before transferring to a −20 °C freezer for longer-term storage.

2.2. DNA Extraction and Quantification, Library Prep, and de novo Assemblies

DNA was extracted from the 14 swab samples and the replicated 14 foot-tissue plugs using Qiagen DNeasy blood and tissue kits in 96-well format. These samples were part of a larger group of 378 mussel tissue samples that were used in a separate population genetic study [19]. No attempts were made to equalize DNA concentrations prior to DNA library construction. However, DNA concentrations of replicated swab and biopsied tissue samples were quantified post-hoc to investigate whether initial concentrations affected the number of reads per individual in the resulting libraries. DNA concentrations were determined using the Qubit^® dsDNA HS Assay (Invitrogen, USA) kit following standard protocols described in the Qubit manual. The DNA concentration of each of the 14 swab and tissue replicates was measured from a single aliquot three separate times, and a mean concentration was calculated for each sample.

Extracted DNA from the broader collection of mussels (n = 378 [19]) was used to create a reduced-complexity genomic library. This library included the 14 tissue and 14 swab replicates utilized here and was generated following modified protocols commonly used in our laboratory group across a wide diversity of species [5,19,20,21,22]. Restriction enzymes EcoRI and MseI were used to digest extracted DNA; EcoRI adapters (i.e., 10–20 base-pair multiplex identifier sequences: MIDs) were ligated onto the resulting fragments, and the 14 swab and 14 tissue replicates were each assigned different unique barcodes. Labeled products were amplified through two rounds of PCR using Illumina primers. PCR products were pooled into a single library and sent to the University of Texas Genomic Sequencing and Analysis Facility (Austin, TX) for sequencing on an Illumina Novaseq SR100 platform.

A number of data processing steps were used to ensure high-quality sequencing data. First, Bowtie v.3 was used to identify PhiX control sequences, and reads that assembled to the PhiX genome were removed [23]. Custom Perl scripts were used to match sample IDs with unique barcode identifiers and remove Mse1 adapters and barcodes from sequence reads. The resulting sequence data were then organized into different assemblies for separate analyses. Since no reference genome is available for F. askewi, three separate de novo assemblies were built using different sets of individuals and/or collection methods. The first reference genome was created using a much larger dataset of tissue-sampled F. askewi individuals (n = 96, [5]). We used this reference to align reads for both tissue and swab replicates, enabling a direct comparison of resulting SNP datasets (i.e., to directly assess coverage, depth and genetic variation of “tissue” and “swab” samples of the same 14 individuals). Two additional reference assemblies were also separately created; one exclusively using tissue samples (tissue-only assembly), and another for the replicated swab samples (swab-only assembly). Each of the three separate reference assemblies (i.e., control reference, tissue-only, and swab-only) followed the same processes using part of the dDocent variant calling pipeline [24]. First, unique reads were identified for each individual, and reads with less than four copies and shared across fewer than four individuals were filtered out. Subsequently, the remaining sequences (i.e., those that had more than four copies and were shared across more than four individuals) were assembled to form each reference using CD-hit by utilizing a similarity threshold of 80% [25,26].

2.3. Number of SNPs and Coverage

Sequence reads were assembled to reference scaffolds using the aln and samse algorithms from BWA (version 0.7.13-r1126; [27]). BCFtools (version 1.9) was used to both identify SNPs as well as to calculate Bayesian genotype likelihoods for each SNP [28]. Biallelic SNPs must have been represented by at least 80% of the individuals in a given assembly and have had a mean sequence depth of ≥2Xn (where n = number of individuals) in order to be included in the respective dataset. A custom Perl script was used to filter out potentially paralogous loci with exceptionally high sequence depth (mean sequence depth > assembly-wide mean + 2*sd). This script also filters loci based on mapping quality (minimum score of 30), as well as the difference in base and mapping quality between the reference and alternative alleles using Mann-Whitney U tests (z-score cutoff = 1.96). Genotype likelihoods were assigned to each SNP for each individual and used to estimate allele frequencies; SNP’s with a minor allele frequency of <0.05 were not included in the final datasets. If there was more than one SNP identified on a contig, only a single SNP was randomly chosen to include for analyses in order to reduce effects of linkage disequilibrium.

2.4. Estimation of Genotype Probabilities

Rather than “calling” SNPs, we utilize genotype probabilities which account for uncertainty associated with sequencing. Genotype probabilities were estimated using Entropy, a hierarchical model that estimates Bayesian-based admixture proportions and genotype probabilities for each individual for a predetermined number of populations (k) [5,29]. Calculation of posterior distributions for k = 2 was done by combining two separate Markov Chain Monte Carlo (MCMC) runs including 50,000 total iterations with a burn-in of 5,000 and sampling every 10th iteration. Genotype probabilities were averaged across both runs after checking for chain convergence with Gelman-Rubin diagnostic statistics.

2.5. Data Analysis

We performed a paired t-test to determine whether there was a significant difference between DNA concentrations of tissue and swab samples. We then used an ANCOVA to assess the relationship between DNA concentrations and the total number of raw reads generated, while accounting for sample type. This ANCOVA was also used to determine whether there was a significant difference in the numbers of reads that each sampling method produced.

We ran Principal Components Analysis (PCA) on genotype probabilities to determine whether SNP alignments to either the sample-type reference or the tissue control reference influenced resulting patterns of genetic variation. In total, we ran four PCAs with: 1) swab data aligned to the swab de novo assembly, 2) tissue data aligned to the tissue de novo assembly, 3) swab data aligned to the control tissue de novo assembly, and 4) tissue data aligned to the control tissue de novo assembly. Rather than focusing on exact genetic relationships among individual samples, analyses were primarily focused on how relationships among individuals changed given the sampling type and the references to which they were aligned. To examine this, we ran Procrustes analysis from the vegan package in Program R and report the Procrustes R, a measure of correlation between two multivariate plots [30]. The first Procrustes analysis was run comparing the two PCAs which were sample-specific; swab and tissue data aligned to their specific de novo assemblies. The second Procrustes analysis compared the two sampling types when aligned to the control reference. We posited that if swab samples produced a SNP dataset with similar genetic patterns to the tissue samples, even when aligned to separate de novo assemblies, the PCAs would be highly (>0.95) correlated.

3. Results

3.1. DNA concentrations

Mean DNA concentrations significantly differed among tissue and swab samples (two-tailed paired t-test, t = 2.171, p = 0.049, df = 13) with tissue DNA concentrations (mean ± 1 SE; 36.2 ng/ml ± 14.6) having 38% greater DNA concentrations than those of swab samples (27.8 ng/ml ± 12.5; Table S1). However, after omitting one swab sample that revealed exceptionally low amounts of DNA (0.02 ng/ml ± 0.03 ng/ml), DNA concentrations did not significantly differ (two-tailed paired t-test, p = 0.0996, t = 1.785, df = 12).

3.2. SNP numbers, Coverage, and Correlations

The mean number of sequencing reads (± 1 S.D.) generated per individual from tissue samples (1,608,863 ± 307,130) was 2.1X more than the number of reads generated from swab samples (518,985 ± 238,072), and this was a significant difference (Table S2; t = 9.693, p < 0.001). DNA concentration was not predictive of the number of raw reads that were generated regardless of the DNA sampling technique (ANCOVA; t = 1.2003, p = 0.284).

When both datasets were aligned to the control-reference, tissue samples produced more SNPs for each individual than swab samples, but at lower coverage. For the tissue dataset, there were 36,764 SNPs identified, each with a mean coverage of 3.64 ± 0.74 reads per SNP, whereas for the swab dataset, an order of magnitude fewer SNPs (n = 3,138) were identified with an average coverage of 4.29 ± 1.94 reads per SNP (Table S3). For the tissue-only assembly (i.e., tissue samples aligned to tissue-only reference), mean coverage of the 8,111 SNPs was 6.79 ± 1.48, whereas for the swab-only assembly (i.e., swab samples aligned to swab-only reference), mean coverage for the 565 SNPs was dramatically higher at 22.64 ± 10.5 (Table S4).

3.3. Comparing PCA Plots

None of the four PCA plots revealed genetic structure among individuals sampled from the two separate locations along the Sabine River (Supplemental Figure S1). When each sample type was aligned to their respective assembly (i.e., tissue data aligned to swab-only reference and swab samples aligned to tissue-only reference), there was a high correlation between the two PCA plots (Procrustes R = 0.982 Figure 1A). Similarly, when tissue and swab data were aligned to the control reference, the PCAs are also highly correlated (Procrustes R = 0.975; Figure 1B).

4. Discussion

Identifying DNA collection methods that minimize mortality and disruption to organisms is important when taxa being studied are rare or of high conservation priority [6,7,8,9,10]. Our results demonstrated that both buccal swabs and tissue biopsy methods produce quality DNA in the Unionid mussel F. askewi. Of the 14 individuals sampled, there was one notable exception where the DNA concentration for a single swabbed sample was near-zero (two of three DNA quantification measurements on the qubit system were “zero” and one measurement was at the lowest possible reading of 0.1 ng/µL). As we did not re-extract the sample, it is unknown whether this was a failure of the swab collection method itself or whether this was simply a failed DNA extraction. It is important to note that, excepting the potentially failed extraction, should the quantity of DNA collected be important for the study being conducted, swabs produce about 20 percent less DNA than tissue biopsies (though this was not significantly less in the current study [p = 0.097]), and this may need to be taken under advisement should other experimental protocols require higher concentrations of DNA.

Both methods of DNA collection studied here have proven to be effective at producing high-quality DNA for Sanger sequencing in other species of freshwater mussels [14,15]. In the current study, GBS sequencing of tissue samples generated twice as many raw reads per individual than swab samples. Tissue samples also generated more SNPs per individual regardless of whether the reads were aligned to the reference assembly or the tissue-only de novo assembly, but the mean read depth of each tissue-based SNP was lower than that of swab samples. Regardless of the sampling method, however, the Procrustes analysis revealed nearly identical relationships among all of the included samples. Thus, the less invasive swab method of DNA collection is likely to be an acceptable sampling method for performing population genomic work in Unionid mussels. However, it should be kept in mind that these results should be interpreted within the context of their sampling distribution. The genetic diversity of the individuals collected from these two closely-located sampling locales (separated by ~97 river km) within the Sabine River is likely lower than what might be expected in individuals collected from a larger sample size across a wider sampling distribution. It is therefore possible that the extra genomic data provided from tissue biopsy samples may provide a substantial benefit when it comes to identifying the pattern and extent of structure of more diverse collections.

We found that visceral swabbing and tissue biopsies reveal similar patterns of genetic variation among individual samples of F. askewi when performing GBS. Such methods are potentially transferable to other mussel systems as well [16]. Future researchers should consider the potential advantages and disadvantages of both methods. Here, tissue samples provided substantially more genomic coverage (i.e., more total SNPs) than viscera swabbing, yet the mean sequence coverage per SNP (and the concomitant increase in accuracy in genotype probabilities or genotype “calls” at each of these fewer SNPs) is higher for the viscera swabbing method. Researchers should not only keep in mind these potential tradeoffs, but also take into consideration the degree to which these different methods have the potential to disrupt growth and survivorship of these species.