Skip to main content

Microsatellite marker development from next-generation sequencing in the New England cottontail (Sylvilagus transitionalis) and cross-amplification in the eastern cottontail (S. floridanus)



The New England cottontail (Sylvilagus transitionalis) is a species of high conservation priority in the Northeastern United States, and was a candidate for federal listing under the Endangered Species Act until a recent decision determined that conservation actions were sufficient to preclude listing. The aim of this study was to develop a suite of microsatellite loci to guide future research efforts such as the analysis of population genetic structure, genetic variation, dispersal, and genetic mark-recapture population estimation.


Thirty-five microsatellite markers containing tri- and tetranucleotide sequences were developed from shotgun genomic sequencing of tissue from S. transitionalis, S. obscurus, and S. floridanus. These loci were screened in n = 33 wild S. transitionalis sampled from a population in eastern Massachusetts, USA. Thirty-two of the 35 loci were polymorphic with 2–6 alleles, and observed heterozygosities of 0.06–0.82. All loci conformed to Hardy–Weinberg Equilibrium proportions and there was no evidence of linkage disequilibrium or null alleles. Primers for 33 of the 35 loci amplified DNA extracted from n = 6 eastern cottontail (S. floridanus) samples, of which nine revealed putative species-diagnostic alleles. These loci will provide a useful tool for conservation genetics investigations of S. transitionalis and a potential diagnostic species assay for differentiating sympatric eastern and New England cottontails.


The New England cottontail (Sylvilagus transitionalis) is a narrow niche specialist that relies on the dense understory vegetation of early successional or shrubland habitats [1]. These ephemeral habitats have declined steeply in recent decades in the Northeastern United States [2,3,4]. As a result, many shrubland species, including the New England cottontail, face severe population declines [5]. Remnant New England cottontail populations are found today in < 14% of the species’ historical range and occur in five geographically and genetically distinct populations located in southern Maine and southeastern New Hampshire, central New Hampshire, eastern Massachusetts on Cape Cod, eastern Connecticut and Rhode Island, and western Connecticut and New York [5]. Within these areas, cottontails face the consequences of population isolation and loss of genetic diversity [6, 7].

As a result of uncertainty in long-term species’ viability, the New England cottontail is a high priority for conservation in the Northeastern United States. To help recover the species, state and Federal natural resource managers collaborated with academicians and other stakeholders to develop and implement a conservation strategy to improve the outlook for the species [8]. As a result, the United States Fish and Wildlife Service determined listing of the species under the Endangered Species Act was no longer warranted [9]. Despite this decision, much uncertainty remains about the species’ viability and population status [7, 10]. To evaluate the effectiveness of the conservation strategy and inform adaptive modifications to improve outcomes, extensive population monitoring and research into population structure and genetic diversity are needed [8]. To this end, our goal was to develop a suite of microsatellite markers to aid future conservation genetics research on the New England cottontail as well as other Sylvilagus species.

Main text


For microsatellite marker development, total genomic DNA was extracted from tissue samples obtained from 13 New England cottontail, one eastern cottontail (S. floridanus), and two Appalachian cottontail (S. obscurus) using the DNEasy Blood and Tissue kit (Qiagen, Germantown, MD). We chose to sequence multiple species to maximize the number of microsatellites available for screening in the target species S. transitionalis. Sequencing libraries were created from these genomic DNA samples for sequencing on the Ion Torrent PGM and Ion Proton (ThermoFisher Scientific, Frederick, MD). For each sequencing library, the extracted DNA was quantified with a Nanodrop spectrophotometer (ThermoFisher Scientific), and 100 ng of DNA from each individual was used for construction of multiple 200 base pair libraries following the manufacturer’s protocol. For libraries with multiple individuals in one run, equimolar proportions of each library were pooled prior to chip loading.

