Introduction

In environments characterized by abiotic and biotic heterogeneity, the interaction between environment and genotype may result in life history divergences among different populations of the same species due to phenotypic plasticity and local adaptation (Endler 1977). Environmentally-induced divergent selection may drive the evolution of traits that increase the relative fitness of a population in its local habitat more than other habitats (Williams 1966; Kawecki and Ebert 2004). Teasing apart intrinsic and extrinsic sources of variation in life history traits is the first step to understanding local adaptation, which is not only important for testing hypotheses of adaptations shaped by specific environmental factors, but also for identifying the selective forces at work (Williams 1993; Kawecki and Ebert 2004).

In oviparous species, embryos develop outside the mother's body and are sensitive to environmental change within the nest. Nest environments, including physical factors and predation pressure, can vary considerably among populations of the same species (Doody et al. 2006; McKinnon et al. 2010). This kind of variation may result in the evolution of adaptive embryonic development with respect to their local nest environments. Identifying and investigating local adaptation of embryonic development among populations is an important part of our understanding of the ecological and evolutionary processes operating on organisms (e.g. population dynamics, and life-history evolution); the developmental environment not only determines the survival of embryos, but also has long-term impacts on offspring quality in oviparous species (Deeming 2004; Noble et al. 2018). However, traditional studies on local adaptation mainly focus on post-embryonic stages of development, rather than the embryonic stage (Sears and Angilletta 2003; Kawecki and Ebert 2004; Sanford and Kelly 2011). Recently, however, an increasing number of studies have demonstrated local adaptation in embryonic development among geographically separated populations. For example, amphibian and reptile embryos develop more rapidly in high-latitude populations, compensating for the slower development induced by low temperatures in the field (Laugen et al. 2003; Du et al. 2010; Zhao et al. 2015). Another study showed that turtle embryos are capable of fine-scale thermal adaptation; embryos from a hot nest environment have a higher thermal tolerance than those from a cold nest environment (Weber et al. 2012). These studies were conducted in common-garden laboratory experiments, designed to elucidate how temperature differentials within incubation environments can shape the local adaptation of embryonic developmental rate and thermal tolerance among populations. However, these kinds of laboratory-based experiments do not fully reflect the range of complex environmental variables (e.g. temperature, moisture, oxygen and their interactions) in natural nests that may be responsible for driving the adaptive evolution of embryonic development. Thus, reciprocal egg-swap experiments are essential and offer a more robust way to elucidate local adaptation of embryonic development in response to these diverse variations of nest environments in different populations of oviparous species. However, such studies are currently still rare, partly because nesting environments are very difficult to find in the field for many species (e.g., squamates) (Li et al. 2018a; Tiatragul et al. 2019).

The toad-headed agama (Phrynocephalus przewalskii) is a short-lived desert lizard that is widespread in northern China. This species possesses male-mediated gene flow and has two major lineages separated by the Helan–Yin mountain chain located in the central part of its distributional range (Urquhart et al. 2009). Local adaption is thus expected to occur in populations from these two lineages because the high mountain range may restrict the passive dispersal of lizards and, therefore, gene flow between populations, which is a prerequisite for local adaptation (Kawecki and Ebert 2004).

Here we carried out reciprocal egg-swap experiments to identify the intrinsic (population origin) and extrinsic (nest environment) sources of life-history traits related to developmental success (incubation period, hatching success, and offspring growth and survival) in two populations located on the north and south sides of the Yin mountain chain. We envisage that these embryonic and offspring traits may respond adaptively to specific environments in one of the following four scenarios. Firstly, these traits show local adaptation in all populations: a set of specialist phenotypes each maximizing fitness in one habitat (Kawecki and Ebert 2004). In this case, the fitness of embryos and offspring from either population is greater in the local nest environment than it would be for embryos and offspring from other populations (Fig. 1a). Secondly, these traits show local adaptation only in one population, i.e. there is a greater fitness benefit for embryos and offspring associated with their local habitat for one population, but not for the other population (Fig. 1b). Thirdly, these traits are greater in one population than the other one regardless of nest environment, because a single phenotype, superior in all habitats, has evolved in this population, or females in one population have evolved to allocate greater resources per egg (thereby higher quality offspring by maternal effects) than the other population (Fig. 1c). Lastly, these traits do not differ among populations from either site, because a single generalist phenotype showing similar adaptation to all habitats has evolved in all populations, or plasticity enables the two populations to match the expressed phenotypes with the local conditions (Fig. 1d).

Fig. 1
figure 1

