Introduction

A central concept in life history theory deals with the relative allocation of energy or resources that organisms devote to reproduction to maximize the production of offspring (^{Hirshfield & Tinkle, 1975}; ^{Shine & Schwarzkopf, 1992}). However, as the available energy is limited, this investment in reproduction (aka reproductive effort) is governed by a balance between reproductive expenditure, on the one hand, and the survival and future fecundity of the organisms, on the other (^{Shine & Schwarzkopf, 1992}; ^{Williams, 1966}). For squamates, relative clutch/litter mass (RCM, RLM) has been widely used as a measure of female reproductive effort because most species lack parental care, and so energy allocation can be summarized in terms of the relative proportion of mass offspring production accounts for (^{Bastiaans et al., 2013}). During reproduction nutrients that otherwise might be directed to body growth are used instead for developing eggs and embryos, thus reducing the potential of success for future reproductive events (^{Landwer, 1994}; ^{Shine & Schwarzkopf, 1992}). Furthermore, gravid/pregnant individuals might be more susceptible to predation due to the physical burden it represents and to the increase in exposure during maternal basking (^{Miles et al., 2000}; ^{Swain & Jones, 2000}). Factors known to influence reproductive effort in lizards include phylogeny, morphology, foraging and predator-escape strategies, habitat use, resource availability, and even parasite load (^{Hirshfield & Tinkle, 1975}; ^{Sorci et al., 1996}; ^{Uller & While, 2015}; ^{Vitt & Price, 1982}). Moreover, sexual size dimorphism might evolve in lizards in response not only to intraspecific competition among males, but also in cases in which a larger size in females allows a higher clutch/litter size, or even as a way to reduce intersexual competition through resource partitioning (^{Cox et al., 2007}).

Whereas selection for different levels of RCM or RLM is dependent upon adult mortality (due to the physical and physiological burden gravidity or pregnancy imposes to females), other reproductive characteristics, such as clutch/litter size and mass or egg/neonate size are shaped by the progeny mortality (^{Messina & Fox, 2001}; ^{Vitt & Price, 1982}). Indeed, the probability of parents having at least 1 surviving offspring increases with clutch or litter size (^{Messina & Fox, 2001}); also, it is assumed that larger offspring are under less risk of predation (^{Uller & Olsson, 2010}). These characteristics tend to be positively correlated to body size, with larger size females producing either more eggs/neonates or fewer larger ones (^{Uller & While, 2015}; ^{Vitt & Price, 1982}).

Whatever their reproductive strategy, investment is usually higher in viviparous than in oviparous species because of their longer gestation times and extended physical and physiological burden (^{Dupoué et al., 2015}; ^{Recknagel & Elmer, 2019}). Because of this, viviparous lizards are usually limited to no more than 1 brood per year and have lower litter sizes, both of which translate in a reduced reproductive output (^{Meiri et al., 2011}). However, relatively few studies on lizard reproductive investment have been conducted on viviparous species (^{Bastiaans et al., 2013}; ^{Recknagel & Elmer, 2019}; ^{Rodríguez-Romero et al., 2005}), this in spite of the higher extinction risk they face due to climate change (^{Sinervo et al., 2010}). Even fewer studies have investigated the relationship between pregnancy and body condition in lizards, despite their potentially serious implications in terms of costs of reproduction (^{Boretto et al., 2014}; ^{Itonaga et al., 2012}). Given that the estimation of life-history parameters is useful to evaluate the intrinsic vulnerability of reptiles to population declines and extirpations (^{Govindarajulu et al., 2011}), information on reproductive investment of live bearing lizards could be useful for conservation purposes. Within the genus Plestiodon, viviparity has evolved once, in the brevirostris group (^{Watson et al., 2014}), a clade comprised by 15 currently recognized species distributed in Mexico and part of Central America (^{García-Vázquez et al., 2021}). Most of them are highland inhabitants and share a viviparous mode of reproduction (^{Goldberg, 2002}; ^{Feria-Ortiz et al., 2007}; ^{López-Hernández et al., 2008}). Here, we present information on the reproductive ecology of 2 species of this group: an undescribed species previously assigned to Plestiodon brevirostris (1 population from Tlaxcala) and P. dugesii (2 populations from Jalisco). Our goals were to determine whether the studied species and populations follow different strategies of reproductive investment, the existence of sexual size dimorphism and trade-offs in reproduction, as well as to elucidate the underlying causes of these patterns.

