Turtle soup, Prohibition, and the population genetic structure of Diamondback Terrapins (Malaclemys terrapin)

Diamondback terrapins (Malaclemys terrapin) were a popular food item in early twentieth century America, and were consumed in soup with sherry. Intense market demand for terrapin meat resulted in population declines, notably along the Atlantic seaboard. Efforts to supply terrapins to markets resulted in translocation events, as individuals were moved about to stock terrapin farms. However, in 1920 the market for turtle soup buckled with the enactment of the eighteenth amendment to the United States’ Constitution—which initiated the prohibition of alcoholic drinks—and many terrapin fisheries dumped their stocks into local waters. We used microsatellite data to show that patterns of genetic diversity along the terrapin’s coastal range are consistent with historical accounts of translocation and cultivation activities. We identified possible instances of human-mediated dispersal by estimating gene flow over historical and contemporary timescales, Bayesian model testing, and bottleneck tests. We recovered six genotypic clusters along the Gulf and Atlantic coasts with varying degrees of admixture, including increased contemporary gene flow from Texas to South Carolina, from North Carolina to Maryland, and from North Carolina to New York. In addition, Bayesian models incorporating translocation events outperformed stepping-stone models. Finally, we were unable to detect population bottlenecks, possibly due to translocation reintroducing genetic diversity into bottlenecked populations. Our data suggest that current patterns of genetic diversity in the terrapin were altered by the demand for turtle soup followed by the enactment of alcohol prohibition. In addition, our study shows that population genetic tools can elucidate metapopulation dynamics in taxa with complex genetic histories impacted by anthropogenic activities.


Introduction
Turtle soup was a popular food item in the United States during late nineteenth and early twentieth centuries. Although many turtle species were consumed, diamondback terrapins PLOS ONE | https://doi.org/10.1371/journal.pone.0181898 August 9, 2017 1 / 18 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 bottlenecks would fail to detect population contractions in populations known to have experienced serious decline, due to human-mediated influxes of genetic diversity.

Sampling and population structure
We analyzed six polymorphic microsatellite loci (TerpSH1, TerpSH2, TerpSH3, TerpSH5, TerpSH7, TerpSH8) [17] from 12 sampling localities throughout the Gulf and Atlantic seaboards and sampled a total of 320 individuals. (Fig 1A). Individuals were sampled using blood or leg muscle tissue, and DNA isolated using a phenol-chloroform extraction. This dataset includes sampling from Texas (TX), Florida (FL), South Carolina (SC), North Carolina (NC), Maryland (MD), New Jersey (NJ), and New York (NY) during 2001 ( Fig 1A). This dataset was previously analyzed for tests of allelic richness, heterozygosity, and number of alleles and can be found in [9]. All sampling protocols and specimen handling followed Ohio University IACUC number L02-06, protocol 13-L-023. We demarcated coastal population structure with two methods. The first was a Bayesian analysis in the program STRUCTURE v. 2.3 [18]. The second was a discriminate analysis of  principal components (DAPC) in the R package adegenet v. 1.4-2 [19]. STRUCTURE delineates genetic clusters by maximizing conformity to Hardy-Weinberg equilibrium while simultaneously minimizing linkage disequilibrium among loci for K user-defined clusters, while DAPC utilizes k-means to maximize variation between groups after PCA transformation. For STRUCTURE, we ran K from 1 to 12 (sampling locations), replicating each value of K ten times, each with a random starting seed. Individual runs were composed of a Markov chain Monte Carlo (MCMC) algorithm of 700,000 steps, with the first 50% removed as burn-in. We used the admixture model, the LOCprior, and LOCISPOP prior for all runs. Preferred values of ΔK were computed with the Evanno method [20] in STRUCTURE HARVESTER v. 0.6.94 [21]. Output from STRUCTURE HARVESTER was processed in CLUMPP v. 1.1.2 [22] to control for labeling switching and multimodality. DISTRUCT v. 1.0 [23] was used to visualize the data. We ran STRUCTURE in a hierarchical fashion: first we ran an initial analysis to detect basal levels of structure [20], and then we searched for additional structure within each identified population individually. For DAPC, we first optimized the number of PC axes to avoid under-or over-fitting our population genetic models. We accomplished this with cross-validation, which uses stratified random sampling and divides the dataset into a training set and a validation set. We partitioned the training set to be 50% of the data and the validation set to 50%, and employed 100 replicates. The cross-validation method estimated 60 axes should be retained. The Bayesian information criterion (BIC) was used to choose the optimum number of clusters for DAPC analysis. Four discriminate functions (DF) were retained for each analysis.

