Genetic variation across trophic levels: A test of the correlation between population size and genetic diversity in sympatric desert lizards.

Understanding the causes of genetic variation in real populations has been elusive. Competing theories claim that neutral vs. selective processes have a greater influence on the genetic variation within a population. A key difference among theories is the relationship between population size and genetic diversity. Our study tests this empirically by sampling two species of herbivorous lizards (Dipsosaurus dorsalis and Sauromalus ater) and two species of carnivorous lizards (Crotaphytus bicinctores and Gambelia wislizenii) that vary in population size at the same locality, and comparing metrics of genetic diversity. Contrary to neutral expectations, results from four independent loci showed levels of diversity were usually higher for species with smaller population sizes. This suggests that selective processes may be having an important impact on intraspecific diversity in this reptile community, although tests showed little evidence for selection on the loci sequenced for this study. It is also possible that idiosyncratic histories of the focal species may be overriding predictions from simple neutral models. If future studies show that lack of correlation between population size and genetic diversity is common, methods using genetic diversity to estimate population parameters like population size or time to common ancestor should be used with caution, as these estimates are based on neutral theory predictions.

Introduction

Genetic diversity is essential to the process of evolution, however we are still trying to understand the relative importance of factors that influence genetic diversity in natural populations. At the level of DNA sequence diversity, the Neutral Theory predicts that most molecular genetic variation found in populations is selectively neutral, and is shaped primarily through the random processes of mutation and genetic drift [1]. An important prediction of this theory is that the size of a population will have a large impact on its genetic diversity. Genetic drift has a greater effect on smaller populations, leading to a predicted positive correlation between population size and within-population genetic diversity. Empirical evidence from natural populations has been mixed. In some studies, the neutral theory appears to adequately explain observed patterns of genetic diversity [2-4], see [5] for review, while others have found results which do not fit neutral predictions [6-7], see [8] for review, or in which distantly related species with vastly different effective population sizes have a lower than expected difference in genetic diversity [6]. Natural selection is thought to be the most likely explanation for the deviations from neutrality found in [6-8], especially for genetic regions with low recombination.

Alternative theories suggest that natural selection may play a larger role in determining levels of polymorphism through both background (purifying) selection and genetic hitchhiking. Theories such as the Nearly Neutral Theory [9] and Genetic Draft [10, 11] have been proposed as mechanisms by which purifying or balancing selection may shape patterns of genetic variation in the genome. Because natural selection operates more efficiently in large populations than in small ones these theories predict that, except in very small populations, there will be little or no correlation between population size and within-population genetic diversity. This is true not only for selected loci but also for neutral loci in linkage disequilibrium with selected loci. Therefore, selection is expected to have more extensive effects in mitochondrial DNA, haploid organisms, and regions of low recombination [11, 12], since recombination can separate a mutation with a high selection coefficient from its original allelic background. Genome wide polymorphism data appears to support this theory, suggesting that the effects of selection are greater in species with larger population sizes [13]. The accuracy of the statistical methods used to identify regions affected by selection in the genome remains controversial [5, 14].

Despite the amount of attention the role of population size in population genetics has received, no theory has yet fully explained the patterns of genetic variation seen in natural populations [14-16]. Surprisingly few studies have been conducted on natural populations to test the relationship between population size and genetic diversity. The studies that exist often draw conclusions based on samples taken across a species’ range rather than from a single deme [13, 17, 18] thereby confounding within- and between-population effects of evolutionary forces like drift and selection as well as increasing the likelihood that populations have experienced different historical effects. Differences in conditions across the range of a species could also cause large differences in selective pressures. Clearly these factors make it difficult or impossible to untangle patterns of genetic diversity when samples are taken over large geographic distances.

An ideal method for testing neutral expectations is to compare the genetic diversity of taxonomically related species living sympatrically. By focusing on a taxonomically related group, it is possible to reduce many confounding variables that may differ among taxonomically more diverse contrasts, such as potential differences in DNA mutation rates and constraints on locus function. Studying a single area reduces (but does not eliminate) the possibility that the gene pools of the focal species have been responding to radically different events in recent geologic history, or that patterns caused by population subdivision are being conflated with patterns within a deme.

Our research compares genetic diversity in species of lizards that differ in population size. Our previous studies found interspecific differences in levels of genetic diversity in sympatric populations of lizards in the Mojave Desert [19, 20]. Hague and Routman [20] showed that species with smaller population sizes had generally lower genetic diversities than species with larger population sizes. We surveyed genetic variation in two additional species with smaller population sizes from the same locality, and in this paper we compare their DNA sequence diversity to our previous results and expand on the previous datasets. For four focal lizard species we tested the neutral theory prediction of the positive relationship between genetic diversity and population size using three criteria to rank relative population size: (a) trophic level, (b) habitat specialization, and (c) a combination of published abundance surveys and long term observation of the lizard community at our study site. The above criteria predict that population size, and therefore genetic diversity would rank as follows: carnivore specialist< carnivore generalist< herbivore specialist< herbivore generalist.

Methods

Ethics statement

All procedures involving animals in this study followed ethical and legal guidelines. The research protocol was approved by the San Francisco State University Institutional Animal Care and Use Committee (IACUC; animal protocol #A14-06). All collection took place on federally owned land, under National Park Service Scientific Research and Collecting Permit # MOJA-2014/5-SCI-0002 Study #00270 issued to EJR. None of the four focal species are federally listed as threatened or endangered, or listed as species of special concern in the state of California.

Taxa sampled

The lizard communities of the Mojave Desert of North America make an ideal study system for a study of comparative population genetics. Numerous species are found living sympatrically, occupying many niches within the environment. Many species reach high local abundances while others are much less dense. This study focuses mainly on two members of the family Crotaphytidae (Crotaphytus bicinctores, the Great Basin Collared Lizard, and Gambelia wislizenii, the Long-nosed Leopard Lizard), and two members of the family Iguanidae (Dipsosaurus dorsalis, the Desert Iguana, and Sauromalus ater, the Common Chuckwalla) in the Mojave National Preserve, San Bernardino County, California, USA. While relationships among the pleurodont reptiles are not fully resolved, it is clear that Crotaphytidae and Iguanidae are closely related to one another [21, 22]. These four species are all large-bodied lizards [Snout-vent lengths: C. bicinctores—8.6–11.2 cm, G. wislizenii—8.2-14-6 cm, D. dorsalis—10.1–14.6 cm, S. ater—12.7–22.8 cm] with long generation times, essentially controlling for factors other than trophic level and habitat niche. The two crotaphytids are predators, and often prey on smaller lizards [23, 24], while the two iguanids are primarily herbivorous [25, 26]. The general expectation that carnivorous species will be less common than herbivorous species of similar size has been supported in this instance by previous observations of the Mojave lizard community [27, 28]; authors’ observations). By sampling two species from each trophic level we have added a degree of replication to the study design, although it is confounded with phylogeny in this case. G. wislizenii and D. dorsalis are desert habitat generalists, while C. bicinctores and S. ater are found only among rocks. Because each lizard family/trophic level group contains one species in each category, we do not expect habitat niche to be a main determinant of genetic diversity differences between trophic levels. In addition, these four species are not known to differ in characteristics that should affect the relationship between Ne and actual population size, such as tendency to inbreed, and are relatively well matched in body size and generation time.

Some genetic work has already been conducted on these, or closely related species. Phylogeographic studies found very high genetic diversity across the range of Sauromalus obesus (= ater) [29], and S. obesus and D. dorsalis [30], which was likely attributable to population subdivision. A phylogeographic study of G. wislizenii in the Mojave Desert found high haplotype diversity in a single mitochondrial gene, but this appeared to be based on only a few individuals collected from three distant localities [31]. A comprehensive phylogenetic study of the family Crotaphytidae sampled 408 individuals, with C. bicinctores and G. wislizenii heavily represented [32], but made no attempt to quantify intraspecific diversity. Earlier work at this sitefound that the two large iguanids studied here had less genetic diversity than small-bodied species with very high population densities [19,20].