Sequence reads from each completed run were imported into Qiagen CLC Genomics Workbench (ver. 7) for processing. Reads were length trimmed to a minimum size of 30 bp, and the quality trimming threshold was set to 0.01 corresponding to a Phred score of 20. Sequence reads from each of the 5 separate runs were then screened individually for all possible di-, tri-, tetra, penta-, and hexanucleotide microsatellite repeat motifs with a minimum repeat length of five with the program QDD (ver. 3) [11]. Primers were designed for each putative microsatellite locus within QDD using the integrated PRIMER 3 code [12]. Seventy-five tri- and tetra nucleotide loci identified by QDD with predicted amplified lengths between 100 and 250 bp were selected for screening within a collection of n = 8 wild-caught New England cottontails sampled from across the species range. A universal M13 sequence was added to the 5′ end of either the forward or reverse primer of each primer pair enabling an M13 FAM labeled fluorescent dye complementary to the universal tail to be incorporated into the polymerase chain reaction (PCR) product [13]. Polymerase chain reactions were performed in 25 μl volumes, consisting of 10 ng of DNA, 1 X PCR Buffer (Promega, Madison, WI), 0.25 μM of labeled forward primer, 0.5 μM of unlabeled reverse primer, 0.1 μM of labeled M13, 2.0 mM MgCl2, 0.2 mM of each dNTP, 0.25 units/μl Bovine Serum Albumin (New England Biolabs, Ipswich, MA), and 0.06 units/μl of Taq polymerase (Promega), using the following cycling conditions: 94 °C for 15 min, 29 cycles of 94 °C for 1 min, 58 °C for 45 s, and 72 °C for 45 s, 5 cycles of 94 °C for 1 min, 52 °C for 45 s, and 72 °C for 45 s, all followed by 72 °C for 10 min. PCR products for each locus were electrophoresed separately on an ABI 3130 Genetic Analyzer (ThermoFisher Scientific) automated DNA sequencer. Alleles were called using GeneMapper (ver. 4) (ThermoFisher Scientific) following the protocols described in [14].

Loci that amplified consistently and were easily scoreable on the sample of n = 8 wild New England cottontails were then tested on DNA extracted from an additional sample of n = 33 New England cottontails from a wild population in eastern Massachusetts, on Cape Cod. In addition, DNA extracted from a sample of n = 6 sympatric eastern cottontails, also collected from eastern Massachusetts were tested for cross-species amplification.

Data analyses

Genotype data from the Cape Cod New England cottontail collection (n = 33) were analyzed for null alleles, large allele dropout, and scoring errors using MICRO-CHECKER (ver 2.2.3) [15]. Genetic diversity and heterozygosity were quantified using GenAlEx 6.5 [16, 17]. Exact tests in GENEPOP [18] were used to determine if genotypes at each locus conformed to Hardy–Weinberg equilibrium (HWE). Multi-locus tests of conformance to HWE were completed using Fisher’s method in GENEPOP. Linkage disequilibrium (LD) was tested for all pairs of loci using contingency tables in GENEPOP. All tests of HWE and LD tests in GENEPOP used the default Markov chain parameters. Significance levels for HWE and LD tests were adjusted using the sequential Bonferroni correction [19].

To evaluate the utility of these loci for population genetic studies, multiple analyses were performed with the n = 33 New England cottontail collection. All multilocus genotypes were subjected to analysis via GENECAP [20] to identify matching samples, calculate match probabilities, and estimate the sibling probability of identity (PI sibs ) [21]. To assess the randomness of the collection (e.g., to ensure the collection did not consist of a small number of families), we analyzed for the presence of full-sibling families using the program COLONY v2.0 [22]. Settings for COLONY analyses included the assumption of male and female polygamy, no per locus genotyping error information, no inbreeding, long run length with full likelihood analysis, high likelihood precision, no allele frequency updates, and no sibship prior. Individual New England cottontails were analyzed as offspring without assignment of individuals as candidate males (fathers) or females (mothers), as these data were not available. While the inference of family relationships is weakened in this situation with no sex, age, relationship information, and the assumption of polygamy for both sexes, COLONY is predicted to be more accurate than pairwise estimates of relationships [22]. Statistical significance of any full-sib relationships was assessed through the P values reported by COLONY. Within a sample of individuals taken at random (with respect to kin) from a population, the frequencies of full and half sib dyads can be used to estimate the current effective size (N e) of the population. Therefore, COLONY was also used to estimate N e and associated 95% confidence interval utilizing the estimates from the sibship assignment full likelihood method. To estimate whether the N e has remained constant (i.e., achieved mutation-drift equilibrium; see [23]), we conducted analyses with BOTTLENECK [24], implementing a two-phased model of mutation (5% IAM; 95% SMM; [25]). Statistical significance of the BOTTLENECK analyses was determined with a Wilcoxon signed-rank test performed by the software.

