Population Dynamics of Wild Mongolian Gerbils: Quadratic Temperature Effects on Survival and Density-Dependent Effects on Recruitment

: It has been hypothesized that animal populations respond nonlinearly to the environment, and such responses are important to understand the effects of climate change population dynamics of small mammals in arid environments at northern latitudes. The aim of this study was to test the following hypotheses: (1) that small rodent populations increase as their semiarid habitat conditions improve from low to intermediate levels of temperature or precipitation, and decline beyond the optimum climate because of decreased survival, and (2) that increased population density would result in stronger negative effects on recruitment than on survival. A wild population of Mongolian gerbils ( Meriones unguiculatus ), a granivorous rodent distributed in Inner Mongolia, China, was live-trapped half-monthly between April and October from 2014 to 2017 and the effects of climate and density on their apparent survival probabilities and recruitment rates were estimated using mark-recapture methods. Increased temperatures initially had a positive effect on population growth rates, and then had negative effects on population growth rates primarily, which was mediated by quadratic effects on survival probabilities, further supporting the optimum habitat hypothesis. Moreover, the increases in temperature had a positive effect on the recruitment of gerbils, whereas population density had a more markedly negative effect on recruitment than on survival. The results of this study suggested that the density-dependent feedback to recruitment may be a primary regulatory mechanism of small mammal populations, and the complex responses of populations to temperature, which is a limiting ecological factor, may raise concerns for the fate of populations of small mammals at northern latitudes, in view of the predicted global climate change scenarios.


Introduction
Understanding the differential and integral effects of exogenous (e.g., weather, food, and predation) and endogenous (e.g., population density) factors on consumer dynamics in different ecosystems has been a longstanding interest in population ecology [1,2]. Arid and semiarid grassland or desert environments at northern latitudes are precipitationlimited ecosystems with a comparatively high evaporation mediated by temperature, where resources are highly variable in space and time [3]. Rodents in these regions are ecologically important consumers and are even considered keystone species in some ecosystems [4,5]. Their populations often show pronounced seasonal or annual variation in numbers and structure with a demographic and/or environmental stochasticity [6][7][8]. Demographic stochasticity is determined by stochastic events such as birth and death, immigration and emigration [9]. Birth and immigration are referred to as recruitment in empirical studies, generally. Likewise, death and emigration are referred to as apparent mortality, and the the optimum [41]. However, to our knowledge, limited information is available in regard to the effects of weather and density on the survival and recruitment of wild Mongolian gerbils in the semiarid steppe ecosystem. The analysis of the effects of changes in climate and population density on Mongolian gerbils' populations may contribute elucidating the population dynamics of small mammals inhabiting arid environments within the context of the future climate warming scenarios predicted at northern latitudes.
The present study focused on the tendency and magnitudes of the effects of temperature, precipitation, and population density on the survival and recruitment of Mongolian gerbils using mark-recapture methods [42,43]. Tests were conducted to evaluate the nonlinear effects of temperature and precipitation on the, survival probabilities, recruitment and population growth rates of gerbils, as predicted by the optimum habitat hypothesis. The prediction that increased population density has stronger negative effects on gerbil recruitment than on survival was also tested.

Study Area
Field experiments were conducted in Houhatai country (42 • 23.613 N, 116 • 06.524 E, altitude: 1300 m), in south-central Inner Mongolia, China. The study site is a typical steppe area with a semi-arid climate characterized by warm summers and cold winters. the average monthly temperature ranged from −13.9 • C to 21 • C and the mean annual total precipitation was about 388.6 mm, ranging from 291.9.0 to 454.8 mm during the 2014-2017 period ( Figure 1). Snow cover lasted from mid-or late October to early April the following year. All temperature and precipitation data were obtained from Zhenglan Qi Meteorological Station (ID. 54205), which is located about 25 km south of Houhatai country.

