Microsatellite loci discovery from next-generation sequencing data and loci characterization in the epizoic barnacle Chelonibia testudinaria (Linnaeus, 1758)

Specimen collections

Specimen collections of the Atlantic lineage took place on Nannygoat Beach, Sapelo Island, GA, USA (31.48 °N, 81.24 °W) between 2012 and 2014 under the collection permit of the University of Georgia Marine Institute, and sanctioned by the Georgia DNR Wildlife Services. We chose to collect from the horseshoe crab Limulus polyphemus (Linnaeus, 1758) because it is relatively abundant and easy to sample: each spring and early summer, horseshoe crabs crawl onto beaches to mate and lay their eggs. During this process, we removed one individual of C. testudinaria per host individual with a sharp knife directly on the beach, and preserved it in 95% EtOH immediately after collection. We collected specimens of the Indo-West Pacific lineage in the vicinity of Townsville, Queensland, Australia (23 °S, 143 °E), in September 2012 from green turtles (Chelonia mydas Brongniart, 1800). Working with officials from the Queensland Department of Environment and Heritage Protection, turtles were captured in-water for routine tagging and release. During capture barnacles were collected and immediately preserved in 95% EtOH.

Microsatellite loci discovery

We extracted genomic DNA from the feeding appendage of a single large hermaphroditic C. testudinaria collected from a horseshoe crab with Gentra Puregene Tissue Kit (Qiagen), and measured DNA concentration with a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). Genomic DNA was fragmented into approximately 700 bp lengths (insert size) and shotgun-sequenced on an Illumina MiSeq sequencer (PE250). We quality-checked paired-end reads with FastQC (Andrews, 2015). The software FastqMcf was used to trim adapters, cut low quality ends and remove low quality reads and their mate-pair read (Aronesty, 2011).

We executed the perl script palfinder to identify short sequence repeat regions (Castoe et al., 2012). The script also calls Primer3, version 2.0.0, to identify potential primer pairs that span the repeat region (Rozen & Skaletsky, 2000). The minimum number of repeat units was chosen as in Castoe et al. (2012). A repeat unit, also called kmer, is defined as the length of the short sequence repeat. For example, a dimer would be a repeat of two base pairs (e.g. GC), and a tetramer would be a repeat of four basepairs (e.g. AGGT). Primer3 parameters were the default values. The search resulted in a large number of potentially amplifiable loci (PALs), repeat regions for which primers were identified. We filtered the results by removing all PALs which occurred less than two times and more than the estimated genome coverage in the genomic reads based on the following reasoning: If the number of primer occurrences is low, the primer sequence may contain sequencing error. If the number of occurrences is higher than the expected genome coverage, the primer region may occur more than once in the genome, leading to amplification of multiple loci (genomic regions). Neither of these outcomes is desirable because a good marker occurs only once in the genome, and has a primer sequence that matches the genomic sequence well. We estimated genome coverage by mapping our genomic reads to 52 nuclear single-copy gene fragments available from the acorn barnacle Semibalanus balanoides (Regier et al., 2010) with the Geneious version 6.0.3 Read Mapper, using the default settings (Kearse et al., 2012), calculated the median coverage for each of the 52 gene fragments, and the grand median of all genes. R scripts for screening palfinder output as well as calculating genome coverage are available as Supplementary Information.

Of the filtered PALs, we chose 48 PALs for trial amplification, which differed in kmer length, kmer motif (e.g. AG vs TG) and fragment size. We extracted and amplified DNA of 16 C. testudinaria individuals for trials. DNA was extracted from feeding appendages of barnacles with the Chelex method (Walsh, Metzger & Higuchi, 1991). Trials used the method of Schuelke (2000) to amplify fragments and simultaneously tag forward primers with a fluorescent dye. Loci that amplified and scored consistently in all individuals were fluorescently labeled with 6-FAM, NED or HEX (Applied Biosystems, Custom Oligo Synthesis Center), and used on a larger number of individuals to characterize the microsatellite loci.

Microsatellite loci amplification