Population size

One difficulty in testing neutral predictions is that it is the effective population size that is expected to affect diversity, rather than the current census population size. Effective population size (Ne) is defined as the number of idealized (usually Wright–Fisher model) individuals needed to account for the observed levels of genetic diversity in the actual population. The relationship between Ne and the actual population size depends on many factors, only one of which is the actual number of individuals, and in practice Ne is usually estimated from empirical estimates of genetic diversity (θ) and the mutation rate (μ) (Ne = θ/4μ for autosomal genes). However, this equation assumes neutrality and cannot be used in studies testing whether genetic diversity and population size are correlated.

In this study, we use a combination of trophic level expectation and long-term observation of relative abundance as a proxy for relative Ne. Ecological theory predicts that, all else being equal, herbivores will be more abundant than carnivores of equal body size and that this relationship is likely to be stable over time[33,34]. Population census [27,28,35] and over 20 years of field work at the study site by one of us (Routman) supports the proposition that these carnivorous species are much less abundant than the herbivorous species. Indeed, for the carnivorous species it took over four field seasons to collect samples that were similar in size to those of the herbivorous species, which were collected in one to two field seasons. This difference in time to complete the sample is an underestimate of the relative population sizes, because most of the samples of the herbivorous species were collected at the start of our work on genetic diversity in the Mojave National Preserve. At that time we were simultaneously searching for 16 different species found in several different habitats, whereas the collection of most of the carnivorous species were later in our overall study, when sampling of most of the common species was completed. Thus, in the later part of the study we were focusing our efforts mainly on the carnivorous species.

Support for the idea that, within trophic level, the rock specialists should have lower population sizes than those of habitat generalists because of less available habitat can be inferred from Fig 1. The fact that rocky outcrops cover much less area than the intervening sandy regions, combined with the fact that the habitat generalist also use the rocky habitats (albeit in lower density) essentially insures that, within trophic level, the habitat generalist will have a larger population size.

Fig 1

Collection sites in the Mojave National Preserve for four species of lizard.

Markers indicate individual lizard collection localities, separated by species as follows: Yellow: Crotaphytus bicinctores; Blue: Gambelia wislizenii; Green: Dipsosaurus dorsalis; Red: Sauromalus ater. Because of the map scale, there is considerable marker overlap, e.g. C. bicinctores markers at the Lava site are not visible.

Some survey data from the Mojave National Preserv supports our ranking of population sizes in these species. Data from a can trap grid [35], collected monthly June 1991-May 1993, January 2000-December 2001, and January 2008-June 2018, yielded the following capture numbers for the 4 focal species: D. dorsalis—68, G. wislizenii– 4, C. bicinctores– 2, S. ater—3. These number reflect our rankings with the exception of S. ater, which has a capture frequency closer to that of the two carnivores. However, S. ater are rarely caught in can trap studies, because their preferred habitat is large boulder fields and rocky cliffs, where can traps cannot be used because the cans require burying.

Collection and sequencing

Collection took place primarily in the Cima Volcanic Field (referred to as Lava below) of the Mojave National Preserve, CA. (Fig 1). This locality was sampled in our previous studies [19,20], allowing for an accurate comparison between datasets. Due to the low abundance of the crotaphytid species, we collected some individuals from nearby desert areas. All of the S. ater samples(n = 37), along with 36 D. dorsalis, 11 C. bicinctores, and 27 G. wislizenii were collected at the Cima Volcanic Field. In order to increase sample sizes, two nearby areas were used as secondary collection sites. 21 C. bicinctores were collected from the nearby Zzyzx Road (23 km from the Cima Volcanic Field). Because the Zzyzx Rd. site has a low density of G. wislizenii, we also collected at the area near Kelso Dunes (35 km from the Cima Volcanic Field), which yielded 5 G. wislizenii and 4 D. dorsalis.

Collection took place over four years (2012–2015). Iguanid samples previously collected and sequenced by Hague and Routman [20] were used for this study, and supplemented by additional individuals of both species to increase sample size. Exact locality data for each specimen can be found in S1 Table. We collected lizards of the target species using slip-knot nooses, and sampled a small piece of tissue (0.5 cm) from the tail of each individual, which was then released. Samples were preserved in 95% ethanol in the field, and subsequently frozen.

We sequenced protein coding regions of mitochondrial cytochrome b (cytb) and three autosomal genes: melanocortin 1 receptor gene (MC1R), recombination activating gene 1 (RAG1), and caspase recruitment domain gene 4 (CARD4). Cytb is commonly used for phylogenetic studies of mitochondrial DNA in vertebrates, and because the mitochondrial genome is non-recombining, should be more affected by genetic hitchhiking than the autosomal genes. We sequenced approximately the same segment of each gene for each of the four species to control for possible differences in variability among regions of the gene.

