Connectivity among thermal habitats buffers the effects of warm climate on life-history traits and population dynamics

1. Contemporary climate change affects population dynamics, but its influence varies with landscape structure. It is still unclear whether landscape fragmentation buffers or enhances the effects of climate on population size and on the age and body size of individuals composing these populations. 2. This study aims to investigate the impacts of warm climates on lizard life-history traits and population dynamics in habitats varying in their connectivity. 3. We monitored common lizard (Zootoca vivipara) populations for three years in an experimental system in which both climatic conditions and connectivity among habitats were simultaneously manipulated. We considered two climatic treatments (i.e., present-day climate and warm climate (+1.4°C than present-day climate)) and two connectivity treatments (i.e., a connected treatment in which individuals could move from one climate to the other and an isolated treatment in which movement between climates was not possible). We monitored survival, reproduction, growth, dispersal, age and body size of each individual in the system as well as population density through time. 4. We found that the influence of warm climates on the life-history traits and population dynamics depended on the connectivity among thermal habitats. Populations in warm climates were i) composed of younger individuals only when isolated; ii) larger in population size only in connected habitats; and iii) composed of larger age-specific individuals independently of the landscape configuration. The connectivity among habitats altered population responses to climate warming likely through asymmetries in the flow and phenotype of dispersers between thermal habitats. 5. Our results demonstrate that landscape fragmentation can drastically change the dynamics and persistence of populations facing climate change.Excel/LibreOffice Calc, a text editor and R are needed to open the files and run the scripts.Funding provided by: Agence Nationale de la RechercheCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100001665Award Number: ANR-11-INBS-0001AnaEE-ServicesFunding provided by: H2020 European Research CouncilCrossref Funder Registry ID: http://dx.doi.org/10.13039/100010663Award Number: 817779Funding provided by: Fondation FyssenCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100003135Award Number: Post-Doctoral FellowshipFunding provided by: French Laboratory of Excellence project 'TULIP'*Crossref Funder Registry ID: Award Number: ANR-10-LABX-41Funding provided by: French Laboratory of Excellence project 'TULIP'*Crossref Funder Registry ID: Award Number: ANR-11-IDEX-0002-02Experiments Starting in 2015, we performed a three-year-long experiment on the common lizard (Zootoca vivipara) using the Metatron, an experimental system comprising 48 interconnected semi-natural mesocosms (100 m2 each (Legrand et al., 2012), Fig. 1b) connected by 19-meter-long corridors which allowed us to simultaneously manipulate temperature (with the help of automatic shutters) and connectivity (by opening or closing the doors which separate the mesocosm from the corridor, see details in Fig. 1 and Appendix A, ethics permit number APAFIS#19523-201902281559649 v3). We created two climatic treatments, a present-day and a warm climate, by automatically closing the shutters at ambient temperature thresholds of either 28°C or 38°C (Bestion, Teyssier, et al., 2015). Given that mesocosms are intrinsically warmer than outside, present-day climates matched thermal conditions near the Metatron (meteorological station of Saint-Girons Antichan (Bestion, Teyssier, et al., 2015)). As our treatments depended on outdoor climatic conditions, the generated climate regimes follow day-to-day fluctuations in a coordinated manner so that daily fluctuation and seasonality are efficiently reproduced. The climate manipulations were active and efficient during the summer daytime, and the difference between treatments varied with the weather (Fig. 1c). Consequently, the warm climate treatment was on average 1.4 and 2.6°C warmer (mean and maximal summer daily temperatures) than the present-day climate treatment and the summer temperature slightly differed among years (Fig. 1c). The warm climate treatment matched the scenario SSP5-8.5 for 2041-2060 and SSP2-4.5, SSP3-7.0 and SSP5-8.5 for 2081-2100 (Masson-Delmotte, V., P. Zhai, A. Pirani et al., 2021). Of the 48 mesocosms constituting the Metatron, we used 16 to create eight pairs of mesocosms that combined a present-day climate mesocosm and a warm climate mesocosm with two connectivity treatments, either opening or closing the connecting corridors (Fig. 1a). For four pairs, lizards could move between climate treatments (connected treatment), while movement was prevented for the four remaining pairs (isolated treatment). In the connected treatments, corridors were opened from early March to mid-October each year, except in 2017 when corridors opened in late March, spanning the entire period of lizard activity. Population initialization In 2015, 240 adults and 306 juveniles were released into 16 mesocosms. The individuals were descendants of lizards captured in the Cevennes, France, in 2010 and 2013, maintained in the Metatron (Ariège, France) for several experiments (Bestion et al., 2017, 2019; Bestion, Teyssier, et al., 2015) and intermixed regularly before the present experiment to prevent high levels of inbreeding. Each mesocosm was initially populated in early July 2015 with adult females, males and juveniles, with a later addition of 4 adults (2 males, 2 females) in September 2015. The population size after this initialization period was 10 females, 5 males and 19 ± 1 juveniles, matching densities observed in natural populations. The individuals added in September were not considered in the analyses for the first year. At population initialization only, we split clutches among different mesocosms, and the different treatments to enhance genetic diversity within populations and released juveniles without their mother to prevent kin competition. All the lizards present in the system, individually tagged at birth by toe clipping (see the ethics statement section, ethics permit number APAFIS#15897-2018070615164391 v3), were therefore of known age. We measured their body size (snout-vent length) at release and ensured that there were no differences in age structure and body size between treatments (p-values > 0.63 for the effects of climate, connectivity treatments and their interaction). Population monitoring In May, from 2016 to 2018, before the females started laying eggs, we closed the corridors to recapture all the individuals and brought them to the laboratory. They were identified, measured for body size and mass and maintained in individual terraria (18x35x22 cm for adult females and gravid one-year-old females and 15.5x25x15 cm for males and non-gravid one-year-old females). Terraria contained a 3 cm sterilized litter layer, a petri dish with water, a piece of absorbent paper, and a cardboard and a plastic tube as a shelter. A light bulb (25 W) and an ultraviolet lamp (Zoomed Reptisun 5.0 UVB 36 W) provided heat for thermoregulation and light 6 h per day (from 9:00 to 12:00 and from 14:00 to 17:00). Lizards were lightly sprayed with water three times a day (in the morning, at mid-day, and in the evening) and offered two crickets (Acheta domestica) daily. Females laid eggs in their terrarium and the juveniles were isolated from their mother directly after parturition. They were marked by toe-clipping, measured for snout-vent length (nearest mm), weighed with an electronic scale (Thermo Fisher, 0.01g) and a tail tip of 0.3 cm was collected for paternity analysis (ethics permit number APAFIS#19523-201902281559649 v3). These captures allowed (i) monitoring population size, mean age and body size of individuals that composed the populations through time,(ii) measuring clutch sizes (i.e. number of viable offspring), yearly survival probability, body growth (i.e., difference in snout-vent length from one year compared to the previous one), and (iii), in connected treatments, the dispersal status of each individual every experimental year. Dispersers (respectively, residents) are defined as individuals recaptured after one year in a different (respectively, the same) mesocosm than the previous year. The dispersal probability observed in this experiment was similar to the one of a previous experiment (Legrand et al., 2012) with similar experimental conditions and in which movements were recorded daily, meaning that our one-year movement monitoring provided a good indication of individual dispersal status. Nevertheless, seasonal back-and-forth movements might have occurred, as it has been previously estimated that 47% of individuals that moved at least once moved a second time in a course of a year (Legrand et al., 2012). Early July, all males, females and their clutch were released into the Metatron into their population of origin. We released adult and juvenile individuals back into the mesocosm where they or their mother were captured. To avoid stress-induced dispersal, we closed the corridors for 7 days after release and opened them later. Over the course of the experiment, two populations went extinct in 2016 (one of each climate treatment in the isolated treatment) and one in 2017 (from present-day climate in the isolated treatment). Extinctions were likely due to demographic fluctuations, as extinctions occurred only in isolated populations where population rescue from neighbouring populations was not possible. In 2016, the two extinct populations were reinitialized with the same density, age and sex structure and phenotypic composition as in 2015 using lizards from stock mesocosms that were not the subject of any previous climatic experiment. In 2017, the extinct population was not reinitialized. Statistical analyses General statistical procedure We analyzed the additive and interactive effects of climate and habitat connectivity first on life-history traits, then on individuals' mean age and body size and finally on population size. We further investigated the influence of climate in each connectivity condition, by separately running the models for each connectivity treatment. All the models are summarized in Tables S1 & S2. All continuous variables were centered and scaled in all analyses. We used generalized/linear mixed models with random intercepts to account for the non-independency of the data points. We proceeded in two steps. First, we built full models with all fixed variables and random effects, and the random structure of each model was selected by AIC (Zuur et al., 2009). Second, we built models with the selected random effect(s) and with all possible combinations of fixed effects, including one without any fixed effects (i.e., null model), and ranked them by AIC. We obtained conditional estimates, standard errors, z-values, relative importance (RI), and p-values of all variables that feated in those models that had a delta AIC of less than 2 from the best model using a model averaging procedure (Burnham et al., 2011). Models containing all variables present in the averaged best models were run to calculate the conditional (R2c, effect of the fixed and random effects) and marginal (R2m, effect of the fixed effect) R2 (Nakagawa et al., 2017). When the best model structure was that without random effect, the adjusted R2 was calculated, except for models with zero-inflated distribution, where the pseudo-R2 was calculated. We provided both RI and p-value, and chose to discuss the influence of variables when their RI, p-value and visual pattern on the figures were consistent (i.e. high RI, low p-value and clear visual pattern), without fixing absolute thresholds. Accordingly, when a treatment was retained in the best average model with low RI and high p-value, its effect was interpreted as weak. All analyses were performed using R version 4.0.5 (R Core Team, 2021) and necessitate lme4 (Bates et al., 2015), glmmADMB (Fournier et al., 2012; Skaug et al., n.d.), performance (Lüdecke et al., 2021), MuMIn (Barton, 2020), DHARMa (Hartig, 2021), and emmeans (Lenth, 2021) R packages. Life-history traits All individuals older than one year were considered adults and analyzed together, while younger individuals (hereafter named juveniles) were analyzed separately. We analyzed the clutch size, yearly survival probability and body growth rate (Table S1). The clutch sizes were analyzed separately for each sex because (i) the reproductive strategy of each sex may differ with respect to the experimental treatments and (ii) to avoid testing for four-way interactions (i.e., climate*connectivity*time*sex). The clutch sizes of males were obtained via paternity analyses and all males' results are presented in appendix B. We used generalized mixed models with a binomial distribution for survival, with zero-inflated Poisson distribution for the clutch size and linear mixed models for body growth. All models included climate treatments, connectivity treatments, and time as a continuous variable and their three-way interaction. The models for the clutch size of juveniles did not converge with the three-way interaction due to a low number of reproductive juveniles and so this interaction was removed. Models further included covariates known to influence life-history traits: body size (e.g. Cotto et al., 2015) for all analyses, sex (e.g. Bestion, Teyssier, et al., 2015) for survival and body growth analyses, and birth date (e.g. Bestion, Teyssier, et al., 2015) in Julian days in the analyses on juveniles. Random intercepts included mesocosm identity, individual identity for analyses on adults, and family identity for the analyses on juveniles as siblings were not independent. In the connected treatment, juvenile and adult dispersal statuses (i.e., disperser or resident, N=116 for juveniles and N=113 for adults) were analyzed (Table S1). The models included climate treatments, time, their interaction, body size and the interaction between body size and climate treatments, as body size strongly influences dispersal decisions, costs and benefits (Cote et al., 2007). Random structure only included mesocosm identity because models did not converge when individual identity was included, due to the low number of dispersers. Note that the mesocosm identity of an individual could change from one year to the other if the given individual dispersed. Mean age, body size and population size We used generalized mixed models with Poisson distribution to analyze population size (N=47) and mean age of individuals composing the populations (N=617), and linear mixed models for individuals' mean body size (N=617). For all models, the fixed effects were climate treatment, connectivity treatment, number of years since population initialization (hereafter referred to as "time") and their three-way interaction. The fixed time effect estimates the temporal pattern of population dynamics, rather than controlling for temporal autocorrelation in our data. Nevertheless, we checked for potential temporal autocorrelation in the residuals of the final models (i.e. Durbin-Watson test and ACF plot). All tests revealed no autocorrelation in the residuals. All the data at time 0 (i.e. before climate treatment) were excluded from the analyses and when a population was reinitialized after extinction, time was set at 0. Models for population size also included mesocosm identity as a random intercept to take into account the repeated model structure (Table S2). For individuals' mean age and body size (Table S2), the models were run at the individual level and included mesocosm and individual identities as random intercepts to account for the non-independency of individuals of the same population and the multiple occurrences of individuals over time. We also included sex and population density as fixed effects to control for the indirect effects of the treatments on age and body size structure through their effects on sex ratio and density. The age of the individuals was also included in the body size analysis to disentangle the direct effect of climate on body size from its indirect effect through age structure changes. The same model without age was additionally run (Table S7). Ethics statement: The Station d'Ecologie Théorique et Expérimentale has a national agreement for use of animals in the laboratory (number B09583), and our experiments are made in accordance with French ethics regulations (Ethics permits number APAFIS#15897-2018070615164391 v3 for toe clipping and APAFIS#19523-201902281559649 v3 for other experimental procedures, including the maintenance of lizards in the Metatron). The lizards were initially captured in the wild under license numbers 2010-189-16 DREAL and 2013-274-0002