Hypotheses of adaptive evolution of embryos from different populations (Population A and Population B) in response to local nest environments. a Local adaptation in all populations, b local adaptation in one population, c a single phenotype superior in all habitats (“super phenotype”), and d a single generalist phenotype evolved in all populations. Lines with blue and orange colors indicate different populations

Full size image

Methods

Study species and sites

The toad-headed agama (P. przewalskii), a small oviparous lizard (to 55 mm snout–vent length), is widely distributed in the desert and semi-desert areas of northern China (Zhao et al. 1999). In the reproductive season (May–July), females select open sites with little vegetation to construct burrow nests (Li et al. 2018a). The nest consists of an inclined tunnel and an expanded chamber at the end where females deposit their eggs (Li et al. 2017). Mean clutch size varies among populations, ranging from 2.4 to 3.9 (Qu et al. 2011; Zeng et al. 2013; Sun et al. 2018). Offspring emerge from nests in August and reach maturity in the next summer (Li et al. 2018a).

Our study sites were located at the Gegentala of Siziwangqi (41°79′ N, 111°80′ E; elevation 1390 m) and Shierliancheng of Jungarqi (40°12′ N, 111°07′ E; elevation 1015 m) of Inner Mongolia, 250 km apart (Fig. S1a). Gegentala is located in the north of the Yin mountain range where the habitat is characterized by desert steppes dominated by herbaceous flowering plants of the genus Chenopodium (Fig. S1b). Shierliancheng is located in the south of the Yin mountains and the habitat is predominantly composed of sandy patches and shrub clusters dominated by Artemisia ordosica (Fig. S1c).

Despite significant between-site differences in lizard morphology (Fig. S1d, e), a phylogenetic analysis based on the sequences of ND2 and ND4 genes indicated that the lizards from the northern site (Gegentala) and the southern site (Shierliancheng) belong to the same species of P. przewalskii (Fig. S2).

Maternally selected nest environments

In June 2017, we captured gravid females with large eggs close to oviposition. Each female was given a paint mark on their dorsal surface using marker pens (Deli permanent marker 6824; Deli, China), and immediately released at the location from which they were collected. In the days following we searched for marked females every two hours throughout the day until they were found constructing nests (Li et al. 2018a). In total, we located 9 nests in the northern site (Gegentala) and 14 nests in the southern site (Shierliancheng) (Fig. S3). Once a nest was located, we excavated the nest to determine its internal structure by measuring the depth, width, height, and length of the burrow (Li et al. 2017). These data provided basic information for subsequent artificial nest construction at the northern and southern sites.

Reciprocal egg-swap experiments

In June, 2018, we collected 40 gravid females from the northern and southern populations, respectively. Females were taken back to the laboratory at our field study sites and housed individually in terraria (310 × 210 × 180 mm), with a 30-mm layer of moist sand as the substrate. A 100-W incandescent bulb was suspended over each terrarium to provide thermoregulatory opportunities for the lizards from 08:00 to 17:00 h. Food (mealworms and crickets), and water were provided ad libitum. Each terrarium was checked at least three times a day for freshly laid eggs. Once eggs were found, clutch sizes were recorded and eggs were individually weighed (± 0.001 g). Following oviposition, females were released at the sites from which they had been collected.

It is very hard, if not impossible, to locate enough maternal nests of P. przewalskii for a statistically robust sample size within a reproductive season, so we used artificially constructed nests rather than maternally selected nests for this study. To create reciprocal egg-swap experiments (two populations × two nest environments), we first constructed artificial nests in the natural habitat of lizards in the two field sites respectively, following the methodology of Li et al. (2017). Briefly, we randomly selected open sites with sparse vegetation coverage to construct artificial nests. The artificial nests mimicked the mean depths of the natural nests in the northern (ca. 17 cm) and southern (ca. 15 cm) sites respectively. Then we paired eggs from one northern female and one southern female, with similar oviposition dates, and assigned one egg from each clutch of each population to the artificial northern and southern nests, respectively.

Abiotic characteristics of artificial nests

We measured several characteristics of the physical environment of artificial nests. We quantified vegetation cover (%) within a 1 × 1-m2 quadrat with the burrow entrance as the center and placed an iButton temperature logger (DS1921; Maxim Integrated Products Ltd., USA) at the bottom of each nest to record hourly temperatures in all paired nests (n = 40). In addition, we measured nest moisture and oxygen levels once a week during the incubation period in 20 paired nests. We calculated nest moisture as the water content (%H2O) of underground soil samples (taken from the equivalent depth of the natural nest), just beside each natural nest (within 10 cm). The soil moisture was calculated using the following formula: (wet mass-dry mass)/wet mass × 100%. Nest oxygen content was measured using a portable digital oxygen analyzer (CY-12C, Meicheng Electrochemical Analytical Instruments Factory, Hangzhou, China) (Li et al. 2017).