We used standard tissue extraction methods to prepare DNA from the samples for sequencing (Quick-gDNA™ MiniPrep kit; Zymo Research). PCR was carried out in Accupower PyroHotStart Taq PCR Premix tubes (Bioneer, Inc.) to amplify the selected genes (See S1 File for primers and PCR protocols used). After purification using ExoSAP-IT® (Affymetrix), we sent the resulting PCR product to Elim Biopharmaceuticals (Hayward, CA) for Sanger sequencing in both directions. PCR primers were also used for sequencing except for a few locus-species combinations for which we had to design sequencing primers. We aligned sequence reads and evaluated sequences visually for heterozygosity and accuracy of sequence data using Geneious version 7.1.5 (http://www.geneious.com) [36].

Data analysis

We sequenced a segment of mitochondrial cytb ranging in size from 610–872 base pairs (Table 1). The ends of the alignment for each species were manually trimmed to include only base pairs present in ≥95% of the aligned individuals. Only the ends of some sequences had uncalled bases, involving < 10 bases. Because Arlequin can give nonsensical results for some analyses when missing bases are present at the ends of sequences, we manually replaced any uncalled base with the base present in that position in other individuals having the same sequence for the scored bases (according to its documentation, this is the default assumption built into Arlequin when missing data is incorporated into the analysis). There were no cases in which polymorphism in these regions caused ambiguity regarding which base should be inserted.

Table 1

Population genetics summary statistics, part 1.

Species	Gene	N	# BP	# Variable Sites	# Hap.	Kn:Ks	Kn/Ks ratio
Crotaphytus bicinctores	cytb	32	872	37	7	7:30	0.23
	CARD4	64	799	6	7	2:4	0.50
	MC1R	62	642	16	24	5:11	0.45
	RAG1	64	774	10	11	6:4	1.50
Gambelia wislizenii	cytb	33	610	20	13	3:17	0.18
	CARD4	66	787	15	13	5:10	0.50
	MC1R	66	741	17	22	4:13	0.31
	RAG1	66	804	20	19	11:9	1.22
Dipsosaurus dorsalis	cytb	40	844	1	2	1:0	1.00
	CARD4	80	789	3	4	1:2	0.05
	MC1R	80	755	3	4	0:3	0.00
	RAG1	80	1028	11	13	9:2	4.50
Sauromalus ater	cytb	37	848	13	5	4:9	0.44
	CARD4	74	807	6	8	3:3	1.00
	MC1R	74	813	7	6	1:6	0.17
	RAG1	74	1049	9	7	3:6	0.50

Abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1; N, number of chromosomes sampled; # BP, base pairs; # Hap., number of unique haplotypes; Syn. Mut., synonymous mutations; Nonsyn. Mut., nonsynonymous mutations.

We used the programs PHASE [37,38] and seqPHASE [39] to reconstruct haplotypes from the autosomal genotype data. Haplotypes determined by cloning were available for some heterozygous individuals of D. dorsalis (RAG1 and MC1R) and S. ater (MC1R) [20], and these were used as known haplotypes for PHASE analysis of these species. The best supported haplotype pair was chosen for each individual from the most highly resolved PHASE runs and used for further analysis. As it was not possible to find a best haplotype pair with 90% or greater confidence for each individual, we tested the potential impact of choosing incorrect pairs on our final results. We generated simulated datasets in R [40] using a script that randomly chooses one of the possible haplotype pairs from the “.out_pairs” file for each individual, without taking into account the probability of each haplotype pair being correct (S2 File). A set of 1000 random datasets were generated for each locus-species combination. Batch analyses of these simulated datasets were compared with the results from the analysis using the best supported haplotype pairs.

We used Arlequin (version 3.5.1.2) [41] for the analysis of population level diversity. For each species, we calculated standard measures of population genetic diversity, including haplotype diversity (h), Nei’s θ (θπ, or π), and Watterson’s θ (θS) [42,43]. For all loci, pairwise differences between species in genetic diversity measures were considered significant if the 95% confidence intervals of the two estimates did not overlap by more than half of a one-sided error bar [44]. We used Tajima’s D to test for deviation from a model assuming neutrality and constant population size [43]. Results were checked in DnaSP version 5.10.01 [45], due to some miscalculations we found in the Arlequin results (see below). Pairwise mismatch distributions were also calculated in Arlequin to test for population expansion. We calculated the ratio of nonsynonymous to synonymous changes for each locus-species combination.

Due to the low observed population sizes of some species, we collected individuals over a wider geographic area than had been used for more abundant sympatric species in order to get a sufficiently large sample size. To account for between-location genetic variation, we tested for population structure within the sampling area. Individuals collected from the same localities were assigned to putative separate subpopulations. We calculated FST ([46]; eqn. 9) using estimators of heterozygosity adjusted for sample size [47], and used Arlequin to find ΦST values between the subpopulations [41]. Due to the high levels of diversity in these lizards, we compared FST and ΦST with Jost’s D [48], which estimates the level of differentiation among subpopulations independent of genetic diversity within populations. We calculated Jost’s D in R, using eqn. 12 from [48] (S3 File).

Results

Sample sizes were as follows: 32 C. bicinctores, 40 D. dorsalis, 33 G. wislizenii, and 37 S. ater (= 64–80 copies of each autosomal locus). Grouped by trophic level, 65 individuals of carnivorous species and 77 individuals of herbivorous species were included in the sample. Population genetics descriptive statistics are reported in Tables 1 and 2 and the rank order of diversity for all measures is found in Table 3.

Table 2

Population genetics summary statistics, part 2.

Species	Gene	N	Haplotype Diversity	θS	θπ
Crotaphytus bicinctores	cytb	32	0.750 (±0.131)	10.536 (±7.004)	8.046 (±8.134)
	CARD4	64	0.487 (±0.139)	1.588 (±1.472)	0.877 (±1.460)
	MC1R	62	0.936 (±0.029)	5.307 (±3.706)	4.491 (±5.199)
	RAG1	64	0.703 (±0.108)	2.733 (±2.150)	1.337 (±1.965)
Gambelia wislizenii	cytb	33	0.866 (±0.076)	8.079 (±5.813)	7.861 (±8.589)
	CARD4	66	0.721 (±0.098)	4.005 (±2.822)	3.139 (±3.733)
	MC1R	66	0.922 (±0.031)	4.821 (±3.301)	3.792 (±4.407)
	RAG1	66	0.909 (±0.031)	5.226 (±3.454)	2.699 (±3.293)
Dipsosaurus dorsalis	cytb	40	0.450 (±0.106)	0.278 (±0.546)	0.533 (±1.052)
	CARD4	80	0.561 (±0.094)	0.768 (±0.922)	1.052 (±1.649)
	MC1R	80	0.489 (±0.094)	0.803 (±0.963)	0.750 (±1.340)
	RAG1	80	0.769 (±0.065)	2.161 (±1.626)	3.291 (±3.699)
Sauromalus ater	cytb	37	0.590 (±0.086)	3.672 (±2.836)	3.517 (±4.096)
	CARD4	74	0.502 (±0.125)	1.525 (±1.404)	1.170 (±1.763)
	MC1R	74	0.646 (±0.071)	1.766 (±1.541)	1.066 (±1.649)
	RAG1	74	0.737 (±0.055)	1.760 (±1.414)	1.812 (±2.270)

95% confidence interval in parentheses. All θ values are per base and multiplied by 1000 for clarity. Abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1; N, number of chromosomes sampled, θS Watterson’s estimator of θ, θ, Nei’s estimator of θ.

Table 3

Rank order comparisons of genetic diversity in four lizard species.

	Haplotype diversity				Nucleotide Diversity ΘS
	cytb	CARD4	MC1R	RAG1	cytb	CARD4	MC1R	RAG1
Species
C. b.	2a	4a	1ab	4a	1a	2	1a	2
G. w.	1bc	1abc	2cd	1abc	2b	1	2b	1ab
D. d.	4ab	2b	4ace	2b	4abc	4	4ab	3a
S. a.	3c	3c	3bde	3c	3c	3	3	4b
Pop. size
Small	1.5	2.5	1.5	2.5	1.5	1.5	1.5	1.5
Large	3.5	2.5	3.5	2.5	3.5	3.5	3.5	3.5
Habitat
Rock	2.5	1.5	3.0	1.5	3.0	2.5	3.0	2.0
General	2.5	3.5	2.0	3.5	2.0	2.5	2.0	3.0

Values in table are species’ diversity ranks within a locus, with 1 being the most diverse. Species are listed with the smaller population size species (C. b. and G. w.) first. Within each column, ranks that share a letter are statistically different (α ≤. 0.05). None of the ΘN estimates were significantly different and are not shown. Lower half of table shows the average ranks of the two small population size species and two larger population size species, respectively. Abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1; N, number of chromosomes sampled, θS, Watterson’s estimator of θ; C. b., Crotaphytus bicinctores; G. w., Gambelia wislizenii; D. d., Dipsosaurus dorsalis; S. a., Sauromalus ater.

Mitochondrial DNA

The rank order of haplotype diversities was G. wislizenii > C. bicinctores > S. ater > D. dorsalis,(Table 2, Fig 2) nearly the reverse of the expectation of Neutral Theory based on population size rank. Haplotype diversity of G. wislizenii was significantly greater than that of all other species, while C. bicinctores diversity was significantly greater than that of S. dorsalis but not S. ater.

Fig 2

Haplotype diversity and Watterson’s estimate of θ for four species of lizard from the Mojave National Preserve.

Solid black line separates the low population size carnivores (left) from the higher population size herbivores (right). Letters within bars (shown only on the haplotype diversity graph for readability) signify habitat specialist (S) or generalist (G). Gene abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1Species abbreviations: CRBI, Crotaphytus bicinctores; GAWI, Gambelia wislizenii; DIDO, Dipsosaurus dorsalis; SAAT, Sauromalus ater. Error bars represent 95% confidence intervals.

Watterson’s and Nei’s θ showed patterns similar to that of haplotype diversity for cytb (Tables 2 and 3). C. bicinctores and G. wislizenii had the highest θ values but were statistically indistinguishable. Relative values of Nei’s θ were similar. However, due to the higher standard error associated with this estimator, no significant pairwise differences were found.

Tajima’s D was not significantly different from 0 for any species (α = 0.05), suggesting that a stable population size neutral model is appropriate for this locus (Table 4). Conversely, raggedness index values were significant at an alpha level of 0.05 only for S. ater, suggesting that cytb in the other species have expanding populations. The potential cause for the discrepancy between Tajima’s D and the pairwise mismatch results is discussed in the following section.

Table 4

Tests of neutral model with constant population size.

Species	Gene	D	P-value	R	P-value
Crotaphytus bicinctores	cytb	-0.857	0.177	0.120	1.000
	CARD4	-1.087	0.143	0.120	0.400
	MC1R	-0.460	0.365	0.103	0.000
	RAG1	-1.399	0.074	0.094	0.094
Gambelia wislizenii	cytb	-0.092	0.534	0.040	0.318
	CARD4	-0.635	0.303	0.115	0.169
	MC1R	-0.638	0.288	0.038	0.200
	RAG1	-1.480	0.056	0.085	0.012
Dipsosaurus dorsalis	cytb	1.305	0.081	0.213	0.091
	CARD4	0.692	0.233	0.071	0.763
	MC1R	-0.124	0.478	0.130	0.163
	RAG1	1.411	0.075	0.117	0.046
Sauromalus ater	cytb	-0.134	0.515	0.526	0.005
	CARD4	-0.551	0.340	0.121	0.999
	MC1R	-0.976	0.186	0.170	0.002
	RAG1	0.077	0.415	0.085	0.412

Abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1; D, Tajima's D; R, Harpending's raggedness index. Significant p-values (α = 0.05) in bold.

The ratio of synonymous to nonsynonymous changes was high in every species except D. dorsalis, in which the only variable site had a nonsynonymous mutation (Table 1). Across all four species, 21.1% of base changes in cytb were nonsynonymous.

Autosomal DNA

The length of the gene segment sequenced ranged from 787–807 base pairs in CARD4, 642–813 bp in MC1R, and 774–1049 bp in RAG1. All genetic diversity measurements for autosomal loci were calculated from the most likely haplotype pairs identified by PHASE. As expected, PHASE was generally able to find best haplotypes pairs with higher levels of confidence for locus-species combinations with lower genetic diversity (especially lower heterozygosity). The percentage of individuals with a best haplotype pair supported with 90% or greater confidence ranged from 100% (C. bicinctores CARD4, D. dorsalis CARD4 and MC1R) to 52% (C. bicinctores MC1R and G. wislizenii MC1R). Results from haplotype randomizations showed that using the best pairs from PHASE was unlikely to yield incorrect diversity estimates even if the best pair has low probability (S2 Table). The ranges of haplotype diversity values in the randomized datasets were quite narrow and were within the 95% confidence interval of the best pairs value for each locus in all cases except MC1R in C. bicinctores and G. wislizenii (and even in these cases there was substantial overlap), showing that the choice of haplotype pairs does not have a major effect on the diversity results. Nei’s θ values were either the same for every randomized dataset (in eight of eleven locus-species combinations), or formed a narrow range. Watterson’s θ values did not vary in any randomized dataset because it is invariant to linkage phase. Therefore, Tajima’s D values varied little among randomized datasets.

The average haplotype diversity averaged across species was highest at the locus RAG1 (0.779) and lowest at CARD4 (0.568). Significant pairwise differences among some species were found at each autosomal locus (Table 3). In both CARD4 and RAG1, the haplotype diversity of G. wislizenii was significantly higher than that of the other three species, which were statistically indistinguishable from one another. At MC1R, there was no difference in the very high h values in crotaphytids, but both iguanids were significantly different from one another and significantly lower than the crotaphytids.

For autosomal loci, Watterson’s θ values showed significant differences in the many of the pairwise comparisons between species. At MC1R, both crotaphytids were significantly more diverse than D. dorsalis but not S. ater. At RAG1, G. wislizenii had the highest value and was statistically higher than either iguanid. At CARD4, there were no significant pairwise differences between species. As was the case for cytb, Nei’s θ values showed no significant differences between any species pairs.

Values of Tajima’s D were not significantly different from 0 (α = 0.05) for any locus-species combination, showing no evidence of selection or changes in population size. The highest absolute value of Tajima’s D was found in G. wislizenii RAG1 (-1.480; p = 0.056).

Harpending’s raggedness index results were inconsistent. All species examined had a significant R value at one or more loci (α = 0.05). This test assumes a null hypothesis of an expanding population, and a significant deviation from the null indicates a stable population size. Thus, as a test of population expansion, R is prone to type II errors, and the lack of consistency among the loci of each species may be the result of a lack of statistical power to detect a multipeaked distribution for some loci. Since at least one locus for each species had a significant R value, and no species had significant Tajima’s D values, it seems reasonable to conclude that population sizes for these species have been essentially stable in recent time.

The ratio of synonymous to nonsynonymous base changes varied greatly among loci and species. For most species, the average across loci was consistently around 30% nonsynonymous (29.0% - 31.9%), but D. dorsalis had a much different ratio, with 61.1% of changes being nonsynonymous. When calculated by locus, it was found that the average proportion of nonsynonymous base changes across all species were 36.7% in CARD4, 23.3% in MC1R, and 58.0% in RAG1.

Population subdivision

In order to test whether the high levels of genetic diversity seen in these species (particularly the crotaphytids) were simply the result of combining two genetically distinct subpopulations, we recalculated basic population genetic measurements for each putative subpopulation (defined by collection locality). Overall, haplotype diversities of the potential subpopulations were similar to those of the pooled samples for each species. Haplotype networks visually show the proportion of haplotypes which came from each collection locality (S1 Fig), and suggest very little population subdivision of haplotypes. C. bicinctores had the most differentiation, with significant differences in haplotype diversity between the two sampling locations at two loci (cytb and CARD4). Even at these two loci, however, one of the two subpopulations had a haplotype diversity value statistically indistinguishable from that of the total sample, showing that the high total diversity is not simply an additive effect. In addition, C. bicinctores had the lowest diversity rank at CARD4 despite the differences in haplotype diversity from the two sampling sites. FST values ranged from 0.011 (RAG1) to 0.249 (cytb), and ΦST ranged from 0.083 (RAG1) to 0.177 (cytb) (Table 5). ΦST values were significantly different from zero at all loci in C. bicinctores (α = 0.05). Jost’s D ranged from 0.044 (RAG1) to 1.000 (in cytb, where no haplotypes were shared between the two subpopulations). Effective number of subpopulations, ΔST, ranged from 1.039 for RAG1 to 2.000 for cytb.

Table 5

Measures of population subdivision.

Species	Gene	FST	ΦST	Jost's D	ΔST
Crotaphytus bicinctores	cytb	0.249	0.177	1.000	2.000
	CARD4	0.061	0.107	0.136	1.082
	MC1R	0.037	0.172	0.692	1.580
	RAG1	0.011	0.083	0.044	1.039
Gambelia wislizenii	cytb	0.006	0.023	0.069	1.061
	CARD4	0.008	0.000	0.034	1.048
	MC1R	0.009	0.003	0.182	1.225
	RAG1	0.001	0.000	0.024	1.133
Dipsosaurus dorsalis	cytb	0.000	0.000	0.000	1.006
	CARD4	0.023	0.023	0.054	1.047
	MC1R	0.000	0.000	0.000	1.008
	RAG1	0.067	0.000	0.292	1.204

Abbreviations: cytb cytochrome b, CARD4 caspase recruitment domain gene 4, MC1R melanocortin 1 receptor gene, RAG1 recombination activating gene 1.

In D. dorsalis, there were no differences among subpopulation and total haplotype diversity values at any locus. FST values ranged from 0 (cytb and MC1R) to 0.067 (RAG1), and all ΦST values were statistically indistinguishable from zero. Jost’s D ranged from 0 (cytb and MC1R) to 0.292 (RAG1). ΔST ranged from 1.006 for cytb to 1.204 for RAG1.

G. wislizenii samples also lacked any significant differences in haplotype diversity among subpopulations and the total. FST values ranged from 0.001 (RAG1) to 0.009 (MC1R), and no ΦST values were significantly different from zero. Jost’s D ranged from 0.024 (RAG1) to 0.182 (MC1R). ΔST ranged from 1.048 for CARD4 to 1.225 for MC1R.

Effect of trophic level/population size

Genetic diversity patterns between species differing in population size were nearly the opposite of what would be expected under neutral theory. Point estimates of the two small N species were ranked first and second most diverse for two of 4 genes (haplotype diversity) and for all 4 genes for ΘS. The average rank of the two small N species was greater than or equal to that of the large N species in all cases (Table 3). An alternative way of grouping the focal species is by habitat preference. Within a desert ecosystem, D. dorsalis and G. wislizenii are habitat generalists, and are often found in the open, while C. bicinctores and S. ater are saxicolous, found primarily among large rocks. Although rocky habitat is common at the study site, in general, it is a small fraction of the total habitat and we would expect habitat generalists to have larger overall populations than habitat specialists. The results do not show a consistent relationship between habitat specialization and genetic diversity (Table 3).

Comparison with other species

We were able to compare haplotype diversity found in the present study with that found in other, more abundant lizard species from the same locality, sampled in previous studies (Fig 3; CARD4 was not sequenced in previous studies [19,20]). To give an indication of the relative sizes of these species, we present the capture numbers from the same can trap study [35] cited above (Callisaurus draconoides—92, Coleonyx variegatus 108 and Uta stansburiana—1,529). The expectation was that the four larger bodied species studied here would have lower genetic diversity than all of the smaller, more common lizards studied. This was true for cytb, the only mitochondrial locus sequenced, where nearly all of three smaller species had significantly greater haplotype diversity than the four larger species examined here (the one exception was C. variegatus, with a haplotype diversity statistically indistinguishable from G. wislizenii).

Fig 3

Haplotype diversities of 4 focal species compared to more abundant lizard species from Mojave National Preserve.

Error bars represent ± 1 standard error. Abbreviations: cytb, cytochrome b; CARD4, caspase recruitment domain gene 4; MC1R, melanocortin 1 receptor gene; RAG1, recombination activating gene 1.

This was not the case for the autosomal loci, however. The haplotype diversity of the two crotaphytids at MC1R was statistically indistinguishable from the values found in U. stansburiana and C. variegatus. At RAG1, levels of haplotype diversity in C. bicinctores, D. dorsalis, and S. ater are all indistinguishable from C. draconoides as well as one another (RAG1 was not sequenced in U. stansburiana). Interestingly, haplotype diversity at this locus in G. wislizenii was significantly greater than that of C. draconoides, one of the most abundant lizards at this locality.

Discussion

Our results do not conform to neutral expectations regarding the relationship between genetic diversity and population size. Among the four focal species examined for this paper, the species with the smallest population sizes (based on observation of current densities and the proxy criterion of trophic level) have diversity levels higher than those of species with larger population sizes. The conclusion from comparing the 4 uncommon species surveyed for this study to the previously published results [19,20] is the same as the “matched comparison” within the 4 focal species of this paper. There is no clear relationship between the rank order of local population size. Indeed, these diversity levels are close to those of the lizard species with the very highest population sizes (Fig 3). A different population size proxy for the four study species would have been to categorize them as either habitat generalists (Dipsosaurus dorsalis and Gambelia wislizenii) or saxicolous (specialized for living around rocks; Crotaphytus bicinctores and Sauromalus ater). The generalist species would be predicted to have a larger population size, since more usable habitat is available to them, and therefore they should have higher genetic diversity under the neutral theory. However, the data show no evidence for this pattern either. The contrast in Table 3 between habitat generalists and specialists shows no consistent relationship across genes.

Although our results are more consistent with the predictions of selective models, there was no strong evidence for selection at the loci sequenced. This does not preclude the possibility that the diversity at these loci was shaped by selection at other, linked loci. The one exception was the gene RAG1, which may be under selection in these lizards. The highest proportion of nonsynonymous base changes was found at this locus, suggesting that there may be functional differences among the copies sequenced. Similar patterns were seen in C. draconoides and C. variegatus, where the percentage of nonsynonymous mutations ranged from 3.3–26.4% in cytb and MC1R, but were at 57.9–58.8% in RAG1 [20].

We attempted to select markers that had as little linkage among loci as possible. Our autosomal loci BLAST to different chromosomes on the Anolis carolinensis genome (RAG1- Chromosome 1, CARD4—Chromosome 6, MC1R –unplaced genomic scaffold). If genetic draft were responsible for the patterns seen here, selective sweeps would likely have been pervasive, affecting much of the genome, rather than a few isolated events in the lineages studied. Interestingly, the variability at cytb, the mitochondrial locus, was not less than that of the autosomal loci. Unless there are great differences in mutation rate, it would be expected that genetic draft would have a greater effect on diversity of mtDNA because of the strong linkage disequilibrium among mitochondrial loci [7]. The mutation rate of mtDNA is expected to vary greatly among taxa [2] however, and we only sequenced one mitochondrial locus. Associative overdominance has been suggested as another mechanism by which genetic diversity in small populations may be greater than predicted by neutral theory [49].

Selection is not the only possible explanation for the lack of correlation between population size and genetic diversity in these lizard populations. Despite the fact that all the populations were sampled from the same small geographic area, it is possible that the different species experienced different phylogeographic or demographic events, such as vicariance or bottlenecks, or that they have very different mutation rates. However, Tajima’s D showed no evidence of a bottleneck in any species. Pairwise mismatch analysis was unable to detect the multipeaked mismatch distribution characteristic of a stable population, but given our sample sizes this is likely to be lack of statistical power (Type II error) rather than positive evidence of a recent population expansion.

Our research was relatively limited in the number of species and loci surveyed, and the noise from coalescent variance(the among-locus variance in time to most recent common ancestor of a sample of DNA sequences from a locus) could mask the relationship between population size and genetic variation. An interesting next step for research with this study system would be to use next generation sequencing to sequence full genomes of all individuals. Full genome data would greatly increase the amount of data available to draw conclusions from, and provide several additional advantages. It would be possible to get a clearer view of how mitochondrial DNA (and other non-recombining segments like sex chromosomes) and nuclear genomes differed in genetic diversity, and intraspecific variation in coding genes could be compared with that found in putative noncoding regions. It would also be possible to get good estimates of the effects of coalescent variance on the variation in genetic diversity among loci [50]. A recently published paper by Grundler et al. [51] took advantage of a long term census study of an Australian lizard community to examine the relationship between population size and genetic diversity using next generation sequencing data. They found a weak relationship between abundance and genetic diversity. They postulate that “nucleotide diversity is heavily influenced by factors other than census population size, or that ecological sampling in this community is unable to capture true population size.” Interestingly, Grundler et al. found that genetic diversity was more strongly correlated with occupancy than with abundance. Assuming that occupancy is correlated with population connectivity, gene flow differences may be having an effect on genetic diversity by raising the effective population size for local population with higher gene flow rates. Our data are somewhat consistent with this idea, because G. wislizenii has the largest range of the four focal species and the lowest Fst values in our two location comparisons. As a habitat generalist, G. wislizenii may have a greater propensity for gene flow. Migration from outside populations may compensate for lower local populations sizes for some species. This may be why G. wislizenii has genetic diversities approaching those of the three species with much greater local population size (Fig 3), although gene flow is unlikely to explain the relatively high genetic diversity of C. bicinctores. As a rock specialist, we expect C bicinctores to have lower gene flow from neighboring populations because of intervening unsuitable habitat, and this is consistent with the higher Fst values for this species (Table 5).

Although our study does not by itself provide a definitive answer about the processes shaping genetic diversity, it is of interest because its results deviate from neutral expectations. It does not appear that population size has had a large effect on genetic diversity in this system, and that suggests that selection, demographic, and/or mutational differences may be playing an important role. If additional studies show a similar lack of relationship between population size and genetic diversity for matched sets of species, caution should be used when implementing algorithms to estimate population parameters like effective population size, time to most recent common ancestor, bottleneck effects, etc., as these methods assume relatively simple neutral models.

Specimen collection localities.

Collection locality for each lizard sampled, by species.

(DOCX)

Click here for additional data file.

Haplotype randomization results.

Genetic diversity measures are shown for best supported pairs, along with the range of measures among simulated datasets.

(XLSX)

Click here for additional data file.

Haplotype networks.

Haplotype networks for all species, showing the haplotypes by collection locality. Within each network, circles are proportionate to number of copies of a haplotype (networks not on same scale). Hash marks between haplotypes represent the number of mutational steps.

(PDF)

Click here for additional data file.

Supplementary methods–Primers and PCR protocols.

(DOCX)

Click here for additional data file.

R script for haplotype randomizations.

(DOCX)

Click here for additional data file.

R script for calculation of FST and Jost’s D.

(DOCX)

Click here for additional data file.

27 Jun 2019

PONE-D-19-14658

Genetic variation across trophic levels: a test of the correlation between population size and genetic diversity in sympatric desert lizards

PLOS ONE

Dear Dr. Routman,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

==============================

I have now received reviews from two reviewers. The comments are very positive from the reviewers and I agree with them. However, there are several valid concerns raised by the reviewers, and I would suggest you to address them thoroughly before we can accept the manuscript. Therefore my decision is a major revision.

I agree with the reviewers that the manuscript needs a  (i) map with details of sampling locations and habitat types, (ii) explicitly stated hypothesis and predictions of this manuscript, (iii) changing some tables to figures for better presentation of the results, and (iv) better presentation of abundance and habitat specialization proxies. In addition, I would suggest to clearly and explicitly state the expectations as mentioned by one of the reviewers. For example, make a statement at the end of the introduction that “We tested neutral vs selection theories for relationship between genetic diversity and (a) body size, (b) population size, (c) trophic level, and (d) habitat specialization. Our expectation is that ….. under population size, genetic diversity would be …Carnivore specialist< carnivore generalist< herbivore specialist< herbivore generalist. “

I suggest to refrain from reference to ‘our lab’ and replace it with our previous research (citation) or previous research (citation). Also, please be consistent in how you use and present the species names: and Genus species, or G species, or G.species, or Gs: choose one and be consistent with it please.

Congratulations on the positive reviews, I hope you will be to address the issues highlighted by the reviewers. Details on submitting the revisions are enclosed.

============================

We would appreciate receiving your revised manuscript by Aug 11 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Trishna Dutta

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We noted in your submission details that a portion of your manuscript may have been presented or published elsewhere. Please clarify whether this publication was peer-reviewed and formally published. If this work was previously peer-reviewed and published, in the cover letter please provide the reason that this work does not constitute dual publication and should be included in the current manuscript.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: This is an interesting study that explores the relationship between population size and genetic diversity with four sympatric species- two herbivorous and two carnivorous lizards. I have a few broad and some specific comments.

Introduction

Generally, the introduction seems to be guided by the results that the authors got. If the motivation behind the study was to test whether neutral or selective forces shape genetic variation, this system and study design may not be most suited to address that.

In the introduction, the authors could focus more on the factors affecting genetic diversity in natural populations. It would be useful to have more information about both the intrinsic and extrinsic factors that affect genetic diversity. Authors mention the intrinsic factors of body size, mutation rates, generation times, etc. But extrinsic factors, such as extent and continuity of available habitat, do not find any mention.

66-67 Please elaborate on/ give examples of the confounding variables that you mention

Methods

91-92 It would be good to mention the body sizes here. People unfamiliar with the species and the region might find it useful.

Results

It would be nice to have a map showing study area with sample locations of the different species in different colours

It would be nice to compare generalist and specialist species within a trophic level instead of pooling them across the two levels as done in table 3.

A lot of information in the first paragraph under the ‘Mitochondrial DNA’ subheading can be moved to methods.

283-288 This information should go under the collection section in methods. Please also mention the distance between the different collection sites. Why was only one species sampled from Zzyzx? Other species do not occur in that area or were they not sampled for a reason?

Discussion

361-363 ‘The generalist species would be predicted to have a larger population size, since more usable habitat is available to them, and therefore they should have higher genetic diversity under the neutral theory. However, the data show no evidence for this pattern either.’ The authors cannot ignore the trophic level while making this statement. It is not okay to club the carnivore and herbivore generalists since population size is impacted by a combination of both- specialization and trophic level. Does the statement imply that the generalist herbivore in their study had a larger population size than the specialist carnivore? If yes, then its okay to club them in the same category and comment on expectations based on neutral theory. Whereas if they expect the population size in the following order: Carnivore specialist< carnivore generalist< herbivore specialist< herbivore generalist, the above statement cannot be made. The authors say nothing about the order in which they expect the population sizes to be for their study species. It should be explicitly mentioned in the introduction.

G. wislizenii has the lowest genetic differentiation and the highest genetic diversity based on multiple estimators. Does it also, among the study species, have a more continuous distribution and larger range?

Grundler et al (2019) does not find a reference in your paper. This recent study looks at the relationship between genomic diversity and abundance, occupancy and habitat specialization using genomic data from 30 species of lizards from arid Australia.

Reviewer #2: This study presents valuable data on genetic diversity, with an interesting ecological component that seeks to shed light on the still equivocal relationship between heterozygosity and population size. As such, it represents an important contribution to the field. However, there are several major points that I feel need reworking. My recommendation is publication after major revision, which should address the following issues:

1. The bulk of the data are presented in tables, which in my opinion lessens both the digestibility and impact of the results. The single figure that is presented could be improved with a key that indicates population size category of each of the species, such that readers can immediately pair heterozygosity with abundance/trophic category.

2. There are issues with the presentation of abundance data, which I feel impedes reproducibility of the study. The abundance data on which the authors base population size classification include a census study that I was unable to access; observational anecdotes (albeit over a long period of time); and coarse proxies. Proxies can be both necessary and valuable, but should be supported by more data or literature reference than found here. Similarly, observations of abundance should be supported with stronger justification as is done in Hague and Routman (2016) (frequently referenced in the current manuscript).

3. Because of the above two issues, I found it necessary to read Hague and Routman (2016) in conjunction with the manuscript in question. Much of the data and ecological context for the current study is presented in the 2016 manuscript. The two studies together present an interesting picture of the ambiguity of the relationship between genetic diversity and population size, and highlight the need for further study of ecological traits and genetic diversity. The current study seems to be a logical continuation of the 2016 project, and Figure 1 is a good example of the integration of data from each of these projects. However, the current manuscript could be improved by elaborating on the combined results in the discussion.

Specific comments in attached document.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: PONE-D-19-14658_reviewer.docx

Click here for additional data file.

Submitted filename: PLOS_review_trophic_gendiv.docx

Click here for additional data file.

2 Sep 2019

28 August 2019

Dear Dr. Dutta:

We would like to thank you and the reviewers for your careful review of our manuscript. We have revised the manuscript incorporating most of the comments from the reviews. This cover letter details the responses, organized by reviewer. The reviewers’ comments are in quotes and smaller font, followed by our response.

Reviewer 1

“Generally, the introduction seems to be guided by the results that the authors got. If the motivation behind the study was to test whether neutral or selective forces shape genetic variation, this system and study design may not be most suited to address that.”

This is the original motivation for this study, which is part of a larger study on multiple species in the same area. This study adds to work we have previously published (Micheletti et al. 2012, Hague and Routman 2016) that examines genetic diversity in sympatric populations of lizards. We are a Master’s only institution, and to accommodate the shorter time period that those students matriculate, individual projects need to be smaller than one would expect from a Ph.D. institution. We feel that the total system is actually quite appropriate for this question, but of course the comparison of genetic diversity to population size is only part of testing neutrality vs. selection, and we anticipate using next generation sequencing to test for the fingerprint of selection in future work. It should be noted that the main stimulus for theorists to develop selective theories for DNA sequence diversity is the supposed lack of correlation, predicted by neutral theory, between population size and genetic diversity when comparing taxa from very different parts of the tree of life. We are simply arguing that comparison of populations of closely related species from identical areas would be a better test of neutral theory. Such comparisons should be done by many different labs for taxonomically divergent groups before a definitive conclusion can be drawn about genetic diversity vs. population size.

“Authors mention the intrinsic factors of body size, mutation rates, generation times, etc. But extrinsic factors, such as extent and continuity of available habitat, do not find any mention.”

“66-67 Please elaborate on/ give examples of the confounding variables that you mention”

The paragraph the reviewer is referring to is from Line 65-70. This is our attempt to justify 1) using taxonomically related species, and 2) using sympatric species, based on the criticisms from the previous paragraphs. We suggest that “focusing on a taxonomically related group” eliminates confounding variables that the reviewer correctly identifies as “intrinsic”. But we also suggest that studying species from a single area reduces some “extrinsic” forces like different geological histories or levels of gene flow (lines 69-71). We have altered the language to make this clearer.