Genomic DNA was extracted from feeding appendages of barnacles with the Chelex method (Walsh, Metzger & Higuchi, 1991). PCR amplifications were performed in 20 ul volumes containing final concentrations of 1× PCR buffer (Bioline), 5% bovine serum albumin 10 mg/mL (Sigma), 200 mM each dNTP, 2 nM MgCl, 0.5 mM each primer, 0.5 units of Promega GoTaq DNA Polymerase, and 1 μl template DNA. PCR conditions were as follows: 4 min initial denaturation, followed by 40 cycles with 45 sec denaturing at 94 °C, 60 sec annealing at 55 °C, 60 sec extension at 72 °C and a final extension time of 10 min. The PCR were carried out in a MJ Research PCR Engine. HiDi and ROX500 size standard were added to each sample, and fragment length analysis was carried out at the Georgia Genomics Facility on an ABI 3730xl. Peaks were called and binned with the microsatellite plugin of Geneious version 8.1 (Kearse et al., 2012).

Microsatellite loci characterization

We first characterized the microsatellite loci on 42 individuals from the Atlantic population. We inspected peak calls for fragment size consistency, using the R package MsatAllele (Alberto, 2009). MsatAllele plots peak calls of a locus in histogram form, facilitating visual binning of alleles. If bins could not be clearly assigned, the locus was excluded from the subsequent analysis. We tested whether loci were in Hardy-Weinberg Equilibrium (HWE) by using 1999 Monte Carlo permutations, as implemented by the function hw.test in the R package pegas (Paradis, 2010). Significance values were adjusted for multiple comparisons based on the method of Holm (1979). We recorded the number of alleles, range of fragment sizes, and allelic richness of each locus. The frequency of null alleles was computed based on the method of Brookfield (1996). Genotyping error rates were calculated by repeating genotyping for all individuals. After characterizing the loci in the Atlantic population, we amplified the loci for 24 individuals from Queensland, Australia, to assess if the loci could be used in cross-lineage analysis. Characterization on individuals of the Indo-West Pacific population were largely the same as for the Atlantic population, but we did not assess genotyping error rates. A R script detailing these analyses is available as Supplementary Information.

Microsatellite marker development

The MiSeq run generated 15,324,079 paired-end reads (35–251 bp long) with 81.05% > Q30. Raw reads are available in NCBI’s short read archive (Study accession: PRJNA310774, run accession: SRR3144544). After quality control, 13,498,280 paired-end reads (19–251 bp long) remained, for a total of 6.2 Gb. The median genome coverage was 8× (min = 3, max = 24) for 52 nuclear single-copy gene fragments. The palfinder script detected 629,990 microsatellite repeat regions, of which 29,627 (5.38%) were potentially amplifiable loci (PALs) with forward and reverse primer. A summary of detected microsatellite repeat regions is available as Supplementary Information. A list of all detected microsatellite regions (with primer sequences) is available on figshare (DOI: 10.6084/m9.figshare.2070070). After removing PALs with more than eight or less than two occurrences of either forward or reverse primer in the sequence read data, 17,265 PALs remained. We chose 48 perfect repeat loci for trial amplification, which differed in kmer length and repeat motif, but were otherwise chosen at random. Of those 48 loci, 12 loci amplified and scored consistently throughout the trials, and were tagged with fluorescently labeled dye (Table 1).

Microsatellite marker characterization

We genotyped 42 individuals successfully at more than half of the 12 consistently scoring loci. Visual inspection of peak call histograms revealed that peak calls of Ctest2 did not have clearly defined bins, and were excluded from subsequent analyses. The number of alleles of the 11 scorable loci ranged from six to thirty (Fig. 1). Microsatellite genotype and collecting date for each individual are available as Supplementary Information. For the Atlantic population, four loci were not in HWE, and showed homozygote excess (Table 2). The estimated percentage of null alleles ranged from 0–23%. Allelic richness ranged from 3.65–20.7, and genotyping error rates ranged from 0–7% (Table 2).

All 11 loci amplified in some of the 23 Australian individuals, and had at least two alleles (Figure 2). No loci amplified in all individuals, and no individual failed to amplify all loci. For the Pacific lineage of C. testudinaria, two loci deviated significantly from HWE (Table 3). The locus Ctest7 showed heterozygote excess, even though not at a significant level (p-value = 0.13). Allelic richness of the Australian population was overall lower than in the Atlantic population. The percentage of null alleles ranged from 0–33%.