Multiple Paternity in Garter Snakes With Evolutionarily Divergent Life Histories.

Many animal species exhibit multiple paternity, defined as multiple males genetically contributing to a single female reproductive event, such as a clutch or litter. Although this phenomenon is well documented across a broad range of taxa, the underlying causes and consequences remain poorly understood. For example, it is unclear how multiple paternity correlates with life-history strategies. Furthermore, males and females may differ in mating strategies and these patterns may shift with ecological context and life-history variation. Here, we take advantage of natural life-history variation in garter snakes (Thamnophis elegans) to address these questions in a robust field setting where populations have diverged along a slow-to-fast life-history continuum. We determine both female (observed) and male (using molecular markers) reproductive success in replicate populations of 2 life-history strategies. We find that despite dramatic differences in annual female reproductive output: 1) females of both life-history ecotypes average 1.5 sires per litter and equivalent proportions of multiply-sired litters, whereas 2) males from the slow-living ecotype experience greater reproductive skew and greater variance in reproductive success relative to males from the fast-living ecotype males despite having equivalent average reproductive success. Together, these results indicate strong intrasexual competition among males, particularly in the fast-paced life-history ecotype. We discuss these results in the context of competing hypotheses for multiple paternity related to population density, resource variability, and life-history strategy. © The American Genetic Association. 2021.

Reproduction requires energetic investment by both females (e.g., number and size of offspring) and males (e.g., territorial defense, mate guarding, sperm size, and quality; reviewed in Schwarzkopf and Shine 1991). Understanding energetic trade-offs in these components of reproduction is the foundation of life-history theory (Roff 1992, Stearns 1992). Therefore, differences in reproductive strategy and cost between the sexes are informative traits to study when unraveling the mechanisms behind life-history evolution (Glaudas et al. 2020, reviewed in Wedell et al. 2006). This is because differential investment in reproduction among individuals in a population can result in some individuals having disproportionately more offspring than others (known as reproductive skew; Clutton-Brock 1989). This increased variance in reproductive success can be either between sexes or among individuals of the same sex and may result in trade-offs with other life-history characters (e.g., a negative relationship between reproductive effort and longevity).

Multiple paternity—when more than one male contributes genetically to a single reproductive bout of a single female—is widespread across the animal kingdom (Taylor et al., 2014). The causes and consequences at the individual and population levels are myriad and likely depend on ecological and life-history contexts. For example, multiple paternity can result either in greater reproductive skew or in a more uniform distribution of males contributing to the next generation. In some vertebrate populations, reproductive skew is exaggerated in promiscuous mating systems due to variance in the occurrence of multiple paternity (e.g., in swordtail fish, Luo et al. 2005) or a small number of males dominating a harem (e.g., in elephant seals, Hoelzel 1999). In other contexts, multiple paternity may reduce—instead of increase—variance in male reproductive success (Pearse and Anderson 2009). At the individual level, fitness benefits of multiple mating can be direct. For example, female Texas field crickets (Gryllus texensis) gain benefits to immune function (Worthington and Kelly 2016) and reproductive output (Worthington et al. 2015) from multiple matings. Fitness benefits of multiple mating may also be indirect and extend across generations (Jennions and Petrie 2000). For example, in dark-eyed juncos (Junco hyemalis), offspring produced from extrapair offspring enjoy higher reproductive success compared with offspring produced by a social pair (Gerlach et al. 2012). And in marbled salamanders, offspring from multiply-sired clutches have higher survival to metamorphosis (Croshaw et al. 2017). Differences in the benefits of multiple matings between males and females can lead to sexual selection and increased reproductive skew, either within or between sexes. Increased reproductive skew can, in turn, decrease the effective population size, and ultimately the ability of a population to respond to selection due to increased inbreeding (Charlesworth 2009). Thus, understanding the specific context and outcomes of multiple paternity sheds light on the interplay between reproductive life-history traits, reproductive strategy, and their evolution.

In snakes, multiple paternity has been identified in nearly every taxa tested and is therefore likely the ancestral reproductive system in snakes (reviewed in Rivas and Burghardt 2005; Uller and Olsson 2008; Voris et al. 2008; Wusterbarth et al. 2010; Jellen and Aldridge 2011; Meister et al. 2012). Past work in snake species has examined the potential drivers of differences among populations in the frequency of multiple paternity (Garner et al. 2002; Prosser et al. 2002; Friesen, Kerns, et al. 2014; Stedman et al. 2016), but we still lack a thorough understanding of how mating systems, life histories, and sexual conflict interact to shape variation in multiple paternity and consequent male reproductive success. In resource-rich environments where most reproductively mature females can dedicate resources to offspring, male–male competition (e.g., combat or sperm competition) may result in relatively fewer males, compared with females, contributing to the next generation (reviewed in Friesen et al. 2020), for examples in snakes, Prosser et al. 2002; Blouin-Demers et al. 2005). In contrast, unpredictable environments and those with limited resources may prompt females to mate with multiple males to increase the genetic diversity of her offspring and thus increase the number of males siring offspring in a population, a strategy generally referred to as bet-hedging (Yasui 1998, 2001; Calsbeek et al. 2007; Wapstra and Olsson 2014). Resource limitation could also decrease the number of males reproducing by increasing variance in male condition and via increased male-male competition.