“91-92 It would be good to mention the body sizes here. People unfamiliar with the species and the region might find it useful. “

We have added body sizes here.

“It would be nice to have a map showing study area with sample locations of the different species in different colours.”

We have added a map to show the desert region (satellite view) and the collection localities of each species.

“It would be nice to compare generalist and specialist species within a trophic level instead of pooling them across the two levels as done in table 3.”

This comparison can be made from Table 3, but we have made it more explicit by showing the results for each species/gene combination as a Figure. (We agree with all the reviewers that showing the results as only Table 3 is not the best way to convey the comparisons.). The pooled comparison is to show that even if we are wrong about the relative population sizes and both of the habitat specialists have lower population sizes, population size still does not explain the differences in genetic diversity.

“A lot of information in the first paragraph under the ‘Mitochondrial DNA’ subheading can be moved to methods”

We have moved some of this info to Methods, leaving only aspects which were done because of the specific results that were obtained.

“283-288 This information should go under the collection section in methods. Please also mention the distance between the different collection sites. Why was only one species sampled from Zzyzx? Other species do not occur in that area or were they not sampled for a reason?”

We only collected Crotaphytus from Zzyzx because the species is uncommon enough that we had to collect from a different locality to get a sample size that approached that of the two herbivorous species. Gambelia was also less common than the herbivorous species, but is not common at Zzyzx. That is why we collected Gambelia at Kelso Dunes to help fill out this sample. When the study started, we were expecting that Gambelia and Crotaphytus would be less variable than the more common herbivorous species, even when combining samples from different areas, and therefore the combined samples would be a conservative test. This turned out not to be true, so we needed to do the analysis with and without combined samples, as described in the manuscript in the Results, subsection Population subdivision.