Historical and contemporary gene flow
We used MIGRATE v. 3.6.5 [24] to estimate historical gene flow levels (M = m h /μ: proportion of migrants per generation, scaled by mutation rate). For these analyses we treated sampling localities as populations with the exception of South Carolina. Because South Carolina lacked structure among its sampling localities (Fig 1G; see results), we treated all South Carolina sampling localities as a single population. For each population, we randomly subsampled 40 individuals for gene flow estimates. If a population had fewer than 40 individuals, we included all individuals in gene flow estimates. For each run, we slice sampled the posterior distribution with a MCMC of 5,000,000 steps with the first 50% removed as burn-in. Each MCMC consisted of four statically heated parallel chains sampled every 500 iterations. Five independent replicates were run, totaling 25,000,000 MCMC steps. Estimates of θ (= 4N e μ: effective population size, scaled by mutation rate) were modeled under a uniform distribution bounded between 0.0001 and 100, and historical gene flow estimates were bounded between 0 and 2000.
For estimates of contemporary gene flow (m: proportion of migrants per generation), we used BAYESASS v. 3.0 [25]. We ran 10 independent MCMC simulations with random starting seeds [26] for 30,000,000 steps, and sampled every 3,000 steps. We discarded the first 50% as burn-in. We used TRACER v. 1.5 [27] to visualize MCMC simulations, and used R scripts to calculate a Bayesian deviancy measure to determine the run that best fit the data [28,29].