Live Trapping and Individual Identification
Mongolian gerbils were live-trapped approximately half-monthly or biweekly from 29 April to 24 October 2014; from 4 May to 19 October 2015; from 27 April to 20 October 2016; and from 12 May to 17 September 2017. During these periods each trapping session lasted 3-5 consecutive days. Mongolian gerbils remain active mainly under the snow cover during winter (W. Liu, personal observation), therefore, no animals were trapped from November to March to avoid trap mortality due to severe weather.
Gerbils were captured using wire mesh live traps (28 cm × 13 cm × 10 cm) baited with fresh peanuts and adopting a concentric circle trapping method [24]. To enhance the likelihood of successful trapping, a trap station was set in three or four concentric circles within a burrow system to cover most of the range. Each trapping week, a total of approximately 150-380 traps were placed on the trapping plot, with 24-28 traps per burrow system and the animals were captured during the diurnal or crepuscular activity periods. The traps were set at 05:00-06:00 in the morning from May to August and were checked every The trapping plot used in this study covered a 2 ha (100 m × 200 m) grassland dominated by Leymus chinensis, Artemisia sieversiana, Thalictrum petaloideum, Stellera chamaejasme, Klasea centauroides, and Aster altaicus, which provided food or cover for the gerbils. Other sympatric small mammals found in the area included the Daurian ground squirrel (Spermophilus dauricus), Daurian pika (Ochotona dauurica), and striped hamster (Cricetulus barabensis). In this region, the potential predators of Mongolian gerbils were the steppe polecat (Mustela eversmanii), corsac fox (Vulpes corsac), and some raptors, such as the common kestrel (Falco tinnunculus) and upland buzzard (Buteo hemilasius). No livestock grazed on the study site during the study period [43].

Live Trapping and Individual Identification
Mongolian gerbils were live-trapped approximately half-monthly or biweekly from 29 April to 24 October 2014; from 4 May to 19 October 2015; from 27 April to 20 October 2016; and from 12 May to 17 September 2017. During these periods each trapping session lasted 3-5 consecutive days. Mongolian gerbils remain active mainly under the snow cover during winter (W. Liu, personal observation), therefore, no animals were trapped from November to March to avoid trap mortality due to severe weather.
Gerbils were captured using wire mesh live traps (28 cm × 13 cm × 10 cm) baited with fresh peanuts and adopting a concentric circle trapping method [24]. To enhance the likelihood of successful trapping, a trap station was set in three or four concentric circles within a burrow system to cover most of the range. Each trapping week, a total of approximately 150-380 traps were placed on the trapping plot, with 24-28 traps per burrow system and the animals were captured during the diurnal or crepuscular activity periods. The traps were set at 05:00-06:00 in the morning from May to August and were checked every 1-2 h until about 11:00. They were then closed between 11:00 and 15:00 to avoid trap mortality due to the heat, and trapping was resumed at 16:00 and continued until 19:00. During September and October, trapping started between 06:30 and 07:30, and continued until 17:30. The traps were also checked every 1-2 h during this period [24,43].
All the captured gerbils were immediately toe-clipped to allow permanent identification, and they were sexed and weighed to the nearest 0.1 g using an electronic balance (JY/YP 6001, Yueping Scientific Instrument Co., Ltd., Shanghai, China). The ID number, trap location (i.e., burrow system ID), sex, body mass, and reproductive condition (e.g., reproductive males with scrotal testes and visible ventral scent glands or reproductive females with a bulging abdomen, enlarged nipples surrounded by white mammary tissue, or an opened pubic symphysis) were recorded [24]. The gerbils' age was estimated based on body mass and the developmental stage of the ventral sebaceous gland: juveniles were identified by a body mass < 30 g and no sign of ventral gland; subadults by a body mass of 30-50 g and a ventral active gland wider than 4.2 mm; and adults by a body mass > 50 g [24,36]. The captured gerbils were released at the same trap station immediately after handling or necessary behavioral tests. The above animal procedures and methods were approved by the Animal Care and Use Committee of the Institute of Zoology, Chinese Academy of Sciences (Ethical Inspection License No: IOZ13047).