Materials and methods

The undescribed species belongs to a taxon previously assigned to Plestiodon brevirostris, which is now known to be a species complex comprised of 4 lineages, 1 distributed in the states of Tlaxcala and Puebla; a second one, in Guerrero, and 2 more from Oaxaca (^{Pavón-Vázquez et al., 2018}). However, none of these taxa have been so far properly described. Given that the type locality of P. brevirostris is indicated by ^{Günther (1860)} as “Oaxaca (Mexico)”, and as the studied population belongs to the Tlaxcala-Puebla lineage, we will refer hereafter to this taxon as Plestiodon sp.

Plestiodon sp. is a skink up to 75.5 mm snout-vent length (SVL) (see supplementary data 5 in ^{Pavón-Vázquez et al. [2018]}). It has an asynchronous fall reproductive cycle, with vitellogenesis occurring between June and October, and ovulation and fertilization in November, so that gestation extends from then until spring (^{Chávez-López, 2012}). The population of Plestiodon sp. studied was from La Malinche, a 4,461 m dormant volcano located in the states of Tlaxcala and Puebla, in central eastern Mexico. Fieldwork at this locality was carried out on the eastern slope of the volcano (19°14’05” N, 97°56’24” W; 2,700 to 3,000 m elevation) in an abrupt transition zone between the Pinus-Quercus forest and induced grasslands. Climate at this altitudinal level is temperate sub-humid with day temperatures ranging from 12 to 18

°C throughout the year (^{Villers-Ruiz, 2006}). Despite its status as a protected natural area, La Malinche has been heavily affected in the last years by land use changes, illegal logging, and forest fires (^{Piñaza-Soto & Hernández- Hernández, 2011}), and more recently by the subterranean channelling of some natural streams, which together have eroded the microhabitat opportunities for Plestiodon sp.

Plestiodon dugesii is a skink of up to 69 mm SVL, distributed in part of the states of Guanajuato, Michoacán and Jalisco in western Mexico, at elevations above 2,000 m. (^{Dixon, 1969}; ^{Flores-Villela & Santos-Barrera, 2007}). Information about the biology of this species is scarce. However, it is under special protection according to the Mexican law NOM-059-SEMARNAT-2010 (^{Semarnat, 2019}), and is included as vulnerable in the Red List of the International Union for Conservation of Nature (IUCN) (^{Flores-Villela & Santos-Barrera, 2007}). We sampled 2 populations of this species, both from the state of Jalisco. The first one was from the surroundings of the cliff El Tecolote, in the municipality of Mazamitla (19°54’43” N, 102°59’30” W; 2,400 m elevation); the second population was from the locality of Ferrería de Tula, in the municipality of Tapalpa (20°03’46” N, 103°42’39” W; 2,500 m elevation). Both study sites have a temperate sub-humid climate with summer rains (^{Ruiz-Corral et al., 2012}), with a mean daytime temperature of 19.6 °C (Mazamitla) and 16 °C (Tapalpa). Although vegetation type at these sites corresponds to pine-oak forest, we primarily encountered lizards using rotten logs at Mazamitla and rocks at Tapalpa, even though both microhabitat types were present at both localities.

We conducted samplings between March and April 2013 and 2014; however, for Plestiodon sp. we only obtained data for 2014. Fieldwork consisted in the active search and capture by hand of lizards in potential microhabitats, such as leaf litter and under rotten logs and rocks. We focused on pregnant females, which we recognized by noticing abdominal distension and through palpation (^{Shine, 1980}), but during the same samplings we also collected non-pregnant females and males for in situ measurements to make comparisons on body size dimorphism (see below). We sexed individuals through manual hemipenial eversion.