“361-363 ‘The generalist species would be predicted to have a larger population size, since more usable habitat is available to them, and therefore they should have higher genetic diversity under the neutral theory. However, the data show no evidence for this pattern either.’ The authors cannot ignore the trophic level while making this statement. It is not okay to club the carnivore and herbivore generalists since population size is impacted by a combination of both- specialization and trophic level. Does the statement imply that the generalist herbivore in their study had a larger population size than the specialist carnivore? If yes, then its okay to club them in the same category and comment on expectations based on neutral theory. Whereas if they expect the population size in the following order: Carnivore specialist< carnivore generalist< herbivore specialist< herbivore generalist, the above statement cannot be made. The authors say nothing about the order in which they expect the population sizes to be for their study species. It should be explicitly mentioned in the introduction.”

We have rewritten this to be clearer that we were referring to a possible alternative contrast based on habitat specialization only, not the original hypothesis that carnivores are less common than herbivores and that within this contrast, habitat specialists are rarer than generalists. That is, we did not get the result that the observed-to-be-rarer species, the carnivorous lizards, were less diverse than the more common, herbivorous species. So one way to get our unexpected result would be if the habitat specialists as a group were actually less common than the habitat generalists. But that contrast did not match neutral theory explanations either. In fact, our genetic diversity results are not consistent with any contrast based on proxies for population size and neutral theory.