We allowed eggs to develop in the artificial field nests for 31 days, after which we excavated the nests and recorded egg survival. Viable eggs were transferred to the laboratory at our field sites where they were incubated individually in plastic cups (200 mL) containing a 5-cm-thick sand substrate with a water content of 4%. We sealed the cups with plastic wrap to minimize evaporative water loss and maintained them at a constant temperature of 26 °C in an incubator (PT2499; Exo-Terra, Canada).

We checked the eggs in cups twice a day for hatching. Newly emerged hatchlings were then measured (snout–vent length [SVL] ± 0.01 mm) and weighed (body mass [BM] ± 0.001 g), and individually toe-clipped for identification. A total of 87 hatchlings (49 from the northern population and 38 from the southern population) successfully hatched from 40 artificial nests in each site. These hatchlings were then transferred to our laboratory in Beijing and randomly assigned into five terraria (100 × 40 × 37 cm, L × W × H). The hatchlings from the two populations were reared together, with 17–18 individuals in each terrarium. Food (mealworms) and water were provided every day ad libitum. We provided hatchlings with 10 h of heating for behavioral thermoregulation (from 08:00 to 18:00 h, simulating the summer photoperiod at the collection sites) by placing a 100-W incandescent bulb above one side of each terrarium. We measured locomotor performance of the young lizards 10 days after hatching in a temperature-controlled room (36 ± 2 °C). Prior to testing, we acclimated the hatchlings in the room for approximately 1 h. We encouraged hatchlings to run along a 1-m racetrack by touching their tails slightly with a soft paintbrush. Locomotor trials were recorded using a Sony DCR-SR220E digital video camera (Sony, Japan) mounted above the racetrack. We tested the locomotor performance of each hatchling twice with an interval of one hour between trials. The locomotor trial videos were analyzed later using AVS Video Editor software (Online Media Technologies Ltd, UK) to measure sprint speed during the fastest 20 cm interval, and average speed over the length of the 1-m racetrack. On 22 August 2018, we re-measured the hatchlings (SVL and BM) to calculate post-hatching growth rates. In addition, we monitored hatchling survival by recording the date of death for hatchlings until the end of October 2018.

Statistical analyses

We used Shapiro–Wilk normality test and Levene’s tests to check the data for normality and homogeneity of variance. Data were normalized by exponential or log-transformation when necessary. We used t tests to compare between-population differences in the microhabitat variables of natural and artificial nests, female body sizes and reproductive traits. We added correlated factors (tested by Pearson correlation analysis) as covariates into linear mixed models and generalized linear mixed models, and then compared these models using the likelihood ratio test to choose the best fitted model for subsequent data analysis. The between-population difference in egg mass was evaluated using a linear mixed model, with population as a fixed effect, maternal identity as a random effect and maternal snout–vent length as a covariate. We used generalized linear mixed models to detect the effects of population and nest environment and their interaction on the incubation period and hatching success of eggs using a Poisson distribution and a binomial distribution, respectively, with maternal identity as a random factor, and egg mass as a covariate if necessary. The effects of population and nest environment and their interaction on hatchling morphology were tested by linear mixed models, with maternal identity as a random factor and egg mass as a covariate. We used linear mixed models to detect the effects of population and nest environment and their interaction on growth rates and locomotor performance, with maternal identity as a random factor, and Cox regression to detect the effects of population and nest environment on hatchling survival.

Results

In their natural habitat, females from the northern population constructed deeper nests with a longer and overall larger burrow than those from the southern population (Table S1). Additionally, females from the northern population were larger in body size, and produced larger clutches with bigger eggs than females from the southern population (Table S2), but the between-population difference in egg mass disappeared after we had removed the effect of maternal body size statistically (F1,77 = 1.057, P = 0.307). The artificial nests we constructed successfully mimicked maternally constructed nests in the field; the northern nests were deeper, had higher vegetation coverage and moisture levels, but lower temperatures (mean, daily maximum and minimum temperatures) than the artificial southern nests (Table 1).

Table 1 Microhabitat variables of artificial Phrynocephalus przewalskii nests in the northern (Gegentala) and the southern (Shierliancheng) sites of Inner Mongolia, China
Full size table