The populations of western terrestrial garter snakes (Thamnophis elegans) surrounding Eagle Lake in Lassen County, CA, offer a well-established study system to test for differences in multiple paternity and reproductive skew in the context of the evolution of life-history strategies. These discrete populations of snakes are dichotomous in habitat type, found either around the rocky shore of Eagle Lake, with constant food and relatively high levels of predation, or in higher elevation mountain meadows, with fluctuating levels of food and water that are dependent upon annual snow melt (Miller et al. 2011; Sparkman et al. 2013; Gangloff et al. 2020). Divergence in life-history characteristics between the 2 types of habitats has resulted in 2 distinct life-history ecotypes of T. elegans (Bronikowski and Arnold 1999; Addis et al. 2017), which differ in morphology (Manier et al., 2007), behavior (Gangloff et al. 2017), and physiology (Robert and Bronikowski 2010; Palacios et al., 2012; Schwartz and Bronikowski 2013; Gangloff et al., 2015). Most relevant to this study are the striking differences in reproduction between the 2 ecotypes (Sparkman et al. 2007). In lakeshore habitats, snakes having a faster pace-of-life (“L-fast” ecotype), grow quickly, achieve larger asymptotic adult sizes, reach sexual maturity earlier, and reproduce more often and with larger litters relative to slower pace-of-life meadow snakes (“M-slow” ecotype). Mating occurs primarily in spring, but can occur throughout the active summer season (personal observation). Females gestate over the summer and give birth to young on a single day in autumn. As with other species of garter snake, the females are capital breeders and depend on the stored resources from the previous year for reproduction in a current year (Rossman et al. 1986; Gregory and Skebo 1998; Madsen and Shine 1999). During gestation, females must actively thermoregulate to maintain optimal body temperature for embryonic development (Arnold and Peterson 2002; O’Donnell and Arnold 2005) and many stop eating midway through gestation (Bronikowski and Arnold 1999). This is evident in energetic trade-offs, for example pregnant females display lower T-cell proliferative ability and lower counts of white blood cells than non-pregnant females (Palacios and Bronikowski 2017). Thus, reproduction is costly for females because of energetic allocations to vitellogenesis and developing embryos, reduced foraging and ingestion capacity, impaired locomotor ability (Seigel et al. 1987; Van Dyke and Beaupre 2011), and lower adaptive immunity.