“G. wislizenii has the lowest genetic differentiation and the highest genetic diversity based on multiple estimators. Does it also, among the study species, have a more continuous distribution and larger range?”

We have included additional discussion of the possible role of gene flow in the Discussion. G. wislizenii does have a larger range than the three other focal species, although this may be an artifact of lower taxonomic/biogeographic attention that this species has received. For example, a mtDNA study (COIII) by Orange et al showed that the Mojave/Great Basin/Colorado plateau populations are separated from the Chihuahuan Desert populations by 6.1% sequence divergence. This suggests that gene flow from the Chihuahuan populations is unlikely to be influencing diversity in the more western populations. If the Chihuahuan Desert populations were designated as a separate species, there would not be as much difference in range sizes. The difficulty is that funding thorough multi-locus biogeographic studies is not a priority, so these studies may never be done.

“Grundler et al (2019) does not find a reference in your paper. This recent study looks at the relationship between genomic diversity and abundance, occupancy and habitat specialization using genomic data from 30 species of lizards from arid Australia. “

I only became aware of the Grundler et al. paper after we submitted the manuscript. We have incorporated some discussion of their paper.

Reviewer 2

“The bulk of the data are presented in tables, which in my opinion lessens both the digestibility and impact of the results. The single figure that is presented could be improved with a key that indicates population size category of each of the species, such that readers can immediately pair heterozygosity with abundance/trophic category”

