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.