We temporarily placed caught pregnant females in separate plastic containers (20 × 12 × 12.5 cm) to which we drilled holes for aeration. Then we transported the lizards to nearby laboratory facilities to conduct the captivity period: the Laboratorio de Herpetología II of the Instituto de Biología of the Universidad Nacional Autónoma de México (Mexico City), in the case of Plestiodon sp., and the Centro de Estudios en Zoología of the Universidad de Guadalajara (Zapopan, Jalisco) in the case of P. dugesii. Once there, we housed them in larger individual plastic enclosures (43.18 × 16.51 × 28.3 cm), with peat moss as substrate and pine bark as shelter and provided them with water and food (mealworm larvae and domestic crickets) ad libitum. Based on the available information of the thermal biology of other Plestiodon species (^{Youssef et al., 2008}), the enclosures were kept at 25-30 °C and under natural photoperiod. For each female, we recorded SVL and axilla- groin length (AGL), as these morphometric variables are known to influence litter or offspring characteristics in lizards (^{Griffith, 1994}; ^{Recknagel & Elmer, 2019}), using a digital calliper (± 0.01 mm; AutoTec®). We checked, at least twice a day, every container until parturition, and immediately weighed the female using a digital scale (± 0.01 g; Grobet®) and proceeded to count, weigh, and measure each of the neonates (using the same calliper and scale).

The field and experimental protocols were reviewed and approved by the Dirección General de Vida Silvestre (DGVS), Secretaría de Medio Ambiente y Recursos Naturales (Semarnat) through the scientific research permit FAUT-0074, granted to Fausto R. Méndez-de la Cruz. After completion of data collection in laboratory, we released all the females, together with their offspring, at the exact sites of capture.

To calculate RLM, we divided litter mass (excluding amniotic fluids and extraembryonic membranes) by post- parturition female mass (PPFM), as other traditional methods of calculation that include the female mass in the numerator and denominator, such as that proposed by ^{Tinkle (1972}; clutch or litter mass/mass of the gravid or pregnant female), are prone to statistical bias (^{Bastiaans et al., 2013}; ^{Rodríguez-Romero et al., 2005}; ^{Shine, 1980}). Previous studies have used body condition estimates to assess the effect of the energy status of pregnant females on their offspring (^{Recknagel & Elmer, 2019}). Here we calculated the body condition index of non-pregnant females, post-parturient females, and neonates, for which we used the scaled mass index proposed by ^{Peig and Green (2009)}, according to the following equation:

where M_i and L_i were defined as body mass (g) and SVL (mm) of the individual i, respectively; L₀ was the mean value of SVL of the species (using mean adult SVL for non-pregnant and post-parturient females and mean neonate SVL for the offspring; data from this study); and b_SMA, the scaling exponent obtained through dividing the slope of an ordinary least squares linear regression of M on L by the Pearson’s correlation coefficient. This index is size-independent and has been shown to be more suitable than other body condition indexes based on ratios or on residuals from regressions (^{Labocha et al., 2014}). We assumed that lower values of body condition index (BCI) in post-parturient females or offspring would indicate a detrimental effect of pregnancy or a reduced capacity to allocate enough resources to developing embryos, respectively (^{Bonnet et al., 2002}; ^{Recknagel & Elmer, 2019}).

Before any statistical procedure, we tested each variable for normality and homoscedasticity through Kolmogorov- Smirnov and Levene’s tests, respectively, considering an α level = 0.05 (^{Zar, 2010}). When data did not satisfy the statistical assumptions for parametric tests, we performed the non-parametric equivalent analyses (^{Dytham, 2011}). In the case of P. dugesii, data came from 2 different populations and sampling years, so we first compared them through Student’s t tests and Mann-Whitney U tests to discard intraspecific or temporal effects on the female and litter characteristics. As sexual dimorphism in body size might reveal the existence of fecundity selection (^{Cox et al., 2007}), we compared SVL between males and females through Student’s t tests and Mann-Whitney U tests. Also, we compared the BCI of post-parturient versus non-pregnant females of each species through Mann- Whitney U tests.