We have added figures to correct this, as stated under “Reviewer 1”

“There are issues with the presentation of abundance data, which I feel impedes reproducibility of the study. The abundance data on which the authors base population size classification include a census study that I was unable to access; observational anecdotes (albeit over a long period of time); and coarse proxies. Proxies can be both necessary and valuable, but should be supported by more data or literature reference than found here. Similarly, observations of abundance should be supported with stronger justification as is done in Hague and Routman (2016) (frequently referenced in the current manuscript).”

The abundance data are not great. The published survey (Persons and Nowak, 2007) is mostly consistent with my 20+ years of observation at this locality, but they did not show some very abundant taxa, such as Coleonyx variegatus, as being among the most common species, suggesting a collecting method bias or year to year variation in apparent abundance. We found a publication that includes the Wallace can trap data from the Zzyzx area, and we have added that to our revision. However, can traps are not effective at collecting Sauromalus, so the data include only 2 specimens of this species. We have added discussion of the fact that it took 4 years to collect the carnivore samples and much of that collecting was focused specifically on the carnivore species, whereas the herbivore sampling took less than two years and was done early in the study, when we were collecting all species of lizards and therefore not always in the best microhabitats for the two herbivore species that are the focus of this manuscript.

“Because of the above two issues, I found it necessary to read Hague and Routman (2016) in conjunction with the manuscript in question. Much of the data and ecological context for the current study is presented in the 2016 manuscript. The two studies together present an interesting picture of the ambiguity of the relationship between genetic diversity and population size, and highlight the need for further study of ecological traits and genetic diversity. The current study seems to be a logical continuation of the 2016 project, and Figure 1 is a good example of the integration of data from each of these projects. However, the current manuscript could be improved by elaborating on the combined results in the discussion.”

We have elaborated on the combined results in the discussion.

“L41-43 – what are the hypotheses for deviations from expectation in these studies? Good to have context here about why we might not expect things to correlate”

We added a statement that the authors credit natural selection for the deviations from neutral expectations. However, the possibility of coalescent variation (variation among loci in time to most recent common ancestor, which can vary enormously even for different loci sampled from exactly the same population) and variation in mutation rates is also present.

L59 – not sure about the use of “haphazard” here; larger scale analyses are important for other reasons, especially considering impacts of population connectivity on genetic diversity

We eliminated the term haphazard.

L117-118 – references for this?

We have added two references for this statement.

L118-120 – References for census studies are inaccessible. I feel more detail is needed here to justify the validity of qualitative data – observations on abundance across seasons? Over periods of how long? Etc. Even a figure plotting the cited census data with genetic diversity would be valuable.