Thamnophis elegans is nonterritorial with a polygynandrous mating system and no parental care, and, therefore, no post-parturition cost of provisioning on either sex. Therefore, we expect levels of multiple mating to be high if costs of mating for females are low. If females are not exerting choice for sires (used here to include both mate choice and/or cryptic choice via in utero sperm selection), we would expect the level of multiple paternity to also be high. The only previous study of multiple paternity in this species detected paternity from up to 3 males in a single litter (N = 6 litters; Garner and Larsen 2005). If encounter rates are equivalent between ecotypes and there is no female choice, we predict increased multiple paternity in the L-fast ecotype due to their larger litter sizes simply because numerically, more ovulated eggs yield increased potential for fertilization. That this is not generally seen in reptiles (Uller and Olsson 2008; Jellen and Aldridge 2011) suggests that females are exerting some control over number of sires (either through pre- or postcopulatory mechanisms. As well, a positive correlation between litter size and number of sires could further arise from sexual selection for males to mate with larger, and therefore more fecund, females (Garner et al. 2002; Becher and Magurran 2004; Garner and Larsen 2005; Neff et al. 2008; Jellen and Aldridge 2011), thus increasing the number of sires contributing to the largest litters. Differences in litter sizes but not levels of multiple paternity would suggest that females with smaller litters are employing strategies to maintain levels of multiple paternity despite fewer opportunities for fertilization (i.e., fewer offspring). Using these garter snake populations with divergent life-history ecotypes, we specifically test for differences in the occurrence of multiple paternity, the average number of fathers per litter, male reproductive success, and female reproductive success. To quantify the consequences of these life-history strategies, we compare 3 measures of reproductive skew: between sexes within a population, among individuals within each sex, and among sires within litters. We provide an empirical test of association between life-history traits and mating systems in natural populations.

Materials and Methods

Populations and DNA Sampling

We collected tissue samples from adult snakes from 3 L-fast and 3 M-slow populations surrounding Eagle Lake in Lassen County, CA, from 2006 to 2008 (Figure 1, Table 1). Snakes were sexed, weighed, and measured (snout-to-vent length [SVL] in mm) at the time of capture. In 2006, 56 gravid females (L-fast: N = 30, M-slow: N = 26) were returned to Iowa State University and their litters were born in captivity. Neonates (N = 463 live-born, 12 still-born) were sexed, weighed, and measured (SVL) within 24 h of birth. Additionally, muscle tissue from the tail tip was sampled from each neonate for genetic analyses. Populations are designated by letter (M for M-slow, L for L-fast) and number combinations previously used in Bronikowski and Arnold (1999), with the addition of an L-fast population L4 (Figure 1, Table 1). In addition to these gravid females, we collected a tissue sample from additional adult males and females from each of the 6 populations to estimate population-level genotype frequencies (N = 213 additional males and females).

Figure 1.

Map of population locations of Thamnophis elegans around Eagle Lake, Lassen County, CA, with M and L indicating replicate Meadow-slow (blue) and Lakeshore-fast (red) populations. Numbering from Bronikowski and Arnold 1999. Map created by Katie Fetting using ArcGIS 10.3 (ESRI, Redlands CA).

Table 1.

Sample sizes and genetic diversity measures computed from the 8 microsatellite loci in Thamnophis elegans

Habitat	Sample sizes					Diversity measures
Population	NAdults	NMales	NFemales	NLitters	NOffspring	GD (SD)	HO (SD)	NA (SD)	PA	FIS	PIC (SD)	PEXC
L-fast
Total/Avg	104	33	71	30	323	0.71 (0.05)	0.71 (0.04)	6.2 (2.2)	6		0.68 (1.5)
L2	56	27	29	13	154	0.66 (0.04)	0.65 (0.02)	7.0 (2.4)	1	0.008	0.69 (1.4)	0.003
L3	7	0	7	4	26	0.77 (0.05)	0.74 (0.06)	5.0 (1.4)	1	0.045	0.67 (1.7)	0.004
L4	41	6	35	13	143	0.71 (0.04)	0.69 (0.03)	6.6 (2.4)	4	0.026	0.67 (1.4)	0.004
M-slow
Total/Avg	165	87	78	26	140	0.71 (0.05)	0.71 (0.02)	7.8 (2.9)	4		0.64 (1.4)
M1	63	39	24	8	52	0.71 (0.04)	0.69 (0.02)	7.9 (2.8)	1	0.019	0.65 (1.4)	0.003
M2	52	25	27	9	33	0.70 (0.05)	0.68 (0.02)	7.1 (3.0)	2	0.026	0.61 (1.4)	0.003
M3	50	23	27	9	55	0.72 (0.05)	0.70 (0.02)	7.9 (2.4)	1	0.031	0.67 (1.4)	0.003
Overall Total/average	269	120	149	56	453	0.71 (0.05)	0.71 (0.03)	7.0 (2.6)	10		0.66 (1.5)

GD, Nei’s unbiased genetic diversity; HO, observed heterozygosity; NA, average number of alleles; PA, number of private alleles; FIS, system of mating inbreeding coefficient; PIC, polymorphic information content, a measure of the information content of the loci for paternity analysis related to expected heterozygosity; PEXC, is the probability of not excluding a random candidate male as the father if the mother genotype is known (calculated in cervus).

Genotyping

We purified DNA from tissue samples using a standard salt extraction method modified from Sunnucks and Hales (1996). Purified DNA was diluted to 25 ng/µL to achieve consistent genotyping results. Eight highly variable microsatellite loci reported in previous publications (Garner et al. 2004; Manier and Arnold 2005; McCracken et al. 1999; Prosser et al. 1999) were modified slightly for this study (Table 2). An M13-tail was added to one primer from each locus to allow tagging with M13 fluorescently labeled primers (Boutin-Ganache et al. 2001). The loci were organized into 3 multiplexes for amplification, using a Type-it Microsatellite PCR kit (Qiagen, Venlo, The Netherlands). All loci were amplified in 7-µL reactions using 25 ng DNA, 0.03–0.65 nM of each locus-specific reverse and forward primer, 0.5–0.6 nM of the M13 primer labeled with a fluorescent dye, and 3 mM MgCl2 (in Type-it Multiplex PCR Master Mix, Qiagen). Amplification was performed with a Bio-Rad iCycler (Bio-Rad Laboratories, Hercules, CA) at 95 °C for 5 min to denature, and then a cycle repeated 30 times of 95 °C for 30 s, 57 °C for 1 min 30 s, and 72 °C for 30 s, and lastly an extension of 30 min at 60 °C. After amplification, 0.5 µL of each sample from each multiplex (1.5-µL total PCR product per sample) was mixed with 8.5-µL H2O to create a 1:20 dilution. 1.5 µL of each of these diluted samples was electrophoresed at the Iowa State University DNA Facility using the ABI 3100 Genetic Analyzer (Applied Biosystems, Foster City, CA). To minimize scoring biases and errors, we used GeneMapper (Applied Biosystems) to define allelic bins. The electropherograms were scored automatically by GeneMapper and subsequently verified or corrected by eye, independently by 2 researchers (M.B.M., T.S.S.).

Table 2.

Properties of the microsatellite loci used in this study

Locus	Multiplex	Final concentration (nM)	Allele size (bp)	Total alleles	Avg HO	Avg HE	Avg PIC	Null alleles	Reference
Nsμ10	3	F = 0.03 R = 0.45	140–162	7	0.626	0.618	0.54	0	Prosser et al. (1999)
Nsμ2	2	F = 0.03 R = 0.45	181–185	5	0.612	0.630	0.55	0	Prosser et al. (1999)
Nsμ3	3	F = 0.03 R = 0.35	167–219	14	0.849	0.859	0.82	0	Prosser et al. (1999)
TS10	1	F = 0.03 R = 0.35	147–187	11	0.815	0.814	0.77	0	Manier and Arnold (2005)
TS2	2	F = 0.03 R = 0.15	136–169	10	0.669	0.678	0.63	1	McCracken et al. (1999)
Te1Ca18	1	F=0.03 R = 0.28	96–131	11	0.633	0.615	0.57	0	Garner et al. (2004)
Te1Ca2	1	F = 0.03 R = 0.50	232–266	13	0.760	0.781	0.74	0	Garner et al. (2004)
Te1Ca3 a	2	F = 0.04 R = 0.65	91–123	9	0.579	0.692	0.63	3	Garner et al. (2004)
Total/average				80	0.693	0.711	0.66	0.5

For the concentration used, F refers to the forward primer, and R refers to the reverse primer. HO, observed heterozygosity; HE, expected heterozygosity; PIC, polymorphic information content. See text for details.

aExcluded from paternity analysis due to rejection of Hardy–Weinberg equilibrium and presence of null alleles (see text for details).

Loci Verification and Population Genetic Measures

The genotypes of field-caught adults (N = 269, including the N = 56 pregnant females) were analyzed in the program MicroChecker (Van Oosterhout et al. 2004) to estimate null allele frequencies and to test for large allele dropout and the Wahlund effect. Pairwise estimates of linkage disequilibrium were calculated in GenePop (Raymond and Rousset 1995). We tested for deviations from Hardy–Weinberg equilibrium (HWE) for each locus by population combination in Arlequin v.3.11 (Excoffier et al. 2005) using an Exact test with Markov chain lengths of 10 000 with 1000 dememorization steps. Significance values were corrected for multiple testing using a sequential Bonferroni correction (Rice 1989). The Te1Ca3 locus was dropped from the subsequent population structure and paternity analyses due to the presence of null alleles and significant deviation from HWE (see Results).

For each population, estimates of genetic diversity were averaged over all loci. These estimates included observed and expected heterozygosities (Avg HO and Avg HE, respectively), unbiased genetic diversity (Nei 1987), the size range of the alleles in base pairs (bp), average number of alleles, total number of alleles, and number of private alleles. Polymorphic information content (PIC; Bolstein et al. 1980) and paternity exclusion power were calculated in Cervus 3.0 (Kalinowski et al. 2007). We used Arlequin (Excoffier et al. 2005) to calculate pairwise FST values to estimate genetic distance among the populations and Fstat (Goudet 2001) to estimate the inbreeding coefficient (FIS).

Paternity Analysis

We successfully assigned paternity to 411 of the 463 (89%) of the live-born offspring. The majority of offspring were genotyped from each litter (minimum 58% of offspring in a given litter, mean 86.3%). Restricting analyses to litters where 100% of babies were sampled (N = 22 litters) does not qualitatively change results (analyses not shown), therefore we present the results with this 58% cutoff. The minimum number of fathers that contributed to each litter was determined using Colony (Wang 2004). This program estimates the likelihood that offspring are related to each other (i.e., are siblings), as well as the likelihood that any of the genotyped adults are the mothers and fathers of the offspring. For litters where not all offspring were successfully genotyped, the actual number of fathers may be higher. The data set included genotypes of known mother/offspring relationships and genotypes for adult males serving as potential fathers, without any designation of relatedness. Each population was analyzed separately using population-specific allele frequencies calculated from the adult samples in that population. We estimated rates of errors (i.e., null alleles, mutations, genotyping error) and mismatches between mother and offspring alleles using MicroChecker. Using the allele frequencies, as well as the known mother–offspring relationships, Colony generated expected genotypes of the males contributing to each offspring in a litter with a 95% confidence level. Colony settings were as follows: polygamous mating system, dioecious and diploid species, medium run length, full-likelihood analysis method, medium likelihood precision, and no sibship size prior. To assess the accuracy of our results and to identify the effects of different error rates on the number of fathers assigned to a litter, simulations were conducted in Colony utilizing varying genotyping error rates (0–2%), as well as different run lengths (“medium” to “very long”) and precision levels (“medium” to “high”). Our paternity assignments were robust to these variations.

Multiple Paternity and Reproductive Success

We used a chi-square test to assess differences in the levels of multiple paternity between ecotypes (binary response variable: Y/N). Because the possible number of fathers contributing to a litter ranged from 1 to 4 in this study (see Results), to test for differences in the number of males contributing to each litter between ecotype we utilized a multiple ordinal logistic regression, including in the model population nested within ecotype and the number of offspring sampled within each litter. To account for differences in litter size between ecotypes (L-fast > M-slow), we conducted a permutation analysis and restricted the litter size for L-fast dams to the median litter size of the smaller M-slow litters (median = 5). For each L-fast dam, we made pseudo-litters by sampling her offspring with replacement to fill this litter of 5, keeping dam and sire associations, and thus making no assumptions about sperm precedence or competition. We repeated this procedure for the sampled population 999 times and calculated the average number of sires for each truncated L-fast litter. We then compared the distribution of paternity in the populations of pseudo-litters with the observed values. If the observed value of sires per litter differs from the number of sires of truncated litters, this would suggest that females may actively employ strategies to alter levels of multiple paternity beyond a sampling effect due to differences in litter sizes.

Reproductive success was measured in females as total live litter size and was measured in males as the number of offspring each male sired. Male reproductive success is necessarily limited to the offspring sampled and may therefore not represent the total number of offspring in the population fathered by a single male. We used ANCOVAs to test the effects of ecotype and population nested within ecotype on female and male reproductive success. The analysis of female reproductive success also included the covariate of body size (SVL). Statistical analyses were performed in SAS 9.1 (SAS Institute, Cary, NC).

Estimates of Reproductive Skew

To evaluate how the reproductive traits within each ecotype could affect the distribution of reproductive success among different groups of individuals, we calculated 3 estimates of reproductive skew within each population or ecotype: (1) To test for differences in the distributions of reproductive success, we employed a Kolmogorov–Smirnov test (Massey 1951) at 3 levels: between males and females within ecotypes, among males between ecotypes, and among females between ecotypes. (2) To test for skew of reproductive success among breeders within a sex we used an Index of Variability (“I-V”; Araki et al. 2007), the variance in reproductive success among individuals of the same sex divided by the mean reproductive success for that sex in a population. Based on the expected Poisson distribution, a value of one indicates that reproductive success was evenly spread among the individuals within that sex in that population, whereas an inflated value indicates one or more members of that sex had disproportionately higher reproductive success. (3) We used a skew estimate from Neff et al. (2008) to assess the distribution of sires contributing to each litter (within-litter skew) as an indicator of sperm selection and/or cryptic female choice. Values significantly different from zero, as assessed with a one-tailed t-test, indicate a skew among sires contributing to each litter.

Results

Of the 168 pairwise within-population loci comparisons, 12 were found to have linkage disequilibrium between pairs of loci, but only one (between TS10 and Te1Ca18 in M2 population) was significant after sequential Bonferroni correction. Because these loci only showed linkage in one population, we assume it is an artifact and not actually due to physical linkage. All loci in all populations were in HWE except for Te1Ca3, which had significantly reduced heterozygosity relative to the expected in both L4 and M2 populations. MicroChecker predicted the presence of null alleles (i.e., an allele that is not detectable due to mutation) in loci within 4 populations used in this study, with 3 of the 4 instances attributed to Te1Ca3. Therefore, we did not include the Te1Ca3 locus in the analyses of parentage in Colony or population structure in Arlequin. The TS2 locus showed evidence of null alleles in only one population (M1), and this error estimate was incorporated into analysis for this population at this locus in Colony.

The populations had similar measurements of genetic diversity, inbreeding coefficients (FIS), and PIC; thus, the loci were equally suitable for determining paternity across all the populations (Table 2). The global FST for these populations was 0.068 (P < 0.0001); of the 15 pairwise FST values, 12 were significant at P < 0.05 and 11 were still significant after sequential Bonferroni corrections (Table 3). This genetic differentiation among populations suggests it is unlikely that males from one population could be potential fathers in other populations. This is supported by our mark-recapture data over the past 40 years, where we have captured only a few female migrants and no male migrants (A.M.B., unpublished data).

Table 3.

Population genetic structure among Thamnophis elegans populations

	L2	L3	L4	M1	M2
L3	0.001
L4	0.019*	0.00
M1	0.076*	0.046*	0.075*
M2	0.067*	0.042*	0.072*	0.003
M3	0.065*	0.039*	0.069*	0.008*	0.007

Pairwise FST values (calculated using 7 loci) below diagonal. Asterisk indicate significant values after sequential Bonferroni correction. See Supplementary Table S1 and Supplementary Figure S1 for additional structure information.

We genotyped an average of 86.3% of the offspring from the 56 litters (range: 58–100% genotyped, combining both ecotypes). Multiple paternity was detected in all 6 populations. The maximum number of fathers identified in a single litter by Colony was 4, which was found in 2 litters, one M-slow (population M2) and one L-fast (population L4). A total of 68 sires contributed to these 56 litters, as identified by unique paternal genotypes contributing to offspring within each litter. Only 6 sires had known identities from being caught in the field (i.e., 6 of the 120 field-caught males included in the Colony analysis). Of the 411 analyzed offspring, 34 were assigned to these 6 males with known identities.

The number of multiply-sired litters did not differ between the ecotypes: 50% of litters in L-fast, 38% of litters in M-slow ( = 0.36, P = 0.55; Table 4). The mean (± SD) number of sires for each ecotype was 1.60 (range: 1–4) sires per litter for L-fast populations and 1.46 (range: 1–4) sires per litter for M-slow populations, and ecotypes did not differ (multiple ordinal logistic regression: Wald’s = 0.024, P = 0.88; Tables 4 and 5). Neither the sampled litter size nor maternal size affected the number of sires in each litter and the lack of difference between ecotype holds if we remove these covariates from the model. The mean number of sires per litter within individual populations ranged from 1.25 to 1.77, and did not differ significantly among populations (multiple ordinal logistic regression: Wald’s = 3.63, P = 0.46; Tables 4 and 5). The permutations restricting the large L-fast litters to smaller M-slow sizes showed that, L-fast litters would not differ in number of sires form observed M-slow litters (i.e., 1.39 vs. 1.46, P = 0.136; Figure 2).

Table 4.

Components of reproductive skew among populations of Thamnophis elegans for males and females

Ecotype	Population	Number of breeders	Reproductive successa	% multiply-sired littersa	I-Va	Skew within litter	Observed sires per litter	Operational sex ratio
L-fast
	Dams	30
	L2	13	11.84	38%	0.91		1.46	0.62
	L3	4	6.50	50%	1.49		1.50	0.86
	L4	13	11.00	62%	0.45		1.77	0.25
	Dam Avg.		10.77	50%	0.98		1.60	0.55
	Sires	39				0.12***
	L2	16	8.69		4.20	0.07*
	L3	5	4.20		2.55	0.06
	L4	18	6.78		2.72	0.17**
	Sires Avg.		7.23		3.59
M-slow
	Dams	26
	M1	8	6.50	25%	1.10		1.25	0.40
	M2	9	3.67	44%	0.27		1.67	0.68
	M3	9	6.11	44%	0.22		1.44	0.97
	Dam Avg.		5.38	38%	0.82		1.46	0.71
	Sires	29				0.07**
	M1	8	6.00		0.81	0.10
	M2	10	3.00		0.89	0.04
	M3	11	4.64		0.96	0.08*
	Sires Avg.		4.45		1.15

Number of breeders in the study for each population; least-squares means of average male and female reproductive success (see text for details); percent of litters with multiple paternity; Index of Variability (I-V) among breeding individuals of the same sex in a population; skew among sires within a litter; average number of observed number of sires; operational sex ratio (values less than 1.0 are female biased). Asterisk denote significant estimates/comparisons (*P < 0.05, **P < 0.01, ***P < 0.001).

aNote that ecotype averages are indices calculated with values for all individuals, not arithmetic means of population values.

Table 5.

Analysis results for number of sires per litter, female reproductive success, male reproductive success, and female size in Thamnophis elegans populations

Source of variation	Number of sires per litter (>58% sampled; N = 56)	Female reproductive success	Male reproductive success	Female size (SVL)
Ecotype
Test statistic	χ12 = 0.024	F 1,49 = 3.31	F 1,62 = 3.53	F 1,50 = 94.26
P-value	0.88	0.075	0.065	<0.0001***
Litter size sampled
Test statistic	χ12 = 1.76	—	—	—
P-value	0.19	—	—	—
Population (ecotype)
Test statistic	χ12 = 3.63	F 4,49 = 4.90	F 4,62 = 1.91	F 4,50 = 0.87
P-value	0.46	0.0021**	0.12	0.49
Body size
Test statistic	χ12 = 0.0038	F 1,49 = 6.66	—	—
P-value	0.95	0.013*	—	—

“—” indicates effect not included in model (see text for details). Asterisk denote significant factors (*P < 0.05, **P < 0.01, ***P < 0.001).

Figure 2.

Distribution of average number of sires per litter for L-fast populations from simulations restricting L-fast litter sizes to M-slow median litter size of 5. Frequency is the number of times the mean number of sires was observed out of 999 iterations. The red line indicates the observed average number of sires for M-slow litters (1.46 sires per litter). The mean number of simulated L-fast fathers was not significantly different from observed M-slow values (P = 0.136).

Female reproductive success differed between ecotypes, with L-fast females having larger litters on average than M-slow litters (raw data quartiles: L-fast = 9, 11, 12, 17; M-slow: 4, 5, 6, 12). This is primarily due to body size differences, as ecotypes differed marginally in female reproductive success after correction for body size (LSmean ± SE: L-fast = 8.8 ± 0.6 offspring, M-slow = 6.8 ± 0.7, Table 5, Figure 3A; body size LS mean ± SE for M-slow females 473 ± 8.2 mm and L-fast females 590 ± 8.8 mm; Table 5; see also Bronikowski and Arnold 1999 and Sparkman et al. 2007). Finally, litter size differed among populations, and it is noteworthy that female reproductive success for one L-fast population was equal to that of one of the M-slow populations (Table 4). Similar to females, reproductive success of males, measured as the total number of offspring each male sired, contrasted between ecotypes with a trend for higher success in the L-fast ecotype (LSmean ± SE: 6.6 ± 0.8 offspring per male) compared with the M-slow ecotype (4.5 ± 0.8 offspring per male; Table 5, Figure 3B).

Figure 3.

Distribution of reproductive success of sires and dams in wild populations of Thamnophis elegans: comparing ecotypes within each sex (A and B) and comparing sexes within each ecotype (C and D). For males, the number of offspring signifies the total number of offspring that can be attributed to a single male of all of the offspring sampled within a given population. For females, the number of offspring is her liveborn offspring. (A) Significant distribution differences between M-slow and L-fast dams (P < 0.0001); (B) significant distribution differences between M-slow and L-fast sires (P = 0.010). (C) Significant distribution differences between sexes within the L-fast ecotype (P = 0.001), whereas sexes did not differ in the M-slow ecotype (D).

When we compared the distributions of reproductive success within each sex between ecotypes, we found that the L-fast males and females had significantly right-shifted distributions (i.e., more individuals exhibiting high reproductive success) compared with M-slow males and females (Kolmogorov–Smirnov test, males: D = 0.40, P = 0.010; females: D = 0.77, P < 0.0001; Figure 3A and B), which is in agreement with the trend for differences in mean reproductive success within each sex. Within ecotypes, the distribution in reproductive success differed between males and females in L-fast populations (Kolmogorov–Smirnov test, L-fast males and females: D = 0.47, P = 0.001) but not in M-slow populations (D = 0.23, P = 0.49; Figure 3C and D). At the ecotype level, the Index of Variability (I-V), which tests the variance among individuals of the same sex within each population, was inflated above one in males from both ecotypes, but not in females (Table 4). This indicates that males of both ecotypes had greater among-individual variation in reproductive success than if reproduction followed a Poisson distribution, with L-fast males demonstrating much greater skew than M-slow males. In females, reproductive success was more evenly spread across individuals. At the population level, the test for sire skew within litters indicated that 3 populations were significantly skewed (2 L-fast populations and 1 M-slow population), which means that within a given litter, some males fathered more offspring than others. However, the 2 ecotypes were not significantly different from one another in within-litter reproductive skew (Table 4).

Discussion

We found that multiple paternity was common in all 6 sampled populations of both life-history ecotypes. Although L-fast females have larger litter sizes and reproduce more frequently than M-slow females, both the number of sires within a litter and the frequency of multiply-sired litters did not differ between ecotypes. Based on life-history and morphological differences between the ecotypes, we had predicted the L-fast ecotype to have a higher frequency of multiple paternity due to a positive association between litter size and the potential for more sires to contribute to a litter. However, the generally higher measures of reproductive skew among sires (both within and among litters in the L-fast ecotype) suggest that M-slow and L-fast individuals are using different tactics to achieve reproductive success. At this time, we cannot assess whether this is due to courtship behavior of males or receptivity behavior of females. Furthermore, our estimates of reproductive skew suggest the hypothesis that M-slow dams are actively engaging in reproductive strategies to increase the level of multiple paternity within litters. However, our permutation analysis suggests that the mean number of sires contributing to litters would not differ between ecotypes if litters sizes were equivalent. Both female strategy and litter size likely contribute to observed patterns of multiple paternity and reproductive skew between the ecotypes.

Differences in M-slow female mating strategy, including both pre- and postcopulatory mechanisms, may compensate for small litter sizes with the outcome of similar multiple paternity levels and number of sires per litters as seen in L-fast populations. This increase in multiple paternity in M-slow females could result from several non-mutually exclusive possibilities: increased mating events and decreased precopulatory female choice relative to L-fast females, increased intrasexual competition among males resulting in greater variance in male reproductive success (and therefore more males producing no offspring), increased cryptic female choice, or potential female bet-hedging in the stochastic habitats of the M-slow populations. Work in other systems demonstrates the potential role of the environment in shaping multiple paternity patterns. For example, the increased level of relative multiple paternity in M-slow populations could result from shorter breeding seasons in the cooler meadow habitats (e.g., in watersnakes: Prosser et al. 2002). On the contrary, the relatively high levels of multiple paternity in M-slow snakes are counter to studies in birds and fish that found more extrapair matings in shorter-lived species and increased levels of multiple paternity in populations experiencing higher predation (Arnold and Owens 2002; Neff et al. 2008). This highlights the need to incorporate both life-history patterns and a characterization of environmental conditions in revealing the drivers of multiple paternity.

Multiple Paternity in a Life-History Context

Previous studies on multiple paternity in garter snakes found positive relationships between litter size and number of sires within a litter (Garner et al. 2002; Garner and Larsen 2005), but this pattern is not universal in snakes (Blouin-Demers et al. 2005; Uller and Olsson 2008). In contrast to these previous studies in garter snakes, we found no association between the number of sires with either larger dams or larger litter size (Table 5). Furthermore, the results of our simulations show that differences in female reproduction between the ecotypes are due to difference in both strategy and litter size. Decreasing the litter size of L-fast dams does not affect the level of multiple paternity. That we did observe increased reproductive skew in both L-fast males and females compared to M-slow males and females suggests that the relationships among morphology, life-history, and multiple paternity is influenced by ecological or physiological variables (Uller and Olsson 2008). We speculate that mating with multiple males confers greater indirect benefits to females and reduces male reproductive skew in M-slow populations, though future work is needed to test these hypotheses directly. Previous work in reptiles points to fitness benefits for females with multiple mates, though the mechanisms conferring these benefits have yet to be thoroughly explored (Olsson et al. 1997; Madsen et al. 2002, 2005; Blouin-Demers et al. 2005). The long-lived females of the M-slow populations may enjoy any combination of possible benefits of multiple matings: “trading-up” (i.e., mating with a successive male after already mating, if the new male appears to be of higher quality); sperm competition (Parker 2020); cryptic female choice; fitness increases through the production of sons that achieve increased fitness (“sexy sons”); and bet-hedging (i.e., mating with several males in order to increase the genetic diversity of offspring; reviewed in Uller and Olsson 2008; Wapstra and Olsson 2014; Friesen et al. 2020). Bet-hedging, in particular, has been found to be a viable explanation for multiple paternity in either small or environmentally variable populations (Yasui, 1998, 2001; Calsbeek et al. 2007; Sarhan and Kokko, 2007; Makinen et al. 2007), such as the M-slow populations in this study (Miller et al. 2011). Additionally, there is some evidence for cryptic female choice in garter snakes (Friesen, Kerns, et al. 2014). However, female garter snakes can store sperm for at least a year (Friesen, Mason, et al. 2014), which means that apparent female choice within a season is not discernable from among-season sperm storage and usage (Friesen, Kerns, et al. 2014; Friesen et al. 2016). Moreover, males deposit significantly fewer sperm with successive matings (Friesen, Uhrig, et al. 2014). Future work in garter snakes can be directed toward specific tests of these contrasting, but not exclusive, hypotheses of sperm competition versus female choice.

Multiple Paternity in a Demographic Context

In the absence of strong evidence to show that females benefit from mating multiple times, levels of mating may simply result from population density and encounter rate (Uller and Olsson 2008, Wells et al. 2017). Although we did not have a reliable method for determining the population density or encounter rates at the time of mating, we were able to calculate operational sex ratios during the active season to estimate one aspect of population demographics. In both ecotypes, more adult females than males were observed over 4 years of capture/mark/recapture fieldwork, but in L-fast populations, this female bias was far more pronounced. The average male to female operational sex ratio was 0.55 in L-fast populations and 0.71 in M-slow populations, indicating a strong bias toward female captures (Table 4). In 2006, the year in which the gravid females in this study were collected, our estimates are similar to the longer-term averages (L-fast OSR = 0.45, M-slow OSR = 0.65). These estimates of sex bias in the adult populations may result from males suffering higher mortality rates associated with more frequent long-distance movements, as identified in other garter snake systems (e.g., Bonnet et al. 1999; Shine et al. 2001a). The female-biased sex ratio may also be compensated for by the temporal distribution of mates in the garter snake mating system, which provides a strongly male-biased sex ratio during the primary time of mating. For example, in many garter snake populations, males congregate around hibernacula and compete by scramble competition for emerging females. The males stay at the hibernaculum for up to several weeks, while the females disperse in a few days (e.g., Shine et al. 2001b).

Conclusions and Future Directions

Our findings suggest that other factors beyond operational sex ratios and male/female encounter rates might underlie the unexpected equivalence in degree of multiple paternity found in these divergent populations of garter snakes. Sperm competition, female reproductive strategies (i.e., bet-hedging, cryptic female choice), behavioral differences in males from L-fast populations (Gangloff et al. 2017), or a combination of any of these mechanisms might offer a more complete representation of the complex interactions that contribute to multiple paternity in these populations. In a sex-specific transcriptome scan of these populations, loci for proteins involved in sperm/egg interactions were identified as being highly variable and potentially under diversifying selection, which could suggest on-going sperm competition (Schwartz et al. 2010). Thus, sperm competition (i.e., sperm morphology, energetics, and swimming performance) is a potential determinant of the skew observed within litters, as has been suggested in other studies (Olsson and Madsen 1998; Friesen, Kerns, et al. 2014; Friesen, Mason, et al. 2014). Future studies utilizing a genome-wide approach could also point to fruitful candidate loci for discerning mechanisms of intrasexual competition and targets of selection in males (see also Levine et al. 2015). Controlled matings in a laboratory setting might also be helpful in determining if mating order has any implications for the reproductive success of males and if male quality (i.e., body condition, sperm quality) has any impact on the survival of the offspring they sire. Additionally, observations of mating behaviors in the field might provide insight into whether females are exerting choice in the form of mating preferences, or if male behavioral aspects (e.g., predator avoidance, mating coercion) are more important in determining the number of matings occurring in populations (e.g., Clark et al. 2014; Lind et al. 2016). Given the prevalence of observed multiple paternity and the known life-history differences among these populations, this system of garter snakes provides a model system to explore these important questions.

Click here for additional data file.