Data analysis
All the statistics were conducted in R (v4.0.2, RCoreTeam 2022). We assessed how corridor quality, drought severity, and number of drought patches influence the persistence ability of metapopulations in three ways. First, we used a generalized linear model (GLM) with a log-link gaussian error to test how those factors influence the time of population extinction. The time of population extinction was defined as the midpoint between the last week which individuals were observed alive in at least one of the four patches of a landscape, and the first time no individuals were observed as being alive. Only three microcosms were alive at the end of experiment (week 16, Fig. 2), but no individuals were observed alive the week after that. Experimental factors (i.e., corridor quality, drought severity, number of drought patches) and all their two-way interactions were included as explanatory variables. As populations were stable and not observed any extinctions until the end of the experiment (Fig. 2), two constant treatments were removed from further analysis.
Second, we used a generalized linear mixed effect model (GLMM) with a zero-inflated negative binomial distribution to investigate how corridor quality, drought severity, and number of drought patches influence the variability of abundance among patches through time. The variability of abundance among patches was defined as the standard deviation of abundances between the four habitat patches. Microcosm identity was included as a random effect to account for the potential differences between microcosms including their starting density. We included time, corridor quality, drought severity, number of drought patches, and their interactions as fixed effects. Of particular interest was how experimental factors (i.e., corridor quality, drought severity, number of drought patches) and time interactively affect the changes in variability, thus we also include their two-way interactions in the fixed effect. Time and experimental factors (i.e., corridor quality, drought severity, number of drought patches) were included as a zero-inflation parameter to account for the overdispersion of zero counts after extinction occurred. Data collected before drought treatment started (i.e., from week 1 to 3) was excluded from model fitting. The model was fitted using “glmmTMB’ package (Magnusson et al. 2017).
Finally, we used nonlinear regression curves with a three-parameter logistic function to fit the changes in total population abundance in each habitat network over time. The model was specified as Y= d / (1+ exp (- b * (X - e))) where Y was total abundance and X was time; d the high asymptote, b the maximum slope and e the time at the maximum slope, were all estimated by models. This is a three-parameter logistic function, and we chose this model due to its simplicity and utility (Kingsland 1982). Models were fitted using ‘drc’ package (Ritz et al. 2015). To explore the effect of treatments on the maximum rate of population change, we extracted the estimated value of b in each model and compared across treatments (Fig. S1). We used a GLM with a gaussian distribution to fit the log-transformed data, where corridor quality, drought severity, number of drought patches, and their two-way interactions was included as predictors.