We have added more information about our justification for relative population sizes, including another reference to Wallace’s can trap study. [Cummings KL, Puffer SR, Holmen JB, Wallace JK, Lovich JE, Meyer-Wilkins K. Biodiversity of amphibians and reptiles at the Camp Cady Wildlife Area, Mojave Desert, California and comparisons with other desert locations. California Fish and Game 2018 Summer 104(3): 129-147].

L124 – how large was the main collection area, and what are overall geographic ranges like for these species? More detail on connectivity of these species with surrounding populations would be useful when considering genetic diversity and population size (especially when there are deviations from expectation as we have here)

We have no data on the connectivity among surrounding populations (although we do Fst analysis of the two species with samples from multiple locations. As noted above, we have added a discussion of the possible role of gene flow to the Discussion.

L147-174 – good use of alternate tests and simulations to validate analyses

Thanks.

L160-162 – there are standard statistical tests to quantify significant difference in pairwise measures that are sounder than a decision made by visual inspection. The citation for this method comes from a psychology article that I am unable to access, and seems, based on the abstract, more relevant for gaining inference from figures and improving research communication in a field that operates under different methods standards.

The standard error approach is not simply visual, and involves a calculation of the overlap. However, it does make assumptions about the shape of the credible interval, especially that the variation is normally distributed. We have calculated the t test suggested by Masatoshi Nei for the haplotype diversity contrasts, and found some additional significant differences. However, because the direction of the difference in diversity is counter to the predictions of neutral theory, our original method is conservative. That is, discovering additional statistical significance does not change the conclusion that the diversities do not match those predicted by neutral theory.

L187:

Table 1 - Dn/ds statistic could be presented as a single column rather than or in addition to #synonymous and #nonsynonymous mutations in separate columns (which makes interpretation of this metric a two-step rather than one-step process)

Done

Tables summarized into figures would be helpful, especially Table 3. Much of this appears in fig 1, but the summary of information could be improved; cannot immediately tell which are small pops, which are large pops, which trophic level the species represent, etc. It’s not a lot to keep track of or to find in the text, but figures should be interpretable by themselves.

Done (see response to Reviewer 1 for the same criticism).

“L236 - Tables 4 and 5 are identical”

Oops! The correct Table 5 has been added.

“L334 – references for eco data? Found in open habitat doesn't necessarily = generalist…if rocky habitat is common then rock specialists could also be common…”small fraction of total habitat” – any kind of quantification available? What is the rest of the habitat characterized by? Ecological specialization is a very interesting and important consideration, and a valuable addition to the study, but it does feel a little too coarse as it stands. It would be helpful to at least include some references to strengthen inferences made here. “

We have added language to our discussion of the collection maps to discuss the fact that the area of suitable habitat for the rock dwelling species is less than that of the species that prefer the open areas. In addition, the habitat generalists also make use of the rock areas in addition to the open areas. Habitat specialization is simply more justification for our ranking of relative population sizes. (I am aware of no population genetics theory that relates the degree of habitat specialization to genetic diversity, other than through its effect on population size).

L386-387 - Reasoning behind “likely lack of statistical power”? Authors have a good number of individuals - do they mean due to number of species? Could use a little clarification

This is not a “number of species” problem. We discuss this a bit in the results section (Lines 270-276 of the original manuscript). This test relies on the ability to differentiate a multi-modal distribution (of pairwise nucleotide mismatches among sampled copies of a locus) from a unimodal distribution. Differentiating distributions require a large number of sampled copies to have strong statistical power. In this case we are trying to infer a population expansion by detecting a unimodal distribution. In the statistical test, the unimodal distribution is the null hypothesis, so that failure to reject the null could be due to population expansion or lack of power.

L395 – need some background here to contextualize coalescent variance

Coalescent variance is the term to describe the variation in TMRCA (Time to Most Recent Common Ancestor) among loci sampled from a common ancestor. As a result, it follows that there can be differences in genetic diversity of a given locus sampled from different populations/species just because the locus has a recent TMRCA in one population and a more ancient TMRCA in another population, simply by chance. I believe that the term “coalescent variance” is in common use in population genetics, but we can include a paragraph like the one in this response if the editor feels that it would be useful. For this revision we have defined it in a few words in the text on its first mention (which is L389 in the original manuscript) as “coalescence variation”, which I have changed to “coalescent variance” for consistency.

L407 – It would be valuable to situate the results of this study better in the context of previous results – Hague and Routman (2016) showed clear relationship between diversity and pop size but now we have conflicting relationship when we throw smaller pops into the mix – what could be going on here? This is especially interesting and merits additional study with all species to look for ecological correlations – and so should be mentioned with more detail in the discussion

We have elaborated on the discussion to discuss the other species for which we have data and the paper mentioned by reviewer 1.

Associate Editor Dutta

“I agree with the reviewers that the manuscript needs a (i) map with details of sampling locations and habitat types, (ii) explicitly stated hypothesis and predictions of this manuscript, (iii) changing some tables to figures for better presentation of the results, and (iv) better presentation of abundance and habitat specialization proxies. In addition, I would suggest to clearly and explicitly state the expectations as mentioned by one of the reviewers. For example, make a statement at the end of the introduction that “We tested neutral vs selection theories for relationship between genetic diversity and (a) body size, (b) population size, (c) trophic level, and (d) habitat specialization. Our expectation is that ….. under population size, genetic diversity would be …Carnivore specialist< carnivore generalist< herbivore specialist< herbivore generalist”

Items i, iii, and iv are address in responses to Reviewers 1 and 2. We have added language to address explicit hypotheses and expectations.

“I suggest to refrain from reference to ‘our lab’ and replace it with our previous research (citation) or previous research (citation).”

Done.

“Also, please be consistent in how you use and present the species names: and Genus species, or G species, or G.species, or Gs: choose one and be consistent with it please.”

We have corrected any deviations from the traditional usage of scientific names: the first use of the name is written out fully and thereafter the genus name is abbreviated as a single letter and only the species name is written out fully. The exception is for readability of the new figures: we adopted the use of a 4 letter code for the species that consists of the first two letters of the genus and species names, e.g. CRBI for Crotaphytus bicinctores. This is defined in the figure legends. If the editor would like us to be consistent across text and figures, we could easily adopt the 4 letter codes in the text as well.

Sincerely,

Eric Routman

References:

Hague MTJ, Routman EJ. Does population size affect genetic diversity? A test with sympatric lizard species. Heredity. 2016 Aug 26;116:92-98.

Micheletti S, Parra E, Routman EJ. Adaptive color polymorphism and unusually high local genetic diversity in the Side-Blotched Lizard, Uta stansburiana. PLoS One. 2012 Oct 25;7(10):e47694.

Submitted filename: Rutherford et al cover letter for revision.docx

Click here for additional data file.

4 Oct 2019

Genetic variation across trophic levels: a test of the correlation between population size and genetic diversity in sympatric desert lizards

PONE-D-19-14658R1

Dear Dr. Routman,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Trishna Dutta

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Dear Dr. Routman,

Congratulations on submitting a well-addressed reply answering the issues raised by the two reviewers and myself. I find that you have addressed a majority of the issues satisfactorily, and I am happy to accept this manuscript for publication in PLOSONE.

My only minor comments are:

(1) I encourage you to insert an inset in Fig1 (map of the study area), and adjust marker transparency so that the underlying collection locations can be visualized.

(2) Page 19, L 368: First word should be Indeed, but is ndeed. Please correct this.

Thank you and Congratulations

Regards,

Trishna Dutta

Reviewers' comments:

22 Nov 2019

PONE-D-19-14658R1

Genetic variation across trophic levels: a test of the correlation between population size and genetic diversity in sympatric desert lizards

Dear Dr. Routman:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Trishna Dutta

Academic Editor

PLOS ONE