We analysed the litter and neonate characteristics through a series of general linear models. First, to compare RLM and litter size between the 2 species, we performed analyses of covariance (ANCOVA) with female SVL, AGL, and BCI as covariables, discarding non-significant factors in a stepwise procedure (^{Bastiaans et al., 2013}). Secondly, we performed Pearson’s linear correlations to assess the relation between litter characteristics (litter size, litter mass, and RLM) with the female characteristics (SVL, AGL, PPFM, and BCI). Finally, to determine the existence of trade-offs in offspring size versus number, we assessed the relation between litter size and mean SVL, mean mass, and mean BCI of neonates with female characteristics (SVL, AGL, PPFM, and BCI) as covariables, through an ANCOVA (^{Bastiaans et al., 2013}). We carried out the statistical tests in SPSS 15.0.1 (SPSS Inc., 2006) and SigmaPlot version 11.0 (Systat Software, 2008).

Results

We collected a total of 25 individuals of Plestiodon sp. (11 males, 7 non-pregnant females, 7 pregnant females) and 54 of P. dugesii, of which 33 were collected in 2013 (14 males, 7 non-pregnant females, 12 pregnant females) and 21 in 2014 (5 males, 7 non-pregnant females, 9 pregnant females). In laboratory, parturition of Plestiodon sp. occurred from mid-April through the end of May; in the case of P. dugesii, births occurred only during May. Minimum size (SVL) of pregnant females was 55.01 mm in Plestiodon sp. and 55.79 mm in P. dugesii, while the size of non-pregnant females ranged from 44.04 to 63.39 mm in Plestiodon sp. and 42.04 to 53.61 mm in P. dugesii. As there were no statistical differences in SVL, litter size, RLM, and BCI of P. dugesii by population or year, we pooled together the data: SVL by population: (males: t = -0.501, df = 13, p = 0.625; females U = 48.000, p = 0.448); litter size by population: (U = 24, p = 0.100); RLM by population: (U = 39, p = 0.669); RLM by year: (Mazamitla: t = 1.258, df = 4, p = 0.277; Tapalpa: t = 1.171, df = 13, p = 0.263); BCI of post-parturient females by population: (t = -0.009, df = 19, p = 0.993). Reproductive characteristics of the 2 species are summarized in Table 1.

Males of Plestiodon sp. (SVL x = 61.486 mm, SD = 4.462 mm, n = 11) were larger than males of P. dugesii (SVL x = 52.983 mm, SD = 5.503 mm, n = 19) (t = 4.353, df = 28, p < 0.001); however, there were no differences in SVL between females of the 2 species (Plestiodon sp.: x = 57.151 mm, SD = 7.660 mm, n = 14; P. dugesii: x = 56.536 mm, SD = 15.282 mm, n = 35; U = 210, p = 0.445). We found no evidence of sexual size dimorphism in either of the 2 species (Plestiodon sp.: t = 1.664, df = 24, p = 0.110; P. dugesii: t = 0.333, df = 53, p = 0.097). Post-parturient females had a higher BCI than non-pregnant females in both species (Plestiodon sp.: U = 7, p = 0.026; P. dugesii: U = 1, p < 0.001; Fig. 1).

Figure 1 Scaled mass body condition index of non-pregnant females (NPF) and post-parturient females (PPF) of Plestiodon sp. and P. dugesii. Boundaries of the boxes indicate the 25^th and 75^th percentiles; solid line, the median, and whiskers, the 10^th and 90^th percentiles. 

Regarding the comparison of litter size, there were no differences by species (F = 2.529, df = 1, p = 0.124) after discarding SVL, AGL, and female BCI, whose effects were non-significant (SVL: F = 0.088, df = 1, p = 0.769; AGL: F = 1.414, df = 1, p = 0.247; female BCI: F = 3.135, df = 1, p = 0.090). In the case of RLM, the effects of female SVL and AGL were also not significant (SVL: F = 0.755, df = 1, p = 0.394; AGL: F = 0.158, df = 1, p = 0.695), so we discarded them. However, in the final model female BCI had a significant effect on RLM (F = 5.133, df = 1, p = 0.032), with higher levels of BCI resulting in lower RLM, although without statistical differences between species (F = 0.223, df = 1, p = 0.641) (Table 2).