Results and discussion

Within the 44,084,987 million sequence reads analyzed across five runs of the Ion Torrent PGM and Ion Proton platforms, 1,420,351 million were identified by QDD software as containing microsatellites using our chosen filtering criteria. Of the 75 loci selected from among these sequences for screening in the sample of n = 8 New England cottontails, all but three amplified consistently. An additional 24 loci were monomorphic. Of the remaining 48 loci, 17 had two alleles, 16 had three alleles, 10 had four alleles, three had five alleles, and two had seven alleles. From this set of loci, we selected the 35 most polymorphic, including all with three to seven alleles and a few with two alleles for additional testing in the wild collection of n = 33 New England Cottontails and n = 6 eastern cottontails from the same population to test for cross-specific amplification. The primer sequences of these 35 loci are provided in Table 1.

Table 1 Characteristics of 35 microsatellite loci in n = 33 Sylvilagus transitionalis and cross-amplification in n = 6 S. floridanus

Microchecker reported no instances of large allele dropout or scoring errors. Analyses of family structure in COLONY indicated the presence of no full sibling dyads, and therefore all individuals from the n = 33 collection were retained for subsequent analyses. Three loci were monomorphic in the Massachusetts wild New England cottontail population, while the remaining 32 loci averaged 3.1 alleles per locus (range 2–6; Table 1). Unique multilocus genotypes were generated for each individual with the probability of 6.4 × 10−8 that two siblings would share identical genotypes (PI sibs ; [21]). Expected heterozygosity (He) ranged widely across loci from 6.0% (StrQ23, StrQ48) to 75.0% (StrQ30) and averaged 47.0% for the Massachusetts wild New England cottontail population. Tests for conformance to HWE revealed that locus StrQ19 was not in HWE (P = 0.016), though this result did not remain significant after sequential Bonferroni correction. In addition, no statistically significant linkage disequilibrium (GENEPOP) was detected after sequential Bonferroni correction.

COLONY estimated the N e (and 95% confidence limits) for the Massachusetts collection to be 48 (30–80), and BOTTLENECK indicated a statistically significant heterozygote excess for the Massachusetts collection (P < 0.001; Wilcoxon signed-rank test), suggesting a recent reduction in the effective population size. These results are consistent with previous work, showing that New England cottontails in eastern Massachusetts have low effective population sizes and have undergone a recent population decline [6, 26].

All but two loci amplified in the eastern cottontail samples, and an additional three loci were monomorphic (including two that were also monomorphic in New England cottontails). The remaining 30 loci had two to six alleles in eastern cottontail. Eight polymorphic loci (StrQ23, Str25, Str41, Str42, Str43, Str44, Str48, and Str63) had species-specific amplification patterns with no overlapping alleles across the screened individuals. In additional, StrQ72 was monomorphic for a different allele in each species. These results suggest a subset of these loci will be useful as a diagnostic screening tool for samples collected during noninvasive genetic monitoring surveys [10, 27]. However, more extensive testing of samples throughout the range where New England cottontail and eastern cottontail are sympatric is needed to confirm these results.

In conclusion, the markers developed in this study could prove useful for future population monitoring and conservation genetics research of New England cottontail. The high discriminatory power of these loci indicate that they will be particularly useful for noninvasive genetic surveys, where robust individual identification from a small panel of markers is required. Further, the cross amplification in eastern cottontails will facilitate use of these markers in this species as well as possibly others in the Sylvilagus genus.