Hatching success was slightly affected by egg mass, with heavier eggs having higher hatching success (z = 1.856, df = 154, P = 0.064). Hatching success was higher in the northern population than the southern population (z = − 1.997, df = 155, P = 0.046), but did not differ between populations after we had removed the effect of egg mass statistically (z = − 1.644, df = 154, P = 0.100). Before and after we had removed the effect of egg mass, hatching success was higher at the northern nest environment than the southern nest environment (before: z = − 1.982, df = 155, P = 0.048; after: z = − 1.967, df = 154, P = 0.049), and was not affected by the interaction between population and nest environment (before: z = 0.960, df = 155, P = 0.337; after: z = 0.952, df = 154, P = 0.341). Further analyses indicated that hatching success of the northern population was higher than that of the southern population in the northern nest environment (z = − 2.042, df = 77, P = 0.041), but not at the southern nest environment (z = − 0.894, df = 77, P = 0.371). The effect of egg mass on hatching success was significant at the northern nest environment (z = 2.199, df = 76, P = 0.028), but not at the southern nest environment (z = 0.737, df = 76, P = 0.461). After we had removed the effect of egg mass statistically, the between-population difference in hatching success at the northern nest environment disappeared (z = − 1.431, df = 76, P = 0.153) (Fig. 2a).

Fig. 2
figure 2

Hatching success and incubation period of eggs from the northern (Gegentala) and southern (Shierliancheng) populations of Phrynocephalus przewalskii from Inner Mongolia, China, when incubated at the northern and southern nest environments, respectively. Data are expressed as means ± SE

Full size image

Incubation period did not differ between populations (z  = 1.497, df = 82, P = 0.134), but was significantly shorter for eggs incubated at the southern nest environment than those from the northern nest incubation environment (z = − 4.410, df = 82, P < 0.001). However, there was no interaction between population and nest environment on incubation period (z  = − 0.236, df = 82, P = 0.813; Fig. 2b).

Hatchlings from the southern population were larger but not heavier and ran faster than those from the northern population. Hatchlings from eggs incubated at the northern site were larger and heavier than those from the southern site, but the locomotor performance of hatchlings was independent of nest sites. In addition, hatchling morphology and running speeds were not affected by the interaction between population and nest environment (Table S3, Fig. 3).

Fig. 3
figure 3

Hatchling body size and locomotor performance of Phrynocephalus przewalskii from the northern (Gegentala) and southern (Shierliancheng) populations of Inner Mongolia, China, when incubated at the northern and southern nest environments, respectively. Data are expressed as means ± SE

Full size image

Hatchling growth rates in SVL and body mass were significantly affected by the interaction between population and nest environment (SVL: F1,44.05 = 7.740, P = 0.008; body mass: F1,49.28 = 6.248, P = 0.016). The southern hatchlings grew faster than did the northern hatchlings when incubated at the southern nest environment (SVL: F1,30 = 21.246, P < 0.001; body mass: F1,30 = 9.631, P = 0.004), but not at the northern nest environment (SVL: F1,45 = 1.891, P = 0.176; body mass: F1,45 = 0.102, P = 0.751; Fig. 4 a,b).

Fig. 4
figure 4

The post-hatching growth of hatchling Phrynocephalus przewalskii from the northern (Gegentala) and southern (Shierliancheng) populations of Inner Mongolia, China, when incubated at the northern and southern nest environments, respectively. Data are expressed as means ± SE

Full size image

Egg mass did not affect the survival rate of hatchlings (χ2 = 1.313, P = 0.252). Survival rate of southern hatchlings was significantly higher than that of the northern hatchlings (χ2 = 4.952, P = 0.026; Fig. 5); however, overall, offspring survival was not affected by nest environment (χ2 = 0.256, P = 0.613) or its interaction with population (χ2 = 0.369, P = 0.543; Fig. 5).

Fig. 5
figure 5

The survival of hatchling Phrynocephalus przewalskii from the northern (Gegentala) and southern (Shierliancheng) populations of Inner Mongolia, China, when incubated at the northern and southern nest environments, respectively. ‘nn’ and ‘ns’ indicate the northern population incubated at the northern nest environment (nn) and the southern nest environment (ns), respectively. Similarly, ‘sn’ and ‘ss’ indicate the southern population incubated at the northern nest environment (sn) and the southern nest environment (ss), respectively

Full size image

Discussion

Our study demonstrates that lizard embryos show local adaptation to their nest environments among populations, adding to the results from previous studies that showed behavioral and physiological responses of embryos to environmental variations within nests (see reveiws by Du and Shine 2015; Du et al. 2019). Below we discuss the adaptive strategies driving embryonic development in response to local environments.