Female SVL has a marginally significant positive effect on neonate mean SVL (F = 4.273, df = 1, p = 0.053) but not female AGL, PPFM, and female BCI; after correcting for their effects, there were interspecific differences in neonate mean SVL, indicating a larger SVL in neonates of Plestiodon sp. than in neonates of P. dugesii (F = 8.630, df = 1, p = 0.009). There were no significant differences in neonate mean mass or BCI by species after accounting for female morphometric measurements (Table 3). Similarly, there were no differences in mean neonate SVL (F = 0.023, df = 1, p = 0.543), mass (F = 0.005, df = 1, p = 0.187), or BCI (F = 1.987, df = 1, p = 0.186) by population in P. dugesii. For Plestiodon sp. there were no significant relationships between any of the reproductive and female characteristics; for P. dugesii, on the other hand, we found negative relationships between litter size and neonate mean mass (r = -0.548, t = -2.541, p = 0.023) and neonate mean BCI (r = 0.549, t = -2.544, p = 0.022; Fig. 2).

Figure 2 Relationship between neonate mean body condition index and litter size in P. dugesii. Solid line depicts the regression line; dashed lines, the 95% interval of confidence, and dotted lines, the prediction interval. 

Discussion

In lizards, sexual size dimorphism has evolved in cases in which: 1) larger size in males increases mating success through an advantage in intrasexual competition, 2) larger size in females determines a higher number of offspring (i.e., clutch or litter size), or 3) differences in body size allow resource partitioning to reduce intersexual competition (^{Cox et al., 2007}). Some oviparous species of the genus Plestiodon, including P. elegans and P. laticeps, exhibit male-biased sexual dimorphism (^{Vitt & Cooper Jr., 1985}; ^{Wei-Guo & Xiang, 2001}). Lack of sexual size dimorphism (i.e., monomorphism) in Plestiodon sp. And P. dugesii might indicate that neither sexual nor fecundity selection or differential use of resources have driven the evolution of size in these species. However, it is known that sexual size monomorphism can occur because of the mutual cancellation of sexual and fecundity selection, a pattern that has been found in other scincomorph lizards (e.g., ^{Dubey et al., 2011}; ^{Ma et al., 2018}). Whatever the cause of this phenomenon, it is likely to be shared by other members of the Plestiodon brevirostris group. This is also consistent with the fact that litter size and mass were not related to female SVL in these species.

Species of the P. brevirostris group exhibit variable patterns of maternal SVL in relation to litter characteristics. For instance, in another population of Plestiodon sp. from Puebla litter size seems not to be influenced by female SVL (^{Chávez-López, 2012}), but in P. copei there is a positive relationship between these traits (^{Guillette, 1983}). Furthermore, ^{Chávez-López (2012)} found, for the referred population of Plestiodon sp. from Puebla, a positive relationship between litter mass and female SVL, a pattern that we did not find for the same species. However, this may be because of our small sample size for Plestiodon sp. Mean values of RLM in Plestiodon sp. and P. dugesii are consistent with those reported for other species of skinks (^{Qualls & Shine, 1997}; ^{Shine, 1980}; ^{Vitt & Price, 1982}), although in some females RLM accounted for about 50% of the PPFM, a high level of reproductive effort considering the active foraging strategy and elongated body plan of these lizards (^{Griffith, 1994}; ^{Vitt & Price, 1982}).

Variation in RLM in both species was influenced by BCI but not by other female characteristics, such as SVL and AGL, which supports its usefulness as an indicator of the energy status of the organisms in a reproductive context. Although in many squamates female body condition tends to decrease because of energy expenditure during reproduction (e.g., ^{Bonnet et al., 2002}), there are cases in which reproductive investment does not impair maternal body condition (^{Itonaga et al., 2012}). In Niveoscincus ocellatus, another viviparous skink, abdominal fat body mass increases with pregnancy (^{Wapstra & Swain, 2001}). However, while N. ocellatus is pregnant in summer, when food resources are more abundant, members of the Plestiodon brevirostris group have a fall reproductive cycle, and therefore gestation occurs during a period of limited food availability (^{Chávez-López, 2012}; ^{Guillette, 1983}). Furthermore, the narrowing of the activity periods in winter due to lower temperatures may cause constraints both in foraging and digestion efficiency of females (^{Hare & Cree, 2016}; ^{Méndez-de la Cruz et al., 2015}). Given the above, it seems unlikely that the higher BCI of post-parturient Plestiodon sp. and P. dugesii is due to an ability to store fat during pregnancy (which would imply increasing either the quantity or quality of the food consumed). Another, more reasonable, explanation would be that the onset of reproduction in females of these species is conditioned to some degree by their body condition, for example, to mobilize the required maternal reserves to sustain vitellogenesis (^{Ramírez-Pinilla et al., 2015}). Nonetheless, we cannot rule out the possibility that the higher BCI of postparturient females is an artifact of the ad libitum feeding in captivity, so these results must be taken cautiously.