These loci have been thoroughly screened in one of the geographically distinct populations of New England cottontails. Given effects of population isolation and genetic drift, it is likely that these markers may have different levels of polymorphism or null alleles in other geographic populations across the species range. Future work should screen the loci carefully in each geographic area prior to embarking on a new study. In addition, while we did not assess the genotyping error rate for these loci, future studies using them could measure this parameter to ensure accuracy of the genotype data.


  1. 1.

    Litvaitis JA. Response of early successional vertebrates to historic changes in land use. Conserv Biol. 1993;7:866–73.

    Article  Google Scholar 

  2. 2.

    Litvaitis JA. Are pre-Columbian conditions relevant baselines for managed forest in the northeastern United States? For Ecol Manag. 2003;185:113–26.

    Article  Google Scholar 

  3. 3.

    Buffum W, McWilliams SR, August PV. A spatial analysis of forest management and its contribution to maintaining the extent of shrubland habitat in southern New England, United States. For Ecol Manag. 2011;262:1775–85.

    Article  Google Scholar 

  4. 4.

    King DI, Schlossberg S. Synthesis of the conservation value of the early successional stage in forests of eastern North America. For Ecol Manag. 2014;324(1):86–195.

    Google Scholar 

  5. 5.

    Litvaitis JA, Tash JP, Litvaitis MK, Marchand MN, Kovach AI, Innes R. A range-wide survey to determine the current distribution of New England cottontails. Wildl Soc Bull. 2006;34:1190–7.

    Article  Google Scholar 

  6. 6.

    Fenderson LE, Kovach AI, Litvaitis JA, Litvaitis MK. Population genetic structure and history of fragmented remnant populations of the New England cottontail (Sylvilagus transitionalis). Conserv Genet. 2011;12:943–58.

    Article  Google Scholar 

  7. 7.

    Fenderson LE, Kovach AI, Litvaitis JA, O’Brien KM, Boland KM, Jakubas WR. A multiscale analysis of gene flow for the New England cottontail, an imperiled habitat specialist in a fragmented landscape. Ecol Evol. 2014;4:1853–75.

    Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Fuller S, Tur A. Conservation strategy for the New England cottontail (Sylvilagus transitionalis). 2012. Accessed 3 April 2017.

  9. 9.

    United States Department of interior (USDOI). US Fish and Wildlife Service 12-month finding on a petition to list the New England Cottontail as an endangered or threatened species. 2015. Accessed 10 January 2017.

  10. 10.

    Brubaker DR, Kovach AI, Ducey MJ, Jakubas W, O’Brien KM. Factors influencing detection in occupancy surveys of a threatened lagomorph. Wildl Soc Bull. 2014;38:513–23.

    Article  Google Scholar 

  11. 11.

    Meglécz E, Costedoat C, Dubut V, Gilles A, Malausa T, Pech N, Martin J-F. QDD: a user-friendly program to select microsatellite markers and design primers from large sequencing projects. Bioinformatics. 2010;26(3):403–4.

    Article  PubMed  Google Scholar 

  12. 12.

    Rozen S, Skaletsky H. Primer 3 on the WWW for general users and for biologist programmers. Methods Mol Biol. 2000;132:365–86.

    CAS  PubMed  Google Scholar 

  13. 13.

    Boutin-Ganache I, Raposo M, Raymond M, Deschepper CF. M13-tailed primers improve the readability and usability of microsatellite analyses performed with two different allele-sizing methods. Biotechniques. 2001;31(1):24–8.

    CAS  PubMed  Google Scholar 

  14. 14.

    King TL, Eackles MS, Chapman DC. Tools for assessing kinship, population structure, phylogeography, and interspecific hybridization in Asian carps invasive to the Mississippi River, USA: isolation and characterization of novel tetranucleotide microsatellite DNA loci in silver carp Hypophthalmichthys molitrix. Conserv Genet Resour. 2011;3:397–401.

    Article  Google Scholar 

  15. 15.

    Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.

    Article  Google Scholar 

  16. 16.

    Peakall R, Smouse PE. GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006;6:288–95.

    Article  Google Scholar 

  17. 17.

    Peakall R, Smouse PE. GENALEX 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics. 2012;28:2537–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Rousset F. GENEPOP’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour. 2008;8:103–6.

    Article  PubMed  Google Scholar 

  19. 19.

    Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.

    Article  PubMed  Google Scholar 

  20. 20.

    Wilberg MJ, Dreher BP. GENECAP: a program for analysis of multilocus genotype data for non-invasive sampling and capture–recapture population estimation. Mol Ecol Notes. 2004;4:783–5.

    Article  Google Scholar 

  21. 21.

    Evett IW, Weir BS. Interpreting DNA evidence: statistical genetics for forensic scientists. Maine: Sinauer Associates, Inc.; 1998.

    Google Scholar 

  22. 22.

    Wang J, Santure AW. Parentage and sibship inference from multilocus genotype data under polygamy. Genetics. 2009;181(4):1579–94.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Davies N, Villablanca FX, Roderick GK. Determining the source of individuals: multilocus genotyping in non-equilibrium population genetics. Trends Ecol Evol. 2001;14:17–21.

    Article  Google Scholar 

  24. 24.

    Piry S, Luikart G, Cornuet J-M. BOTTLENECK: a computer program for detecting reductions in the effective population size using allele frequency data. J Hered. 1999;90:502–3.

    Article  Google Scholar 

  25. 25.

    Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996;144:2001–14.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Papanastassiou K. A landscape genetics approach for comparing connectivity across the range of the New England cottontail. Masters thesis. University of New Hampshire, Natural Resources Department; 2015.

  27. 27.

    Kovach AI, Litvaitis MK, Litvaitis JA. Evaluation of fecal mtDNA analysis as a method to determine the geographic distribution of a rare lagomorph. Wildl Soc Bull. 2003;31:1061–5.

    Google Scholar 