Analysis of Survival and Recruitment
The effects of weather and density on half-monthly apparent survival and recruitment of Mongolian gerbils was evaluated using the Jolly-Seber (JS) model [44,45]. Before the survival probability, recapture probability, and recruitment rate in the JS model were calculated, the Cormack-Jolly-Seber (CJS) model was first adopted to determine the covariates influencing the gerbils' recapture probabilities in order to reduce the number of possible combinations of temperature, precipitation, and density [44][45][46]. The goodness of fit test for the CJS model is based on contingency tables that verify the assumptions of parameter homogeneity across individuals of the same group [47]. The assumptions of no behavioral response (i.e., not trap-shy or trap-happy) and no movement (i.e., no transient individuals) were tested for the CJS model using the 2CT + 3SR tests (male, χ 2 46 = 42.83, p = 0.61; female, χ 2 52 = 63.13, p = 0.14) in the RELEASE program [48]. The CJS model was first built including sex and time (i.e., trapping half-monthly) as categorical variables (hereafter, categoricalcovariate models) to test whether the survival and recapture probabilities of gerbils were sex-dependent and time-dependent. Sixteen categorical-covariate CJS models were built with all possible combinations of the effects of time and sex on the probabilities of survival and recapture (i.e., four different combinations for survival and four different combinations for recapture; see Supplementary Materials Figure S1).
Information-theoretic approaches were used to select the most parsimonious model to infer the effects of covariates [49]. The model with the lowest value of the Akaike information criterion corrected for small sample sizes (AICc) or highest AICc weight (w) was the most parsimonious, and the models with a ∆AICc value < 2 were considered as competing models, in the sense that they competed with the best model with substantial support provided by the statistical empirical data [49]. The variance inflation factors, e.g., median c-hat, were estimated using the most complex model (global model) with the time-sex interaction affecting both survival and recapture probabilities [50]. If the median c-hat value was greater than 1.0, this factor was used to scale the model deviances and compute the quasi-Akaike information criterion corrected (QAICc) for small sample sizes in the model selection [50]. The median c-hat was estimated at 0.928, so as to not scale the model deviances for the CJS model in our analysis.
If the effect of sex on survival or recapture probabilities was included in the most parsimonious model or competing categorical-covariate CJS models, it was also included in subsequent survival analyses. If survival probabilities varied over time in the most competing categorical-covariate model, the same combinations of covariates were used for survival probabilities (hereafter, survival submodels) and recapture probabilities (hereafter, recapture submodels) to model the recruitment rates (hereafter, recruitment submodels) in the Pradel-f JS models [50,51].
To assess the net effects of weather and increases in population density on demographic parameters (e.g., survival, recruitment, and population growth rate), firstly, temperature, precipitation, and population density were included as covariates of survival probabilities with the time or sex categorical covariates for recapture probabilities, as in the most parsimonious categorical-covariate CJS model (hereafter, environmental-covariate models). Subsequently, different JS models were built by including different combinations of weather variables and population density to model recruitment via Pradel-f models or to model the population growth rate (λ) via Pradel-λ JS models [50,51] in Mongolian gerbils as in the most parsimonious environmental-covariate CJS model. Additional survival submodels with weather variables were also built only within the Pradel-f models. If density was included in the recruitment submodel but not in the survival model of the most competing model, it was concluded that changes in population size exerted stronger effects on recruitment than on survival. The variable importance index was used to compare the strength of the effect of a covariate between survival and recruitment. Specifically, the variable importance index of a covariate for survival or recruitment was computed as the sum of the Akaike weights of the models, whose survival or recruitment submodel included the covariate [49]. The greater the index value, the more important the covariate.
Weather variables included the average daily temperature ( • C) and average daily precipitation (mm) of each trapping interval between two successive trapping sessions. Daily temperatures and precipitation were significantly correlated (Pearson's r = 0.46, p = 0.002, Supplementary Materials Figure S1). A quadratic term of temperature or precipitation was included in the models to test the nonlinear effects on survival, recruitment, and population growth rate (λ). The minimum number of alive animals (MNA) and MNA per hectare (MNA/ha) were used to index the population size of gerbils [42] and population density, respectively. Through the AICc and Akaike weight, the best and competing models were selected to infer the effects of temperature, precipitation, and density on recruitment or population growth rate. The CJS and JS models were built in MARK version 8.3, and all the covariates, including weather variables, were included in the design matrix in this software [50].
Additionally, the half-monthly proportions of reproductively active males (PRM) and females (PRF) were calculated by dividing the number of reproducing males or females by the total number of males or females captured in a trapping session, respectively. The mean half-monthly population growth rates (r) were obtained using the formula: r = (15/T) ln (MNA t /MNA t−1 ), where r is the mean half-monthly population growth rate; ln ( ) is the natural logarithm function; MNA t−1 and MNA t are the MNAs at time t−1 and t, respectively; and T is the number of days between t−1 and t. Generalized least square (GLS) models were used to test for quadratic relationships between population growth rates and weather with an autocorrelated error of order 1 [52].

General Population Demography
In

Effect of Climate and Population Density on Survival and Recruitment
In the best categorical-covariate CJS model, both survival probabilities and recapture probabilities were time-dependent and sex-independent (Akaike weight w = 0.998). The second-best model, in which survival probabilities were time-and sex-independent and recapture probabilities were time-independent and sex-dependent, had a ΔAICc of 12.65

Effect of Climate and Population Density on Survival and Recruitment
In the best categorical-covariate CJS model, both survival probabilities and recapture probabilities were time-dependent and sex-independent (Akaike weight w = 0.998). The second-best model, in which survival probabilities were time-and sex-independent and recapture probabilities were time-independent and sex-dependent, had a ∆AICc of 12.65 (Supplementary Materials Table S1 step 1). Then the time-dependent and sex-independent recapture probabilities in the Pradel-f JS models were used to estimate the survival probability submodels and recruitment rate submodels [ϕ(t)p(t)f (t), Akaike weight (w) = 1.000]. The second-best Pradel-f model [ϕ(t)p(t)f(g*t)] had a ∆AICc of 26.69 (Supplementary Materials Table S1, step 2). Half-monthly survival ranged from 0.22 to 1.00, and recruitment rates ranged from 0.00 to 1.38 (Figure 3). Time-dependent and sex-independent recapture probabilities were used in the environmental-covariate CJS models ( Table 1). The most parsimonious of these models had survival probabilities (φ) as a function of population density (i.e., MNA/ha; coefficient = 0.013, 95% confidence interval (CI) = 0.007 to 0.019), precipitation (coefficient = 0.038, 95% CI = −0.107 to 0.182), and the linear (coefficient = 0.113, 95% CI = 0.068 to 1.583) and quadratic (coefficient = −0.002, 95% CI = −0.0045 to −0.0001) terms of temperature (Akaike weight = 0.983). The second-best model had a ΔAICc of 8.96 (Table 1, step1). Therefore, four combinations of covariates were used for the survival submodel and the recruitment rate submodel of the Pradel-f JS models (Table 1, step2): (1) temperature and population density, (2) temperature and precipitation, (3) temperature, precipitation and population density, and (4) a quadratic temperature term as the only covariates, the nonlinear effects of temperature.

Model Likelihood
Par Deviance Step 1: Modelling survival probability, including environmental covariates with half-monthly temporal (t)-dependent Time-dependent and sex-independent recapture probabilities were used in the environmental-covariate CJS models ( Table 1). The most parsimonious of these models had survival probabilities (ϕ) as a function of population density (i.e., MNA/ha; coefficient = 0.013, 95% confidence interval (CI) = 0.007 to 0.019), precipitation (coefficient = 0.038, 95% CI = −0.107 to 0.182), and the linear (coefficient = 0.113, 95% CI = 0.068 to 1.583) and quadratic (coefficient = −0.002, 95% CI = −0.0045 to −0.0001) terms of temperature (Akaike weight = 0.983). The second-best model had a ∆AICc of 8.96 (Table 1, step1). Therefore, four combinations of covariates were used for the survival submodel and the recruitment rate submodel of the Pradel-f JS models (Table 1, step2): (1) temperature and population density, (2) temperature and precipitation, (3) temperature, precipitation and population density, and (4) a quadratic temperature term as the only covariates, the nonlinear effects of temperature. AICc and ∆AICc represent the Akaike information criterion corrected for small sample size and the difference in AICc between the given model and the most parsimonious model, respectively. The models highlighted in bold were selected in each step and were used as the starting point for subsequent steps. Par denotes the number of identifiable parameters.
The best environmental-covariate Pradel-f JS model (∆AICc < 2.0, Table 1, step 2) included linear and quadratic temperature terms in the survival probabilities submodel, and included density and linear temperature in the recruitment submodel. The quadratic effect of temperature explained 51.8% of the variability in survival probabilities with a coefficient of −0.003, and its 95% confidence interval (CI) excluding zero (−0.006 to −0.001), and population density was not included in the best competing survival submodel. In contrast, population density (coefficient = −0.010, 95% CI = −0.016 to −0.004) and only linear temperature (coefficient = 0.047, 95% CI = 0.028 to 0.066) were included in the recruitment submodel of the best competing JS models with a variable importance index of 0.518 (Table 1, step 2). The second-best model included temperature, precipitation, and population size, with ∆AICc and AICc weight (w) measuring 2.04 and 0.186, respectively (Table 1,

Discussion
The present study demonstrated that the population abundance indices and survival probabilities of Mongolian gerbils (M. unguiculatus) distributed in the semiarid steppe of Inner Mongolia varied from 2014 to 2017, (Figures 2a and 3). It was also shown that the breeding season of gerbils was relatively short, beginning in April, peaking between June and July, and terminating by the end of July or early August (Figure 2b), which is in line with the findings of our previous studies on the agro-pastoral ecotone in this region [24,35]. The peak population sizes of the gerbil population declined from 2014 to 2017 (Figure 2a). In particular, our analyses showed that the differential and integral effects of temperature and population density on vital rates were key processes determining the population fluctuations of gerbils. At the same time, this study also provides a nonlinear interpretation of the complex role of climate as a limiting factor, indicating that temperature affected the population growth rates of gerbils through its quadratic effects on survival probabilities, and positive effects on recruitment. (Tables 1 and 2, and Figure 4).
There are two plausible explanations for the observed effects of temperature on the Mongolian gerbil population examined. First, the rise in temperature may significantly disrupt the physiological functions of the animals. Our previous laboratory studies found that temperature extremes significantly reduced body weight, body temperature, metabolic rate, and cognitive abilities in Mongolian gerbils [53][54][55], as well as in voles and rats [14,56]. High temperatures may also increase the risk of disease transmission mediated by decreased immunity [57]. A recent study showed that the gut microbiota was connected with thermoregulation in gerbils in response to changes in temperature [58]. In addition, other studies indicated that the reproductive performance of male gerbils was significantly reduced in hot environments, with a 92% reduction in sperm density and the complete loss of sperm activity [53], which is in line with the reported effects of high temperature on the impairment of spermatogenesis in rat testes [59]. Moreover, the reproductive performance of female rodents (e.g., common voles, Microtus arvalis), indicated by milk production and pup growth, was also shown to decrease in hot environments [60]. In our group previous studies, it was found that, when comparing female gerbils with nearly identical mean litter sizes, those lactating in hot conditions (30 • C) exported 32.0 kJ day −1 less energy as milk at peak lactation than those lactating in cold (10 • C) or warm (21 • C) conditions, with no difference between the latter two groups [61]. At the same time, transferring lactating gerbils from warm to hot conditions resulted in reduced milk production and, on the 14th day of lactation, the litter masses in cold and hot conditions were 12.2 and 9.3 g lower than those in the warm condition, respectively [61]. Therefore, cold or hot temperature extremes may result in physiological susceptible responses, leading to a reduction in the survival rates of Mongolian gerbils.
The second plausible explanation concerns climate warming, which may cause the variation in local rodent populations by altering plant composition. Previous studies have shown that climate change induced the loss of perennial and Cyperaceae species over the years in a tallgrass prairie [62], and even decreased the productivity of plants in grasslands [63]. Changes in the characteristics of grassland vegetation communities could reduce the quantity and quality of the plant species preferred by specific rodents, potentially causing severe impacts on the survival of offspring, and on their body regulation ability [61]. A recent study reported that increasing temperatures may also alter the gut microbes of rodents as they consume different foods, which in turn affects their growth [64]. In the steppe of Inner Mongolia, climate warming could exacerbate the degraded grassland and an additional increase in temperature may further reduce the plant species preferred by rodents [20]. In the present study, considering that precipitation was significantly correlated with temperature (Supplementary Materials Figure S1), increases in one or both of these parameters could initially improve habitat conditions in the examined semiarid areas, and thus enhance survival of Mongolian gerbils, maybe through bottom-up trophic processes. Nevertheless, temperatures above optimal degrees may decrease the demography of small mammals. Zhang et al. (2003) reported that the numerical response of Brandt's voles was mainly determined by the availability of habitats with sparse, short grasses in our study area [20]. Brandt's vole populations do not persist in habitats with high plant cover (>80%) and tall (>17-20 cm) grasses [65]. Chen et al. (2015) also reported that climatic factors, especially the increase in precipitation, could have a nonlinear effect on the survival probability of Daurian pika populations [29]. Compared to Brandt's voles and Daurian pikas, Mongolian gerbils prefer habitats with shorter, sparser vegetation [66]. Additionally, they are territorial group-living granivorous rodents, and the size of a social group is positively related to the area in which food resources can be found [36]. Increased temperature and precipitation, especially in spring and summer, were shown to enhance plant coverage, height, and productivity, thus improving population growth through bottom-up trophic [22]. However, high temperatures and precipitation beyond the optimal condition may result in gerbil habitats being dominated by tall and dense vegetation, and the food availability of soil seed banks was significantly reduced due to seed germination [20]. A reduction in suitable or optimal habitats and/or an increase in the incidence of disease might delay the growth and affect the apparent survival of populations in Mongolian gerbils. These indirect effects of climate may also have significant impacts on the population growth rate of gerbils.
In addition, as increases in population density had stronger negative effects on recruitment than on survival (Table 1), as well as on population growth rate in Mongolian gerbils (Table 2), it was suggested that the population growth rates of small mammals are more sensitive to changes in recruitment than to those in survival. Generally, the recruitment of rodents is contributed to birth (involving new born, infant survival) and immigration [33]. Direct density dependence may reduce body mass, reproduction, or immigration with increasing densities and subsequently stabilize population abundances [32]. Thus, the present study suggested that density dependence served as a regulatory mechanism of population growth rate of gerbils mainly through recruitment, and not survival. Densitydependent feedbacks to the recruitment of gerbils, on the other hand, may also serve as the variation in life strategies with a demographic stochasticity to improve the gerbils' ecological adaption, as observed in other small mammals [29,32]. Long-term and/or multiple-site mark-recapture studies are needed to better understand population dynamic patterns in Mongolian gerbils. The complex responses of populations to temperature, as a limiting ecological factor, may raise concerns for the future fate of gerbils within the context of global climate change scenarios [67,68], which predict potential increases in temperature and precipitation with the marked seasonal change at northern latitude zones.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14080586/s1, Figure S1. Correlation between the average daily temperature and precipitation of each trapping interval between two successive trapping sessions (Pearson's r = 0.46, p = 0.002). Table S1 Author Contributions: W.L. designed the experiment; W.L. and K.D. contributed to the experiment and the analysis of the data; W.L. drafted the manuscript; W.L. and K.D. revised the manuscript. All authors critically revised the manuscript, agreed to be fully accountable for ensuring the integrity and accuracy of the work, and read and approved the final manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request. The raw sequence data are available in the NCBI Sequence Read Archive under the accession number (PRJNA818461).