Numerous studies have addressed the existence of trade-offs between clutch/litter size and offspring size or mass (e.g., ^{Bastiaans et al., 2013}; ^{Olsson & Shine, 1997}; ^{Recknagel & Elmer, 2019}). The negative relationship between litter size and mean neonate mass and BCI in P. dugesii, although moderate, reveals that a larger litter size may reduce in some degree the probability of offspring survival due to an impairment in body condition (^{Uller & Olsson, 2010}). In Plestiodon sp. on the contrary, a larger litter size does not imply a reduction in neonate mass or BCI, and neonates of this species had a larger SVL than in P. dugesii. The above suggests differential pressures on offspring survival between the 2 species, which might be related to the environmental conditions of their habitat (^{Díaz et al., 2005}). For instance, while in both localities of P. dugesii habitat provides plenty of shelter places and a reasonably high offer of food resources, the habitat at the locality of Plestiodon sp. is more degraded, and therefore offspring might be more prone to predation and less able to find prey.

^{Meiri et al. (2011)} proposed to analyse several components of the life histories of lizards in terms of their reproductive output and biomass production (i.e., productivity). They defined “productivity” as the total mass of offspring produced in a year (hatchling or neonate mass x clutch or litter size x number of clutches or litters per year), and “specific productivity” as productivity divided by female mass. Since it is to be expected Plestiodon sp. and P. dugesii to produce 1 litter per year, as occurs with other members of the brevirostris group (^{Chávez-López, 2012}; ^{Guillette, 1983}), productivity of Plestiodon sp. would be 0.85 g/year with a specific productivity of 0.21. In P. dugesii productivity would be 1.001 g/year and specific productivity 0.28. These values are consistent with those reported by Meiri et al. (2011) for other viviparous lizards from temperate environments. Furthermore, the authors found that both measures of biomass production tend to decrease in colder environments because mean annual temperatures constrain the duration of the reproductive season. The latter point is also consistent with the lower values estimated for Plestiodon sp., as this species lives at higher altitude and in a more thermally challenging environment.

Coupled with environmental factors, evolutionary history of lizards might determine patterns of reproductive investment (^{Uller & While, 2015}). Plestiodon sp. and P. dugesii diverged about 10 Ma ago (^{Bryson et al., 2017}), possibly a short time to produce substantial differences in their reproductive ecology, but enough to influence neonate characteristics. It remains to be assessed the extent to which certain evolutionary or environmental factors are responsible for shaping these patterns in species of the Plestiodon brevirostris group. Another aspect that deserves to be addressed is the evolution of parity mode in Plestiodon. Currently, it is considered that viviparity arose only once within the genus, in the P. brevirostris group (^{Watson et al., 2014}); however, recent phylogenies have placed P. sumichrasti, a Central American, lowland oviparous species (^{Miller, 1997}), nested within this clade (^{Brandley et al., 2012}; ^{Bryson et al., 2017}). Whether this represents 2 distinct origins of viviparity in the genus or an evolutionary regression to oviparity, in P. sumichrasti has not yet been tackled.

In summary, our results indicate that in spite of the overall similarities in the reproductive characteristics of Plestiodon sp. and P. dugesii, differences in their reproductive investment, and particularly in the relationship between offspring size and number, exist between these 2 species. Further studies are needed to deepen into the reproductive biology of other viviparous members of the genus Plestiodon, which in comparison to their oviparous congeners have received less attention.