Download references

Authors’ contributions

TLK performed bioinformatics analyses, helped conceive the study, and helped draft the manuscript. He passed away before completion of the final draft. His contributions to the field of conservation genetics are noted, and he will be missed. ME carried out the high throughput sequencing and laboratory screening of the markers and helped draft the manuscript. TJM helped in the development of the study, conducted quality control screening of the markers and helped write the manuscript. AA performed bioinformatic analyses and helped write the manuscript. TPH helped in the development of the study, secure funding for the project, and provided samples. AT helped in the development of the study, assisted with sample collection, and edited the manuscript. AIK helped conceive the study and write the manuscript. All authors read and approved the final manuscript.


Mary Sullivan (University of Rhode Island) helped process New England cottontail samples. We would like to thank the following individuals and their affiliates for providing New England cottontail samples, as the study would not have been possible without their efforts: Eileen McGourty (USFWS), Annie Curtis (MA Army National Guard), David Scarpitti (MA Division of Fish and Wildlife), Howard Kilpatrick (CT Department of Energy and Environmental Protection), Heidi Holman (NH Fish and Game Department), Wally Jakubas (ME Department of Inland Fisheries and Wildlife), Michael McBride and Lou Perrotti (Roger Williams Park Zoo), Jacob McCumber (MA Army National Guard), Tony Perry (Mashpee Wampanoag Tribe). Use of trade, product, or firm names does not imply endorsement by the US Government.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The microsatellite sequences generated from this study are available through NCBI’s GenBank ( and are accessible via the GenBank Accession Numbers KX530819–KX530853.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Tissue samples of wild S. transitionalis were collected under the University of Rhode Island Institutional Animal Care and Use Committee approved protocol number AN11-12-011. One DNA sample from S. obscurus was kindly provided by Marian Litvaitis (University of New Hampshire), and the other S. obscurus tissue sample (specimen 2235-KB) was provided by the Frostburg State University Mammal Museum from a historical sample collected in 1998. Tissue samples of S. floridanus were provided by Anthony Tur of the USFWS, whom obtained the samples from licensed hunters.


Funding for this study was provided by a U.S. Fish and Wildlife Service—U.S. Geological Survey Science Support Partnership Grant.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information



Corresponding author

Correspondence to Aaron Aunins.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

King, T.L., Eackles, M., Aunins, A. et al. Microsatellite marker development from next-generation sequencing in the New England cottontail (Sylvilagus transitionalis) and cross-amplification in the eastern cottontail (S. floridanus). BMC Res Notes 10, 741 (2017).

Download citation


  • Microsatellites
  • Cross-amplification
  • Sylvilagus transitionalis
  • Sylvilagus floridanus
  • Sylvilagus obscurus
  • Next-generation sequencing