Temporal changes in gene flow
Because MIGRATE uses the coalescent, it estimates gene flow over long periods of time, approximately 4N e generations (several thousand years) in the past [24]. Thus, its gene flow estimates pre-date translocation events in the early twentieth century. By contrast, BAYESASS estimates gene flow ". . .over the last several generations" [25], where "several" is commonly interpreted as roughly five generations (e.g. [12,30] With a generation time of 12 years (Roosenburg, unpublished data), the contemporary gene flow estimates of BAYESASS are within the last 60 years or so. These estimates post-date known terrapin translocation events.
To estimate temporal changes in gene flow, we compared historical and contemporary estimates. First, we multiplied the historical gene flow estimates from MIGRATE (M = m h /μ) by a mutation rate (μ) of 2.72 x 10 −3 mutations site -1 generation -1 . This mutation rate was estimated explicitly for the microsatellites used in previous studies [31].

Testing coastal population structure
We compared eight models of gene flow in MIGRATE, generated with the methods described above (Fig 2). Our initial model was a linear stepping-stone (Model A), which restricted gene flow to adjacent populations. We then added or removed gene flow routes from Model A based on results from STRUCTURE, DAPC, and our Δm estimates. Models B-F incorporated gene flow from translocation events with varying connectivity to Florida. Models G and H modeled the Suwannee Seaway. We used approximate Bézier scores and log Bayes factors (LBF) to determine which model best explained coastal population structure.

Bottlenecks
BOTTLENECK v. 1.2.02 [32] was used to test for bottlenecks. Locus mutation was modeled with the stepwise-mutation model (SMM) and the two-phase model (TPM). The TPM modeled 95% of mutations as single-step while multi-step mutations were modeled with 12% variance [32]. Each test consisted of 20,000 permutations. We tested for bottlenecks by sampling locality and by STRUCTURE cluster. We used the Wilcoxon signed-rank test, which detects bottlenecks 25-250 generations in the past [33], and a mode-shift test, which detects bottlenecks ". . .within the past few dozen generations" [34].

Population structure-STRUCTURE
Initial runs of STRUCTURE found ΔK = 2 best describes coastal population structure ( Fig  1B). However, ΔK scores of 4 and 7 are comparable, and we also report them for comparative purposes (Fig 1C and 1D). As ΔK identifies basal levels of hierarchical structure [20], we also tested for substructure within groups for ΔK = 2, which included northern and southern groups. Inference of population substructure within the northern group (NC, MD, NJ, NY) revealed K = 2 or K = 3 subclusters (Fig 1E and 1F). This is due to two solutions of K having overlapping likelihood scores (S1 Fig). If K = 2 is adopted, NJ and NY constitute a cluster ( Fig  1E), while if K = 3 is adopted, they form separate subclusters ( Fig 1F). The southern group (TX, FL, SC) included ΔK = 3 subclusters (Fig 1G), one for each state. We did not detect substructure among SC sampling localities. Thus, in total there are either five or six genotypic clusters, dependent upon K = 2 or K = 3 for north Atlantic terrapins (S1 Fig). Population structure-DAPC BIC indicated the optimum number of clusters is six (Fig 3). DAPC for the six-cluster analysis is shown in Fig 4, and membership probabilities are provided in Fig 5. Mid-and north Atlantic terrapins form clusters 1, 3, and 6; cluster 4 diagnoses Chesapeake Bay (MD) terrapins, and cluster 2 includes populations located in the Gulf (TX, FL). Cluster 5 is composed of SC terrapins and links cluster 4 with clusters 1, 3, and 6. Cluster 2 exhibits no overlap with any cluster.

Historical and contemporary gene flow
The highest levels of historical gene flow were from NJ to NC (Table 1;  Of the 42 gene flow routes, six exhibited increased levels of contemporary gene flow (Table 3; +Δm > 0.010), 22 demonstrated reduced levels of contemporary gene flow (-Δm < -0.010), and 14 were relatively stable over time (-0.010 < Δm < 0.010). Four of the six gene flow routes estimated to have increased contemporary gene flow are found between adjacent populations. The two routes between non-adjacent populations to show higher levels of contemporary gene flow are from TX to SC (Δm = +0.0567) and from NC to NY (Δm = +0.0320).

Models of population structure
As the diamondback terrapin has a linear distribution (Fig 1), we used a stepping stone process as our null model of population connectivity (Fig 2A). To model translocation events, we modified the stepping stone model to include unidirectional gene flow between TX and SC, and NC and NY; these models differed only in their relative isolation of FL (Fig 2B-2F). We modeled the Suwannee Seaway by modeling bidirectional gene flow between TX and SC, with FL either completely isolated (Fig 2G) or a sink population (Fig 2H). We found that the Atlantic  Exchange model (Fig 2E), which modeled translocation events from TX to SC and NC to NY (gene flow between non-adjacent populations), and which included bidirectional gene flow between SC and FL but unidirectional gene flow between TX and FL, vastly outperformed all other models, including the linear stepping-stone model (Bézier scores = -13,7232 vs. -72,456.30). Converting approximate Bézier scores into posterior probabilities confirmed the Atlantic Exchange model as the best-fit model (PP = 100%). All other models had estimated posterior probabilities near 0%.

Bottlenecks
We recovered no evidence for genetic bottlenecks for any sampling locality with the TPM or SMM, although the TPM for SC6 approached significance (Table 4; P = 0.055). Mode-shift tests also failed to detect bottlenecks for any sampling locality. In addition, neither the TPM, SMM, nor a mode-shift test detected bottlenecks for any STRUCTURE cluster (Table 4).

Discussion
Over the last two centuries, the relationship between humans and terrapins has been complex, and the terrapin's population genetic structure reflects this relationship. The demand for turtle soup resulted in historical population contractions and extirpations, and culminated in the construction of terrapin farms [3][4][5][6]. To get flavorful terrapins to market quickly, Texas and Carolina terrapins were hybridized at the North Carolina terrapin farm [6]. Then, in 1920, the enactment of Prohibition restricted access to sherry, which drastically cut demand for turtle soup. Consequently, many terrapins were released into local waters, which promoted population admixture and may have resulted in the reintroduction of genetic diversity. We documented population genetic structure that is consistent with historical accounts of terrapin translocation during the twentieth century [1,3]. We recovered two or three genotypic clusters in the mid-and north Atlantic (MD, NJ-NY) and three genotypic clusters in the    Gulf and southern populations (TX, FL SC), for a total of six genotypic clusters (Fig 1B-1G;  Fig 3). We also recovered the North American Gulf/Atlantic phylogeographic divide previously described in the terrapin and other taxa ( Fig 1C) [8,31,[35][36][37]. Recent work [38] has suggested that sampling localities with uneven sampling may lead to erroneous STRUCTURE results. While this is a concern, we do not think it is affecting our data. The population of concern in our study (Texas) has been recovered in previous studies [9], and is supported by our DAPC results (Figs 6 and 7) and an additional PCA analysis (S2 Fig).
In addition to delineating population structure and quantifying population connectivity, we also found that a modified stepping-stone model with genetic exchange along the Atlantic seaboard and unidirectional gene flow from TX to SC and from NC to NY best describes terrapin population structure (Fig 2E). This model of population connectivity outperformed a linear stepping-stone model, as well as models of the Suwannee Seaway, a natural conduit of gene flow between Gulf and Atlantic populations (Fig 2). Furthermore, we found Florida populations to be divergent from neighboring populations (Fig 1B-1D and 1G; Figs 4 and 6; Table 3), which complements accounts that Florida terrapins were not translocated due to their inferior size and taste [3], as well as other studies that reported FL terrapins to be genetically distinct [9][10][11].
Our Bayesian model comparisons supported bidirectional gene flow between SC and FL, and unidirectional gene flow from TX to FL (the Atlantic Exchange model, Fig 2E), but did not support alternative models of connectivity, such as the Suwannee Seaway (Fig 2).
We found that NC is highly admixed with other mid-north Atlantic populations (Fig 1C-1F), which could be the result of human-mediated gene flow. For example, relative to historical levels of gene flow, we documented increased contemporary connectivity from NC to NY Table 2    Human-mediated gene flow in terrapins (Table 3), which is consistent with accounts of translocation [1]. North Carolina was the location of a terrapin breeding operation [4][5][6], while NY was the location of a large terrapin market [1]. With the volume of terrapins brought to market, it is possible some terrapins escaped or were released into local waters. North Carolina also exhibited increased contemporary gene flow into SC and Chesapeake Bay, MD (Table 3), consistent with a previous study that detected large increases of contemporary gene flow into Chesapeake Bay [12]. Finally, we observed admixture between TX and SC populations (Fig 1D and 1G; Fig 6), as well as increased levels of contemporary gene flow from TX into SC (Table 3), consistent with chronicled translocation and hybridization experiments [1,6,8].
Although we documented population genetic evidence of translocation between some nonadjacent populations, our study failed to find genetic evidence for some known instances of translocation. In particular, we did not detect increases in contemporary gene flow or admixture between TX and NC (Fig 1B-1D; Figs 4 and 6; Table 3). There are several possible explanations: translocation can fail [39,40], and released terrapins from TX may not have successfully interbred with local populations in NC. Alternatively, our sampling in NC may have not included admixed localities (Fig 1A). We also did not find increased contemporary gene flow from MD to NY, although we detected admixture between these populations (Fig 1B-1F; Figs 4 and 6; Table 3). Terrapins from MD may not have been released into local waters in NY, or they may have failed to interbreed. Despite current weak demand for turtle soup, terrapins from the Chesapeake Bay region continue to be sold at markets in New York, often illegally [41,42]. Although our Bayesian model tests cast doubt on the Suwannee Seaway (Fig 2) as a viaduct for genetic transmission in the southern portion of the terrapin's range, it remains possible that long distance (coastal) dispersal and gene flow are structuring a portion of our populations. Despite a variety of potential biogeographic barriers, the Balkan Pond Turtle Human-mediated gene flow in terrapins (Mauremys rivulata) exhibits a paucity of genetic variation across its range [43], suggesting it is capable of transoceanic gene flow. Infrequent but long-range coastal gene flow may also be occurring between Diamondback Terrapin populations. It is well documented that terrapin populations historically underwent severe contractions [2], but we failed to detect any population bottlenecks in any region (Table 4). For example, the decline of terrapins from Chesapeake Bay is documented by shrinking harvests after decades of overexploitation [7], but we did not detect bottlenecks in this region. Indeed, previous work in Chesapeake Bay indicated terrapins exhibit relatively high levels of genetic diversity [9,12]. It is possible that we failed to detect population bottlenecks because tests of heterozygosity excess have weak statistical power [44]. Alternatively, the failure to detect bottlenecks may be the consequence of terrapin translocation events reintroducing genetic diversity into populations. If true, the enactment of Prohibition may have inadvertently benefited the terrapin in two ways. The first was the collapse of the turtle soup market, which slowed the harvesting of natural populations. The second was the closure of terrapin farms and the release of translocated individuals into local populations, which may have reintroduced genetic diversity and increased population viability [45].
Thus, our study suggests that population genetic structure in the diamondback terrapin possesses the signature of historical translocation events. Translocation among natural populations is known to increase levels of population admixture and genetic diversity [46,47]. At least in some cases, increased levels of genetic diversity save populations from the negative consequences of inbreeding depression and lowered mean population fitness [48][49][50][51][52]. However, translocated populations may exhibit high levels of genetic diversity but have low effective population sizes, suggesting several population genetic metrics are required to judge the efficacy of translocation [53]. An alternative outcome of translocation is that it may harm populations by causing outbreeding depression, harm locally adapted populations by moving them away from an adaptive peak [54,55], or introduce diseases [56]. Nonetheless, given the extent of the current biodiversity crisis [57][58][59][60], including increasing population fragmentation [61][62][63][64] conservation-oriented translocation has become increasingly pivotal to maintaining population viability [65].
Historically, terrapins have been divided into seven subspecies [8], but recent genetic studies are incongruent with extant taxonomy [9][10][11]. While our analyses also cast doubt on current terrapin taxonomy, we do not make recommendations for future taxonomic changes. Our study lacked sampling in the Gulf of Mexico, where there may be more population structure than currently documented. Our study also used six microsatellite loci; more thorough geographical sampling and additional loci are required to make recommendations on taxonomic changes.
Our study shows that convoluted genetic histories can be disentangled with modern population genetic tools and that translocation can leave an indelible fingerprint in populations. Anthropogenic influences increasingly disrupt population dynamics [12,61,62,66]; however, the indirect consequences of social and political activities are not always predictable. Our study suggests the population genetic structure in the Diamondback Terrapin may be the byproduct of an interaction between market demand for turtle soup during the late nineteenth and early twentieth centuries, followed by the enactment of Prohibition in 1920, which resulted in the large scale release of captive terrapins into local waters.