Maternal nest selection not only maximizes the developmental success of embryos but can also optimize offspring phenotypes to enhance fitness (Refsnider and Janzen 2010; Mitchell et al. 2013; Reedy et al. 2013). Consequently, both embryo and offspring phenotypes are needed to evaluate the adaptation of embryos to their environments. Consistent with the prediction of local adaption in one population (Fig. 1b), embryos from the northern population are better adapted to their local nest environments than embryos from the southern population, in terms of embryonic survivorship during development. However, the opposite pattern is true in terms of offspring growth (Figs. 2, 4). This seemingly contradictory phenomenon may reflect different strategies for local adaptation developed by the two populations in response to their own environments, which may be due to the existence of divergent selective pressures within different life stages among populations (Salvanes et al. 2004). Embryos from the northern population achieve higher fitness via an increase in embryonic survivorship, whereas embryos from the southern population gain a higher fitness potential via an increase in offspring growth rates. This enables offspring to mature and reproduce earlier, increasing lifetime fitness in terms of fecundity (Stearns and Koella 1986; Li et al. 2018b). Given the male-mediated gene flow that exists between populations in this species (Urquhart et al. 2009), we suspect that local adaptation of embryonic development and offspring growth rates may be achieved via genetic divergence in maternal lineages. Previous studies have shown that embryonic development and growth rates are closely related to mitochondrial functions (Sun et al. 2015), and could thus be maternally inherited through mitochondria. Molecular studies to verify this hypothesis would be of great interest in future research.

The higher survivorship of northern population embryos incubated in their local nest environment indicates that the northern embryos have developed specific strategies in adaptation to the cold and wet nest environments prevalent in the northern locality. The larger egg size in the northern population is likely an adaptive trait that may facilitate successful development in the northern nest environment, because the between-population difference in hatching success disappeared after the effect of egg size has been removed statistically. Differences in nest temperatures and moisture levels between the northern and southern nest environments may also account for the between-population difference in embryonic survivorship, because temperature and moisture can significantly affect hatching success in reptiles (Deeming 2004; Du and Shine 2015). Nonetheless, our experimental design does not allow us to identify the specific environmental force driving this adaptation. Future common-garden experiments are needed to tease apart the role of individual environmental factors of nests on embryonic survivorship, and thereby identify the underlying selective forces (Kawecki and Ebert 2004). Such environmental factors not only include abiotic factors like temperature and moisture (Deeming 2004; Du and Shine 2015), but also involve biotic factors such as predation which contributes to egg mortality in reptiles and birds (Spencer and Thompson 2003; McKinnon et al. 2010).

It is noteworthy that our study cannot exclude the contribution of other factors (e.g., maternal effects, and genetic drift) in addition to local adaptation to the between-population difference in life-history traits. As seen in many other lizards (Angilletta et al. 2004; Du et al. 2005), egg size differed among populations in our case (Table S2). This represents a between-population difference in maternal investment that may subsequently affect embryonic development and offspring traits. Given that egg mass correlated with hatching success and incubation period, the higher hatchability and developmental rate in the northern population could not only be caused by genetic adaptation, but also by between-population difference in maternal investment. However, egg size manipulation by yolk removal in some other lizards have demonstrated that the between-population difference in incubation period is not due to egg size difference (Oufiero and Angilletta 2006; Goodman 2010). Therefore, it would be of great interest to identify the relative contribution of genetic adaptation and maternal effects in shaping the population difference in these fitness-related traits by egg size manipulation experiments in future. In addition, given the genetic difference between the two populations (Fig. S2), the between-population differences in embryonic development and offspring traits might be partly due to lineage-specific evolutionary history (Coyne and Orr 2004; Massetti et al. 2019). Future studies on multiple populations from each geographic region could help to discriminate local adaptation from other driving factors (e.g., genetic drift, founder effects, etc.) underlying the between-population differences. Despite these caveats, our study demonstrated adaptive responses of embryonic development to local environments among lizard populations via reciprocal egg-swap experiments in field nests, adding new knowledge to reptilian embryogenesis, in which previous studies focused on developmental plasticity (Deeming 2004; Noble et al. 2018).

More generally in oviparous species, embryos are relatively immobile, and thus experience strong and geographically divergent selection pressures within the nest. Compared with post-embryonic life-history stages that are affected by a diversity of environmental factors and ecological relationships (e.g., food availability, competition), embryos experience relatively simple environmental conditions that are easier to measure in the field and manipulate in the laboratory (e.g., temperature, moisture) (Du et al. 2010). Combining field investigation and manipulative laboratory experiments enables us to determine the process and underlying driving forces of life history adaption to local environments (Weber et al. 2012; Mitchell et al. 2015; Li et al. 2018a). Therefore, among-population differences in embryonic development provide excellent opportunities for studying the adaptive evolution of life history traits adjusted to local selective forces.