Philippa Kaur

Sampling of coastal communitiesHere, we integrated data from five different projects that had surveyed coastal communities across five countries47,48,49,50. Between 2009 and 2015, we conducted socioeconomic surveys in 72 sites from Indonesia (n = 25), Madagascar (n = 6), Papua New Guinea (n = 10), the Philippines (n = 25), and Tanzania (Zanzibar) (n = 6). Site selection was for broadly similar purposes- to evaluate the effects of various coastal resource management initiatives (collaborative management, integrated conservation and development projects, recreational fishing projects) on people’s livelihoods in rural and peri-urban villages. Within each project, sites were purposively selected to be representative of the broad range of socioeconomic conditions (e.g., population size, levels of development, integration to markets) experienced within the region. We did not survey strictly urban locations (i.e., major cities). Because our sampling was not strictly random, care should be taken when attempting to make inferences beyond our specific study sites.We surveyed between 13 and 150 households per site, depending on the population of the communities and the available time to conduct interviews per site. All projects employed a comparable sampling design: households were either systematically (e.g., every third house), randomly sampled, or in the case of three villages, every household was surveyed (a census) (see Supplementary Data file). Respondents were generally the household head, but could have been other household members if the household head was not available during the study period (i.e. was away). In the Philippines, sampling protocol meant that each village had an even number of male and female respondents. Respondents gave verbal consent to be interviewed.The following standard methodology was employed to assess material style of life, a metric of material assets-based wealth48,51. Interviewers recorded the presence or absence of 16 material items in the household (e.g., electricity, type of walls, type of ceiling, type of floor). We used a Principal Component Analysis on these items and kept the first axis (which explained 34.2% of the variance) as a material wealth score. Thus, each community received a mean material style of life score, based on the degree to which surveyed households had these material items, which we then scaled from 0 to 1. We also conducted an exploratory analysis of how material style of life has changed in two sites in Papua New Guinea (Muluk and Ahus villages) over fifteen and sixteen-year time span across four and five-time periods (2001, 2009, 2012, 2016, and 2002, 2009, 2012, 2016, 2018), respectively, that have been surveyed since 2001/200252. These surveys were semi-panel data (i.e. the community was surveyed repeatedly, but we did not track individuals over each sampling interval) and sometimes occurred in different seasons. For illustrative purposes, we plotted how these villages changed over time along the first two principal components.SensitivityWe asked each respondent to list all livelihood activities that bring in food or income to the household and rank them in order of importance. Occupations were grouped into the following categories: farming, cash crop, fishing, mariculture, gleaning, fish trading, salaried employment, informal, tourism, and other. We considered fishing, mariculture, gleaning, fish trading together as the ‘fisheries’ sector, farming and cash crop as the ‘agriculture’ sector and all other categories into an ‘off-sector’.We then developed three distinct metrics of sensitivity based on the level of dependence on agriculture, fisheries, and both sectors together. Each metric incorporates the proportion of households engaged in a given sector (e.g., fisheries), whether these households also engage in occupations outside of this sector (agriculture and salaried/formal employment; referred to as ‘linkages’ between sectors), and the directionality of these linkages (e.g., whether respondents ranked fisheries as more important than other agriculture and salaried/formal employment) (Eqs. 1–3)$${{{{{{rm{S}}}}}}}_{{{{{{rm{A}}}}}}}=,frac{{{{{{rm{A}}}}}}}{{{{{{rm{A}}}}}}+{{{{{rm{NA}}}}}}},times ,frac{{{{{{rm{N}}}}}}}{{{{{{rm{A}}}}}}+{{{{{rm{NA}}}}}}},times ,frac{left(frac{{{{{{{rm{r}}}}}}}_{{{{{{rm{a}}}}}}}}{2}right),+,1}{{{{{{{rm{r}}}}}}}_{{{{{{rm{a}}}}}}}+,{{{{{{rm{r}}}}}}}_{{{{{{rm{na}}}}}}}+1}$$
(1)
$${{{{{{rm{S}}}}}}}_{{{{{{rm{F}}}}}}}=,frac{{{{{{rm{F}}}}}}}{{{{{{rm{F}}}}}}+{{{{{rm{NF}}}}}}},times ,frac{{{{{{rm{N}}}}}}}{{{{{{rm{F}}}}}}+{{{{{rm{NF}}}}}}},times ,frac{left(frac{{{{{{{rm{r}}}}}}}_{{{{{{rm{f}}}}}}}}{2}right),+,1}{{{{{{{rm{r}}}}}}}_{{{{{{rm{f}}}}}}}+,{{{{{{rm{r}}}}}}}_{{{{{{rm{nf}}}}}}}+1}$$
(2)
$${{{{{{rm{S}}}}}}}_{{{{{{rm{AF}}}}}}}=,frac{{{{{{rm{AF}}}}}}}{{{{{{rm{AF}}}}}}+{{{{{rm{NAF}}}}}}},times ,frac{{{{{{rm{N}}}}}}}{{{{{{rm{AF}}}}}}+{{{{{rm{NAF}}}}}}},times ,frac{left(frac{{{{{{{rm{r}}}}}}}_{{{{{{rm{af}}}}}}}}{2}right),+,1}{{{{{{{rm{r}}}}}}}_{{{{{{rm{af}}}}}}}+,{{{{{{rm{r}}}}}}}_{{{{{{rm{naf}}}}}}}+1}$$
(3)
where ({{{{{{rm{S}}}}}}}_{{{{{{rm{A}}}}}}}), ({{{{{{rm{S}}}}}}}_{{{{{{rm{F}}}}}}}) and ({{{{{{rm{S}}}}}}}_{{{{{{rm{AF}}}}}}}) are a community’s sensitivity in the context of agriculture, fisheries and both sectors, respectively. A, F and AF are the number of households relying on agriculture-related occupations within that community, fishery-related and agriculture- and fisheries-related occupations within the community, respectively. NA, NF and NAF are the number of households relying on non-agriculture-related, non-fisheries-related, and non-agriculture-or-fisheries-related occupations within the community, respectively. N is the number of households within the community. ({{{{{{rm{r}}}}}}}_{{{{{{rm{a}}}}}}}), ({{{{{{rm{r}}}}}}}_{{{{{{rm{f}}}}}}}) and ({{{{{{rm{r}}}}}}}_{{{{{{rm{af}}}}}}}) are the number of times agriculture-related, fisheries-related and agriculture-and-fisheries-related occupations were ranked higher than their counterpart, respectively. ({{{{{{rm{r}}}}}}}_{{{{{{rm{na}}}}}}}), ({{{{{{rm{r}}}}}}}_{{{{{{rm{nf}}}}}}}) and ({{{{{{rm{r}}}}}}}_{{{{{{rm{naf}}}}}}}) are the number of times non-agriculture, non-fisheries, and non-agriculture-and-fisheries-related occupations were ranked higher than their counterparts. As with the material style of life, we also conducted an exploratory analysis of how joint agriculture-fisheries sensitivity has changed over time in a subset of sites (Muluk and Ahus villages in Papua New Guinea) that have been sampled since 2001/200252. Although our survey methodology has the potential for bias (e.g. people might provide different rankings based on the season, or there might be gendered differences in how people rank the importance of different occupations53), our time-series analysis suggest that seasonal and potential respondent variation do not dramatically alter our community-scale sensitivity metric.ExposureTo evaluate the exposure of communities to the impact of future climates on their agriculture and fisheries sectors, we used projections of production potential from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) Fast Track phase 3 experiment dataset of global simulations. Production potential of agriculture and fisheries for each of the 72 community sites and 4746 randomly selected sites from our study countries with coastal populations >25 people/km2 were projected to the mid-century (2046–2056) under two emission scenarios (SSP1-2.6, and SSP5-8.5) and compared with values from a reference historical period (1983–2013).For fisheries exposure (EF), we considered relative change in simulated total consumer biomass (all modelled vertebrates and invertebrates with a trophic level >1). For each site, the twenty nearest ocean grid cells were determined using the Haversine formula (Supplementary Fig. 5). We selected twenty grid cells after a sensitivity analysis to determine changes in model agreement based on different numbers of cells used (1, 3, 5, 10, 20, 50, 100; Supplementary Figs. 6–7), which we balanced off with the degree to which larger numbers of cells would reduce the inter-site variability (Supplementary Fig. 8). We also report 25th and 75th percentiles for the change in marine animal biomass across the model ensemble. Projections of the change in total consumer biomass for the 72 sites were extracted from simulations conducted by the Fisheries and marine ecosystem Model Intercomparison Project (FishMIP3,54). FishMIP simulations were conducted under historical, SSP1-2.6 (low emissions) and SSP5-8.5 (high emissions) scenarios forced by two Earth System Models from the most recent generation of the Coupled Model Intercomparison project (CMIP6);55 GFDL-ESM456 and IPSL-CM6A-LR57. The historical scenario spanned 1950–2014, and the SSP scenarios spanned 2015–2100. Nine FishMIP models provided simulations: APECOSM58,59, BOATS60,61, DBEM2,62, DBPM63, EcoOcean64,65, EcoTroph66,67, FEISTY68, Macroecological69, and ZooMSS11. Simulations using only IPSL-CM6A-LR were available for APECOSM and DBPM, while the remaining 7 FishMIP models used both Earth System Model forcings. This resulted in 16 potential model runs for our examination of model agreement, albeit with some of these runs being the same model forced with two different ESMs. Thus, the range of model agreement could range from 8 (half model runs indicating one direction of change, and half indicating the other) to 16 (all models agree in direction of change). Model outputs were saved with a standardised 1° spatial grid, at either a monthly or annual temporal resolution.For agriculture exposure (EA), we used crop model projections from the Global Gridded Crop model Intercomparison Project (GGCMI) Phase 314, which also represents the agriculture sector in ISIMIP. We used a window of 11×11 cells centred on the site and removed non-land cells (Supplementary Fig. 5). The crop models use climate inputs from 5 CMIP6 ESMs (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL), downscaled and bias-adjusted by ISIMIP and use the same simulation time periods. We considered relative yield change in three rain-fed and locally relevant crops: rice, maize, and cassava, using outputs from 4 global crop models (EPIC-IIASA, LPJmL, pDSSAT, and PEPIC), run at 0.5° resolution. These 4 models with 5 forcings generate 20 potential model runs for our examination of model agreement. Yield simulations for cassava were only available from the LPJmL crop model. All crop model simulations assumed no adaptation in growing season and fertilizer input remained at current levels. Details on model inputs, climate data, and simulation protocol are provided in ref. 14. At each site, and for each crop, we calculated the average change (%) between projected vs. historical yield within 11×11 cell window. We then averaged changes in rice, maize and cassava to obtain a single metric of agriculture exposure (EA).We also obtained a composite metric of exposure (EAF) by calculating each community’s average change in both agriculture and fisheries:$${{{{{{rm{E}}}}}}}_{{{{{{rm{AF}}}}}}}=,frac{{{{{{{rm{E}}}}}}}_{{{{{{rm{A}}}}}}}+,{{{{{{rm{E}}}}}}}_{{{{{{rm{F}}}}}}}}{2}$$
(4)
Potential ImpactWe calculated relative potential impact as the Euclidian distance from the origin (0) of sensitivity and exposure.Sensitivity testTo determine whether our sites displayed a particular exposure bias, we compared the distributions of our sites and 4746 sites that were randomly selected from 47,460 grid cells within 1 km of the coast of the 5 countries we studied which had population densities >25 people/km2, based on the SEDAC gridded populating density of the world dataset (https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev11/data-download).We used Cohen’s D to determine the size of the difference between our sites and the randomly selected sites.Validating ensemble modelsWe attempted a two-stage validation of the ensemble model projections. First, we reviewed the literature on downscaling of ensemble models to examine whether downscaling validation had been done for the ecoregions containing our study sites.While no fisheries ensemble model downscaling had been done specific to our study regions, most of the models of the ensemble have been independently evaluated against separate datasets aggregated at scales down to Large Marine Ecosystems (LMEs) or Exclusive Economic Zones (EEZs) (see11). For example, the DBEM was created with the objective of understanding the effects of climate change on exploited marine fish and invertebrate species2,70. This model roughly predicts species’ habitat suitability; and simulates spatial population dynamics of fish stocks to output biomass and maximum catch potential (MCP), a proxy of maximum sustainable yield2,62,71. Compared with spatially-explicit catch data from the Sea Around Us Project (SAUP; www.seaaroundus.org)70 there were strong similarities in the responses to warming extremes for several EEZs in our current paper (Indonesia and Philippines) and weaker for the EEZs of Madagascar, Papua New Guinea, and Tanzania. At the LME level, DBEM MCP simulations explained about 79% of the variation in the SAUP catch data across LMEs72. The four LMEs analyzed in this paper (Agulhas Current; Bay of Bengal; Indonesian Sea; and Sulu-Celebes Sea) fall within the 95% confidence interval of the linear regression relationship62. Another example, BOATS, is a dynamic biomass size-spectrum model parameterised to reproduce historical peak catch at the LME scale and observed catch to biomass ratios estimated from the RAM legacy stock assessment database (in 8 LMEs with sufficient data). It explained about 59% of the variability of SAUP peak catch observation at the LME level with the Agulhas Current, Bay of Bengal, and Indonesian Sea catches reproduced within +/-50% of observations61. The EcoOcean model validation found that all four LMEs included in this study fit very close to the 1:1 line for overserved and predicted catches in 200064,65. DBPM, FEISTY, and APECOSM have also been independently validated by comparing observed and predicted catches. While the models of this ensemble have used different climate forcings when evaluated independently, when taken together the ensemble multi-model mean reproduces global historical trends in relative biomass, that are consistent with the long term trends and year-on-year variation in relative biomass change (R2 of 0.96) and maximum yield estimated from stock assessment models (R2 of 0.44) with and without fishing respectively11.Crop yield estimates simulated by GGCMI crop models have been evaluated against FAOSTAT national yield statistics14,73,74. These studies show that the models, and especially the multi-model mean, capture large parts of the observed inter-annual yield variability across most main producer countries, even though some important management factors that affect observed yield variability (e.g., changes in planting dates, harvest dates, cultivar choices, etc.) are not considered in the models. While GCM-based crop model results are difficult to validate against observations, Jägermeyr et al14. show that the CMIP6-based crop model ensemble reproduces the variability of observed yield anomalies much better than CMIP5-based GGCMI simulations. In an earlier crop model ensemble of GGCMI, Müller et al.74 show that most crop models and the ensemble mean are capable of reproducing the weather-induced yield variability in countries with intensely managed agriculture. In countries where management introduces strong variability to observed data, which cannot be considered by models for lack of management data time series, the weather-induced signal is often low75, but crop models can reproduce large shares of the weather-induced variability, building trust in their capacity to project climate change impacts74.We then attempted to validate the models in our study regions. For the crop models, we examined production-weighted agricultural projections weighted by current yields/production area (Supplementary Fig. 1). We used an observational yield map (SPAM2005) and multiplied it with fractional yield time series simulated by the models to calculate changes in crop production over time, which integrates results in line with observational spatial patterns. The weighted estimates were not significantly different to the unweighted ones (t = 0.17, df = 5, p = 0.87). For the fisheries models, our study regions were data-poor and lacked adequate stock assessment data to extend the observed global agreement of the sensitivity of fish biomass to climate during our reference period (1983-2013). Instead, we provide the degree of model run agreement about the direction of change in the ensemble models to ensure transparency about the uncertainty in this downscaled application.AnalysesTo account for the fact that communities were from five different countries we used linear mixed-effects models (with country as a random effect) for all analyses. All averages reported (i.e. exposure, sensitivity, and model agreement) are estimates from these models. In both our comparison of fisheries and agriculture exposure and test of differences between production-weighted and unweighted agriculture exposure we wanted to maintain the paired nature of the data while also accounting for country. To accomplish this we used the differences between the exposure metrics as the response variable (e.g. fisheries exposure minus agriculture exposure), testing whether these differences are different from zero. We also used linear mixed-effects models to quantify relationships between the material style of life and potential impacts under different mitigation scenarios (SSP1-2.6 and 8.5), estimating standard errors from 1000 bootstrap replications. To further explore whether these relationships between the material style of life and potential impacts were driven by exposure or sensitivity, we conducted an additional analysis to quantify relationships between the material style of life and: 1) joint fisheries and agricultural sensitivity; 2) joint fisheries and agricultural exposure under different mitigation scenarios. We present both the conditional R2 (i.e., variance explained by both fixed and random effects) and the marginal R2 (i.e., variance explained by only the fixed effects) to help readers compare among the material style of life relationships.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Testing H1 and H2 at community composition levelAs noted above, the simple fact that fungi grow more slowly than bacteria is the basis of the hypotheses that (H1) fungal communities should be more resistant than bacterial communities to drought stress, and (H2) that fungal communities should be less resilient than bacterial communities when the stress is relieved by rewetting18. In addition to growth rate, these two hypotheses may be related to differences in the form of growth between fungi and bacteria. For example, multicellular hyphal growth versus unicellular division or the greater thickness of fungal cell walls as compared to those of bacteria47,48. We tested H1 and H2 at the community composition level by blending the fungal and bacterial datasets generated from the same leaf, root, rhizosphere and soil samples collected from field-grown sorghum that had been either irrigated as a control, or subjected to preflowering drought followed by regular wetting beginning at flowering10,11.We followed the approach of Shade et al.17 to detect resistance and resilience, which had been developed for univariate variables, e.g., richness. For multivariate data, e.g., community composition, we modified it by calculating pairwise community dissimilarity for two groups: within-group (control-control pairs, drought-drought pairs, or rewetting-rewetting pairs), and between-group (control-drought pairs, or control-rewetting pairs). Ecological resistance to drought stress is detected by comparing compositional dissimilarity of between-group pairs (control-drought pairs) against within-group pairs (control-control pairs and drought-drought pairs) for each of the droughted weeks (weeks 3–8). Ecological resilience to rewetting is detected by assessing, from before to after rewetting, the change in the difference of compositional dissimilarity between within-group pairs and between-group pairs. Here, the point just before rewetting was week 8 and the points after rewetting were weeks 9–17. A t-test was used to assess the statistical significance of the differences in resistance or resilience between bacterial and fungal communities at each time point for each compartment.To account for the different resolutions of ITS and 16 S, we compared bacterial 16 S OTUs against both fungal ITS, species-level OTUs as well the fungal family level (Supplementary Fig. 1). The results of analyses using either fungal families or OTUs are consistent. Out of 36 comparisons (15 roots, 15 rhizospheres and 6 soils), different family and OTUs results were detected in four instances. In two of these, significances detected by OTUs were not detected by family (root, weeks 4 and 17) and, in the other two cases, significances detected by family were not detected by OTUs (rhizosphere, weeks 7 and 8). (Fig. 1). We report only results that are consistent at both the species and family levels (Fig. 1).In line with our first hypothesis, H1, we found that the resistance to drought stress for fungal mycobiomes was consistently stronger than that for bacterial microbiomes for weeks 5 in root, weeks 4–6 in rhizosphere, and weeks 4 and 6–8 in soil (Fig. 1, Supplementary Fig. 1 and Supplementary Table 2). In support of our second hypothesis, H2, when the stress of pre-flowering drought was relieved by rewetting, we found that the resilience of the bacterial communities was consistently higher than that for the fungi in weeks 9–16 in root, and weeks 11–17 in rhizosphere (Fig. 1, Supplementary Fig. 1 and Supplementary Table 2).Surprisingly, we found that resilience was stronger for fungal than bacterial communities in the first week (week 9) of rewetting in the rhizosphere (Fig. 1, Supplementary Fig. 1 and Supplementary Table 2). This high resilience of fungi may be associated with the quick growth of sorghum roots when rewetted. The rhizosphere zone around these newly formed roots may be quickly colonized by soil fungi, a community that was weakly affected by drought. This result suggests that re-assembly of the rhizosphere microbial community is more complex than previously expected.The finding that fungal community composition in the soil is not shaped by drought prevented us from further detecting resilience (Fig. 1). Note fungal community in early leaves was excluded from analysis due to the high proportion of non-fungal reads in sequencing11.Testing H1 and H2 at all-correlation levelNext, we moved from the comparison of whole communities to correlation among individual bacterial and fungal taxa to test the hypotheses about resistance, H1, and resilience, H2. As noted above, previous research provided the foundation for the stress gradient hypothesis, which predicts an increase in positive associations in stress32,33,34,35,36,37. Further, ecological modeling predicts that negative associations promote stability40. Concerning specific associations, studies of Arabidopsis and associated microbes reported that positive associations are favored within kingdoms, i.e., within bacteria or within fungi, while negative associations predominate between kingdoms38,39. Given these foundations, concerning H1, we expected an increase in the proportion of positive correlation by drought stress that would be strongest for B-B, followed by F-F, and lastly by B-F; for H2 we expected rewetting to cause a decrease in the proportion of positive correlation, again most strongly for B-B, followed by F-F, and lastly by B-F.Overall, at the all-correlation level, we found no consistent support for the differences postulated for bacterial and fungal responses in H1. For example, strong increases in the proportion of positive correlations under drought could be found in all microbial pairings for some compartments (B-B in leaf and root, F-F in rhizosphere and soil, and B-F in root and rhizosphere) (Fig. 2a, Supplementary Figs. 2, 3). Neither did we find consistent support for the differences ascribed to bacteria and fungi in H2 as the strongest decreases in the proportion of positive correlations during rewetting occurred at F-F in rhizosphere and soil, and B-B in leaf and root (Fig. 2b, Supplementary Figs. 2, 3).Fig. 2: Correlations of microbes in drought stress and drought relief.Estimates of combined correlations (row a) show an increase in positive correlations under drought stress across the four compartments (root, black; rhizosphere, blue; soil, red; leaf, green). Data points underlying the lines in the figure are provided in the alternative version in Supplementary Fig. 2. This result is in line with the stress gradient hypothesis which posits that stressful environments favor positive associations because competition will be less intense than in benign environments32,33,36,37. Note that positive trends in combined correlations can arise in two ways. First, from an increase of positive correlations (row b) that exceeds the rise in negative correlations (row c), e.g., Leaf bacterial-bacterial (Bac-Bac) correlations or rhizosphere fungal-fungal (Fun-Fun) correlations in the drought period (Negative correlations in row C values are multiplied by −1 to facilitate comparison). Second, from a decrease in negative correlations that exceeds a decrease in positive correlations, e.g., root bacterial-bacterial correlations or root bacterial-fungal (Bac-Fun) correlations in drought. Combined (a), positive (b) and negative (c) estimates of correlation (Spearman’s rho, ρ) are given for four compartments (root, rhizosphere, soil and leaf), and three types of correlations (Bacterium-Bacterium, Fungus-Fungus, Bacterium-Fungus). T-tests (two sided) were carried out for linear mixed effect modelling that incorporates link type and compartments as random factors. Detailed distribution densities of correlations are presented in Supplementary Fig. 3. Source data are provided as a Source Data file.Full size imageWe found support for the stress gradient hypothesis because drought increased the relative frequency of positive correlations among microbial taxa (Fig. 2a, Supplementary Figs. 2, 3). The increases were due, largely, to B-B correlations in leaf and F-F correlations in the rhizosphere during drought, when the relative frequency of positive correlations was increased (Fig. 2b, Supplementary Figs. 2, 3) and the frequencies of negative correlations were decreased or weakly increased (Fig. 2c, Supplementary Figs. 2, 3). Less obvious increases in the relative frequency of positive correlations (such as B-B in root, F-F in soil, and B-F in root and rhizosphere) occurred where drought reduced both positive and negative correlations, but the losses of negative correlations exceeded those of positive correlations (Fig. 2, Supplementary Figs. 2, 3).In support of the expectation that correlations would be more negative between taxonomic groups than within taxonomic groups, we found that the relative frequency of positive correlations was generally lower for B-F than B-B and F-F correlations (Fig. 2, Supplementary Figs. 2, 3). Moreover, as ecological modeling has indicated that negative associations should promote stability of communities40, we hypothesize that B-F correlations would be more stable than B-B and F-F networks in response to drought stress. However, we found no support for this hypothesis, as B-F correlations (for example in root) did not always show the least response to drought stress (Fig. 2, Supplementary Figs. 2, 3).Testing H1 and H2 at species co-occurrence levelFor our final test of H1 (resistance) and H2 (resilience) we focused on co-occurrence networks based on significant, positive correlations. These networks have been reported to be destabilized for bacteria but not for fungi in mesocosms subject to drought stress19, and shown to be disrupted for bacteria in natural vegetation studied over gradients of increasing aridity41,42. Using these results as guides, for H1 we expected that drought stress should disrupt co-occurrence networks most strongly for B-B, followed by F-F, and lastly by B-F. For H2 we expected that relief of stress by rewetting should strengthen microbial co-occurrence networks most strongly for B-B, followed by F-F, and lastly by B-F.For this test we constructed microbial co-occurrence networks using significant positive pairwise correlations between microbial taxa, B-B, F-F and B-F, and compared the network complexity between fully irrigated control and drought, and between control and rewetting following drought. In general, we found no consistent support for the difference between bacteria and fungi inherent in H1. Rhizosphere was the one compartment where B-B vertices dropped and F-F vertices rose in response to drought, as expected, but this result was offset in root and soil, where vertices dropped in all networks, B-B, F-F and B-F (Figs. 3, 4; Supplementary Figs. 4, 5). Analysis by co-occurrence networks highlighted the differences between plant compartments. In root drought strongly disrupted networks of B-B, B-F and F-F, but in the other three compartments, network disruption was weaker, and networks were even enhanced by drought for F-F in rhizosphere and B-B in leaf (Figs. 3, 4).Fig. 3: Networks of significant positive cross-taxonomic group correlations (bacteria and fungi).a Fungal operational taxonomic units (OTUs) (blue) and bacterial OTUs (black) are graphed as nodes. Significant positive Spearman correlations are graphed as edges (ρ  > 0.6, false discovery rate adjusted P  More

Population genetic data collection from primary data sourcesFigure 4 describes the overall data collection workflow for the four datasets that comprise CaliPopGen. We first identified literature potentially containing population genetic data for California by querying the Web of Science Core Collection (https://webofknowledge.com/) for relevant literature from 1900 to 2020 with the terms: topic = (California*) AND topic = (genetic* OR genomic*) AND topic = (species OR taxa* OR population*). We included only empirical peer-reviewed literature and excluded unreviewed preprints. In using these search terms, our goal was to broadly identify genetic papers focused on California with population or species-level analyses, while avoiding purely phylogenetic studies or those focused on agricultural or model species. This resulted in 4,942 unique records.Fig. 4Flow chart of the data collection process that generated the CaliPopGen databases.Full size imageWe next screened titles and abstracts to retain articles that: (1) provided data on populations of species which are self-sustaining without anthropogenic involvement; (2) included at least some eukaryote species; (3) included population(s) sampled within California; (4) mentioned measures of genetic diversity or differentiation; and (5) were not reviews (thus restricting our search to only primary literature). We retained 1869 studies after this first pass of literature screening (see Technical Validation for estimate of inter- and intra-screener bias).Our second, more in-depth screening pass involved reading the full text of these 1869 studies. We had two goals. First, we confirmed that retained papers fully met all five of our inclusion criteria (the first screen was very liberal with respect to these criteria, and many papers failed to meet at least one criterion after close reading). Second, we eliminated papers where the data were not presented in a way that allowed us to extract population-level information. For example, many of the more systematics-focused studies pooled samples from large, somewhat ill-defined regions (“Sierra Nevada” or “Southern California”); if such regions were larger than 50 km in a linear dimension, we deemed them unusable for making geographically-informative inferences. Other studies presented summaries of population data, often in the form of phylogenetic networks or trees, but did not include information on actual population genetic parameters and therefore were not relevant to our database. We retained 528 publications after this second pass.From this set of papers, we extracted species, locality, and genetic data for each California population or sampling locality described in each study (Fig. 3A). This included Latin binomial/trinomial, English common name, population identifiers, and geographic coordinates of sampling sites. We also noted population/sampling localities that were interpreted as comprised of interspecific hybrids, and listed both parental species. We collected population genetic diversity and differentiation statistics for each unique genetic marker for each population/sampling locality; as a result, a sampling locality may have multiple entry rows, one for each locus or marker type. Parameters extracted for each population/marker combination include sample size, genetic marker type, gene targets, number of loci, years of sampling, and reported values for effective population size (Ne), expected (HE) and observed (HO,) heterozygosity, nucleotide diversity (π, pi), alleles-per-locus (APL), allelic richness (AR), percent polymorphic loci (PPL), haplotype diversity (HDIV), inbreeding coefficient (e.g. FIS, FIT, GIS), and pairwise population genetic comparison parameters (FST, GST, DST, Nei’s D, Jost’s D, or phi). We note that while there are technical differences between allelic richness and alleles-per-locus, source literature often used the terms interchangeably, and we include the parameters and their values as named in the source. We define marker type as the general category of genetic marker used (e.g., “microsatellite” or “nuclear”), while gene targets are the specific locus/loci (e.g., “COI”). We present these data in two separate datasets, one containing all population-level genetic summary statistics (Dataset 121, see Fig. 3C and detailed description in Table 1) and a second for estimates of pairwise genetic differentiation (Dataset 221, see Fig. 3D and detailed description in Table 2).Table 1 Description of the population genetic data in Dataset 121.Full size tableTable 2 Description of the pairwise genetic distance data in Dataset 221.Full size tableAll genetic data were extracted directly from the source literature. However, we also updated or added to the metadata for these population genetic values in several ways. We included kingdom, phylum, and a lower-level taxonomic grouping for each species (usually class), and updated scientific and common names based on the currently accepted taxonomy of the Global Biodiversity Information Facility22. When geographic coordinates were not provided for a sampling locality, as was frequently the case in the older literature, we used Google Maps (https://www.google.com/maps) to georeference localities based on either in-text descriptions or embedded figure maps guided by permanent landmarks like a bend in a river or administrative boundaries. Because this can only yield approximate coordinates, we recorded estimated accuracy as the radius of our best estimate of possible error in kilometers. If coordinates were provided in degree/minute/seconds, we used Google Maps to translate them to decimal degrees. In cases where coordinates were not provided and locality descriptions were too vague to determine coordinates with less than 50 km estimated coordinate error, we did not attempt to extract coordinates but still provide the genetic data. All coordinates are provided in the web Mercator projection (EPSG:3857). We excluded studies that reported genetic parameter values only for samples aggregated regionally (“Southern California” or “Sierra Nevada”). If marker type was not explicitly included, we classified marker type based on the gene targets reported, if provided.Life history trait data collectionTo increase the utility of CaliPopGen, we also assembled data on life history traits for all animal (Dataset 321) and plant (Dataset 421) species contained in Datasets 121 and 221. We assembled trait data that have previously been shown to correlate with genetic diversity, including those related to reproduction, life cycle, and body size, as well as conservation status (e.g.23,24,25,26,). Life history data were compiled by first referencing large online repositories, often specific to taxonomic groups, like the TRY plant trait database27, and the Royal Botanic Gardens Kew Seed Information Database28. If trait data for species of interest were unavailable from these compilations, we conducted keyword literature searches for each combination of species and life history trait, and extracted data from the primary literature. When data were not available for the subspecies or species for which we had genetic data, we report values for the next closest taxonomic level, up to and including family, as available in the literature.For both animals and plants, we defined habitat types as marine, freshwater, diadromous, amphibious, or terrestrial. Marine species include those that are found in brackish or wetland-marine habitats, as well as bird species that primarily reside in marine habitats. Freshwater species include those that are found in wetland-freshwater habitats, as well as species that primarily reside in freshwater. The diadromous category includes fish species that are catadromous or anadromous. We considered species to be amphibious if they have an obligatory aquatic stage in their life cycle, but also spend a significant portion of their life cycle on land. Terrestrial species were defined as those that spend most of their life cycle on land and are not aquatic for any portion of their life cycle. In a few cases (e.g., waterbirds that are both freshwater and marine, semi-aquatic reptiles), a species could reasonably be placed in more than one category, and we did our best to identify the primary life history category for such taxa. If the taxonomic identity of an entry was hybrid between species or subspecies, this was noted in the speciesID column and no life history data were reported.The CaliPopGen Animal Life History Traits Dataset 321 (description of dataset in Table 3) includes habitat type, lifespan, fecundity, lifetime reproductive success, age at sexual maturity, number of breeding events per year, mode of reproduction, adult length and mass, California native status, listing status under the US Endangered Species Act (ESA), listing status under the California Endangered Species Act (CESA), and status as a California Species of Special Concern (SSC). For some traits, value ranges were recorded–for example, minimum to maximum lifespan. In other cases, we recorded single values and, when available, a definition of this single value, (for example, minimum, average, or maximum lifespan). We report either the range of the age of sexual maturity (minimum to maximum), or a single value, depending on the available literature. For sexually dimorphic species, we report female adult length and weight when available, because female body size often correlates with fecundity. Across animal taxonomic groups, different measures of body size and length measurements are often used, reflecting community consensus on how to measure size. Given this variation, we report the type of length measurement, if available, as Standard Length (SL), Fork Length (FL), Total Length (TL), Snout-to-Vent Length (SVL), Straight-Line Carapace (SLC), or Wingspan (WS).Table 3 Description of the animal life-history data in Dataset 321.Full size tableThe CaliPopGen Plant Life History Traits Dataset 421 (description of dataset in Table 4) includes habitat type, lifespan, life cycle, adult height, self-compatibility, monoecious or dioecious, mode of reproduction, pollination and seed dispersal modes, mass per seed, California native status, NatureServe29 element ranks (global and state ranks, see Table 5 for definitions), listing status under the Federal Endangered Species Act (ESA), and listing status under the California Endangered Species Act (CESA). In contrast to most animal species, plant lifespan was typically reported as a single value. We define life cycles as the following: Annual: completes full life cycle in one year; Biennial: completes full life cycle in two years; Perennial: completes full life cycle in more than two years; Perennial-Evergreen: perennial and retains functional leaves throughout the year; Perennial-Deciduous: perennial and loses all leaves synchronously for part of the year. Some species are variable (for example, have annual and biennial individuals), and in those cases we attempted to characterize the most common modality.Table 4 Description of the plant life-history data in Dataset 421.Full size tableTable 5 Description of the Conservation status (Heritage Rank) from California Natural Diversity Database29.Full size tableBecause of the paucity of data available for chromists and fungi, we did not extract life history trait data for the relatively few species in these taxonomic groups.Data visualization and summaryWe used the R-package raster (v3.1–5) to visualize the spatial extent of the data in CaliPopGen in Fig. 3. Panel (A) shows a summary plot of all unique populations of both the Population Genetic Diversity in Dataset 121 and the Pairwise Population Differentiation in Dataset 221. Panel (B) shows the total number of unique populations in each California terrestrial ecoregion. Panel (C) depicts all data entries of Population Genetic Diversity Dataset 121, summed for each 20×20 km grid cell. Panel (D) shows the density of pairwise straight lines drawn between pairs of localities in the Pairwise Population Differentiation Dataset 221, depicted as the total number of lines per 20×20 km grid cell. The number of populations and species of both Datasets 121 & 221 are summarized for each marine and terrestrial ecoregion in Table 6.Table 6 Summary of total numbers of populations and species per California ecoregion, separately for population genetic and pairwise datasets.Full size table More

Systematic reviewA systematic review (PROSPERO Registered: CRD42020165355) of all cluster-randomised control trials currently published on ITNs [including conventional nets (CTNs), pyrethroid-only long-lasting nets (pyrethroid-nets), and pyrethroid-piperonyl butoxide synergist nets (pyrethroid-PBO ITNs)], IRS or a combination of both interventions was completed to validate an established transmission model for Plasmodium falciparum malaria parameterised using entomological assessment of the interventions. Three search platforms, Web of Knowledge, PubMed and Google Scholar were used and further studies were included from three recent Cochrane reviews that have focused on individual- or cluster- randomised control trials testing either ITNs, IRS or both26,27,28. Our search criteria focused on studies within Africa, and those reporting an epidemiological outcome such as parasite prevalence or clinical incidence in a defined age-cohort. A total of 138 studies were initially identified for further assessment (Supplementary Fig. S2).Those papers identified through the systematic review went through another round of screening to ensure they fell within the scope of the work and were compatible with existing modelling parameterisation. These criteria included (i) the intervention falls within an existing World Health Organization recommendation (so trials, or arms of trials, investigating pyrethroid-pyriproxyfen ITNs29 or insecticide-treated curtains30 were excluded), (ii) the entomological impact of the product had been previously statistically characterised as part of the modelling framework (trials investigating DDT31 or propoxur IRS32 were excluded), (iii) the study was within the Africa continent, (iv) the study randomised interventions in the intervention arm across the community (i.e., interventions were not targeted to individuals or risk groups within the community)33,34,35, and (v) the study was not reporting a cluster-randomised design36. A full description of why studies and arms were excluded is provided in Data S1.1.RCTs can assess the public health impact of interventions using different epidemiological endpoints. The two most common metrics used in malaria RCTs is infection prevalence (generally assessing parasitemia in a particular age group using microscopy or rapid diagnostic tests) or clinical incidence (typically assessed using active case detection in a cohort, which had previously been cleared of infection). These metrics are both equally valid though may give different results. For example, it may be harder to change malaria parasite prevalence with a partially effective intervention in a high-transmission setting (where people have a high chance of being reinfected) compared to a low-transmission setting (where reinfection is less common). Similarly, estimates of clinical incidence will vary depending on the study design and regularity of follow-up. For example, there are practical constraints on the number of times people within an active cohort can be tested. In areas of higher transmission incidence estimates will be greater the more regularly the cohort is tested as people infected multiple times between screening will be less common. This information on the regularity of screening is not always reported making it difficult to adjust models accordingly. It is also important to account for cluster-level effects when interpreting trial results, and this cluster-level data is also mostly unavailable37. The systematic review identified more studies that evaluated interventions in their ability to change malaria prevalence, with 13 out of 14 RCTs showing how the intervention changed parasite prevalence between the study arms compared with 8 RCTs, which reported changes in clinical incidence. Therefore, we focus on prevalence as our metric for epidemiology impact in this framework though note this should be repeated with clinical incidence estimates should more data become available. The final dataset had 73 cross-sectional surveys of prevalence in a defined age-cohort, 37 trial arms from 13 different RCTs.Characterising the entomological impact of ITNs and IRSExperimental hut trials (EHTs) measure the outcome of wild, free-flying, mosquito attempting to feed on volunteers resting indoors in the presence of an indoor intervention38. This includes (i) whether or not a mosquito is deterred away from a hut, which has the intervention (calculated by the number of mosquitoes found in the control hut relative to the intervention hut), (ii) whether the mosquito exits without feeding (repellence, measured as the percentage of alive unfed mosquitoes inside the intervention hut), (iii) the percentage entering the hut that successfully blood-feed, or (iv) the percentage of mosquitoes which die. Intervention efficacy is typically summarised for the intervention huts relative to a no-intervention (or untreated net) control huts, be it induced mortality (the increase in the percentage of mosquitoes dying over a 24-h period) or blood-feeding inhibition (the reduction in the percentage of mosquitoes receiving a blood-meal).EHTs use specially built structures that follow a defined floor-plan and set of specifications. There are multiple designs of experimental hut as they were originally intended to replicate the predominant type of housing found in the local area. We recently conducted a systematic review to capture the average behaviours of mosquitoes across different hut designs19. The two most used huts in Africa are the West African design and East Africa hut39 (a third hut—the Ifakara hut—is not considered here39). The meta-analyses showed that the associations describing the probable outcome of a mosquito feeding attempt (deterrence, repellence, successful feeding, or death) varies according to hut design. It is unclear that hut design best predicts epidemiological impact.Meta-analyses of EHT data have shown how the entomological efficacy of pyrethroid-nets has diminished over time, probably due to the rise of pyrethroid-resistant mosquitoes16,19,40, though there may be some manufacturing changes41. EHTs are conducted throughout Africa but are limited to the sites where the huts are built and cannot directly inform estimates of ITN efficacy outside of these areas. The most widely used quantitative measure for approximating the phenotypic level of resistance in the local mosquito population is the discriminating-dose bioassay. There are two main types of discriminating assays, the WHO susceptibility bioassay and the CDC bottle bioassay42,43. Both these assays measure the proportion of local Anopheline mosquitoes that survive 24-h following exposure to a discriminatory dose of pyrethroid for 60 min. Results from these bioassays are highly variable44 though collating data from multiple tests has shown clear trends over time45. The relationship between the level of resistance in the local mosquito population (as measured in a discriminating-dose bioassay) and the mortality induced by ITNs in EHTs can be used to extrapolate the results from hut trials to other geographical regions16.Modelling rationaleThe two main metrics recorded in EHTs do not capture all entomological impacts of ITNs and IRS. Though useful, induced mortality does not consider the sub-lethal impact of interventions whilst blood-feeding inhibition fails to differentiate between preventing blood-meals and killing mosquitoes, which are likely to have very different epidemiological impacts. Killing mosquitoes reduces the force of infection for users and non-users (through a community effect) so the overall effectiveness of treated nets and IRS will vary according to how abundantly and regularly they are used by the local human population. In addition, the impact of ITNs and IRS is likely to vary between sites because of factors such as the disease endemicity itself driven by societal behaviours, seasonality of transmission and the use of other malaria control interventions, amongst others. This means that raw EHT data is unlikely to directly correlate with the results of RCTs.EHTs are widely used to parameterise malaria transmission dynamics mathematical models46,47,48. These models rigorously quantify the outcome of each mosquito feeding attempt and, by making a limited number of assumptions, can estimate an overall entomological efficacy by combining the impact of the level of personal protection elicited by the intervention to the user and the indirect community effect provided to both users and non-users. Transmission dynamics mathematical models are designed to mechanistically capture the underlying processes governing malaria transmission and so can account for known non-linear processes such as the acquisition of human immunity49,50,51. This enables these models to translate the entomological efficacy quantified in an EHT into predictions of epidemiological impact given the characteristics of the site. Unfortunately, to date, there are no published EHTs that have been conducted alongside RCT evaluation of ITNs or IRS products (and therefore evaluated against the same mosquito population). To overcome this issue we parameterise the models using a meta-analyses of 136 EHT results16,19 collated from across Africa, which quantifies how mosquito deterrence, repellence, successful feeding, or death varies with time since the intervention is deployed and according to the level of pyrethroid resistance in the local mosquito population (as measured by the discriminating-dose bioassay). This approach has been able to recreate the epidemiological impact observed in RCTs evaluating a small number of ITNs15 or IRS products9, but this is the first attempt at using this method to validate the modelling framework against all trials evaluating nets and IRS.There is considerable uncertainty in how the entomological efficacy of treated ITNs varies with the level of resistance in the local population. This is a key relationship determining how field discriminating-dose bioassay data should be interpreted yet it is highly uncertain, with a recent meta-analyses indicating that it is equally well explained by two different functional forms (the logistic or log-logistic functions)19. Similarly, it is unclear whether the epidemiological impact of ITNs or IRS is best captured by all experimental hut data combined (Supplementary Fig. S14C, D)19 or if the meta-analyses should be restricted to just West or East African hut design data alone. To rigorously differentiate between these options six different models are run for each trial arm (n = 37), varying both the relationship between discriminating-dose bioassay and EHT mosquito mortality (either the logistic or log-logistic function) and the data used in the EHT meta-analyses (all data, East or West African design huts). The ability of these models to recreate the observed results is statistically compared and the most accurate selected for the main analyses.Transmission dynamics modelThe malaria transmission model that we use here incorporates the transmission dynamics of Plasmodium falciparum between human hosts and Anopheles mosquito vectors. The differential equations and associated assumptions of the original transmission model52 have been comprehensively reported in the Supplementary Material from Griffin et al.53, Walker et al.54 and Winskill et al.55. The model has been extensively fitted to data on the relationship between vector density, entomological inoculation rate, parasite prevalence, uncomplicated malaria, severe disease and death49,52,53,56,57. Model equations and assumptions are provided in the Supplementary Methods and https://github.com/jamiegriffin/Malaria_simulation. Unless stated (Supplementary Data S1), default parameters are taken from these papers.Data requirements for model simulationThe transmission model can be parameterised to describe the specific ecology of each RCT location using data on the mosquito bionomics, seasonal transmission patterns, historic use of various interventions—principally insecticide-treated ITNs or the residual spraying of insecticides (IRS)—and baseline endemicity. These data are recorded within the research articles reporting the trials at the trial arm level (Supplementary Data S1.2 notes where data are available and which resources were used; Supplementary Data S1.3 lists the key data identified for model parameterisation) and Supplementary Fig. S1 provides a diagram of how they are combined to inform the model.Briefly, the Anopheles mosquito species composition at baseline is used to determine the proportion of mosquitoes with bespoke behaviours that could alter exposure risk to mosquito bites and thus transmission risk. Species-specific mosquito behaviours are parameterised from systematic reviews on anthropophagy, using the human blood index47,58,59, and the proportion of mosquito bites that are received indoors or in bed because this impacts the efficacy estimate for indoor interventions60.Other information that are specific to each trial also help interpret our success at predicting, or not, the observed results of an intervention tested in an RCT; the diagnostic used to measure prevalence or incidence is useful because different tests have different sensitivities61, which can be included in the model framework54. The baseline burden of infection is particularly important to enable the model to be calibrated to the endemicity of the study site by varying the number of mosquitoes per person (the human:mosquito ratio). This is determined by a cross-sectional estimate of parasite prevalence in a defined age-cohort at a particular time of year of the baseline survey.For any location, the current level of endemicity is determined by the historic interventions already operating at the site. Therefore, wherever possible, ITN use and the historic use of sprayed insecticides, as well as the estimated proportion of clinical cases that are drug-treated, are included as baseline parameters.In addition to the waning potency of insecticide active ingredient outlined above, the impact of nets can also wane because of changes in the proportion of people using them. This can be driven by the quality of the product, seasonal patterns in humidity or other social patterns of use62,63,64. Where data are available, this waning adherence to net use is captured by fitting an exponential decay function to the proportion of people using nets measured at cross-sectional surveys throughout the trials:$${{{{{{{mathrm{U}}}}}}{{{{{mathrm{sage}}}}}}}}_{i}={e}^{-{sigma }_{i}t}$$
(1)
where σ is a parameter determining how rapidly people stop using nets in an intervention arm i of the trial and t is time in years. Parameter estimates for pyrethroid-only and pyrethroid-PBO ITNs are provided for different levels of resistance for the 6 potential methods of associating bioassays and using data (Supplementary Data S1.4).The IRS product used is equally important as the entomological impact of different products vary, particularly for pyrethroid-based IRS in the presence of resistant mosquitoes9. Supplementary Data S1.5 show the parameter estimates for products included in the analysis.The seasonality of transmission has been defined previously for each RCT site (at the administration subunit 1 level) using normalised rainfall patterns obtained from the US Climate Prediction Center65. The daily time series are aggregated to 64 points per year for years 2002 to 2009. A Fourier function is fitted to these data to capture seasonality by reconstructing annual rainfall patterns54,66. We deliberately do not match rainfall data from the respective RCTs, which would likely improve the model estimates because we are ultimately testing whether this framework has predictive power across future years or alternative ecologies, where we will not know how rainfall will exactly impact mosquito densities and hence malaria transmission.Statistical analysisThe mean simulated malaria prevalence (matching the age-cohort of the trial) is recorded for all RCT surveys timepoints. This equates to a total of 73 cross-sectional surveys post-implementation. The process was repeated using the 6 different entomological parameter sets (the relationship between bioassay and hut trial mortality and the hut design used to summarise treated net entomological impact). An illustration of the different models and their fit to data is demonstrated in Supplementary Fig. S17 for a recent study trialling pyrethroid-only nets, pyrethroid-PBO ITNs alone or in combination with a long-lasting IRS product in Tanzania5. The difference between the observed and predicted prevalence at each timepoint is shown for all RCTs in Supplementary Fig. S18. A simple linear regression is conducted comparing observed and predicted results are summarised in Supplementary Table 3. Let Xi denote the malaria prevalence predicted by the model at timepoint i while Yi is the observed prevalence. The regression,$${Y}_{i}=m{X}_{i}$$
(2)
for i = 1,…,c + n, where m is the gradient between the observed and predicted result (consistent across studies), c is the number of post-intervention datapoints in the control arms and n is the number of post-intervention datapoints in the intervention arms (c + n = 73 for analyses of all RCTs). Better fitting models have a higher adjusted R2 (adjusted R2 values of one indicate the model is perfectly predicting the trial result) whilst the gradient of the regression m indicates any bias (with value of one reporting the model can predict prevalence equally well across the endemicity range). Results are presented for all ITNs and IRS RCTs and separately for RCTs of different types of (pyrethroid-only ITNs, pyrethroid-PBO ITNs and IRS, Supplementary Table 3). The log-logistic model (results 4–6 in Supplementary Table 3) describing the relationship between bioassay and hut trial mortality consistently fits the data better, with models fit using either all hut trial data or East African design huts having a similar accuracy (adjusted R2 = 0.95). This parameter combination also had the least bias, with the best fit regression line being closer to one.The average efficacy of the different ITNs and IRS combinations was calculated by comparing malaria prevalence for the different trial arms to the respective control arms at matched timepoints following the introduction of interventions. Let ({E}_{{jk}}^{l}) be the relative reduction in the malaria prevalence between the control (k = 0) to intervention (k = 1) arms at matched timepoint j in the same trial for either the predicted (l = Xjk) or observed (l = Yjk) malaria prevalence,$${E}_{j}^{X}=({{X}_{j0}-{X}}_{j1})/{X}_{j0},{{{{{rm{ and }}}}}},{E}_{j}^{Y}=({Y}_{j0}-{Y}_{j1})/{Y}_{j0}$$
(3)
for j = 1,…,n. The goodness of fit for the efficacy estimates is calculated in a similar manner to the prevalence estimates by substituting in ({E}_{j}^{X}) and ({E}_{j}^{Y}) into Xi and Yi in E2, respectively. Models are on average able to estimate the efficacy of the interventions at different timepoints (Supplementary Table 3). Estimates for some timepoints diverge substantially (for example, the study testing conventional nets in the Gambia relative to untreated nets67 measured negative effect in one setting; the treated net arm having more infected children whereas the model predicted a 12.5% reduction due to the CTN (with parameters derived from all EHT data and the log-logistic function, 4 in Supplementary Table 3), Supplementary Data S1.8), but in most studies the trial average (averaged across all timepoints) is remarkably consistent. Accuracy is lower than estimates of absolute prevalence, in part because the difference between the percentage of people slide positive in low-endemicity settings may be relatively modest in absolute terms but might represent a substantial difference as a percentage. It is also important to note that when the models do systematically miss some timepoints, this is consistent across the control and treated arms. For example, in the Protopopoff et al. study in Tanzania5 (Figs. S14 and S17) efficacy is over-estimated in all arms 18 months after the start of the trial, but the relative difference between the arms (in terms of ordering, and the efficacy estimate) is relatively consistent. This indicates that unmeasured factors, such as differences in the timing and duration of the rainy season, may have occurred across all trial arms. As previously, the log-logistic functional form describing the relationship between bioassay and hut trial mortality consistently fits the data better (Supplementary Table 3, options 4 to 6). The models fit describing the entomological efficacy of any net using all EHT data predicts efficacy data better with East African design hut data providing similar accuracy (adjusted R2 = 0.64 vs. 0.62, respectively). Following this we select the log-logistic functional form to describe the relationship between mortality in the discriminating-dose bioassay and EHT and characterise the entomological efficacy of treated ITNs using data from both East and West African design huts for the main analyses (Fig. 2B, C).The ability of the best-performing model (Supplementary Table 3, column 4: log-logistic function and all EHT data) to capture the relative drop in prevalence over time compared to the baseline (pre-intervention) estimate is shown in Supplementary Fig. S19. This value is denoted as ({dot{E}}_{t}^{l}) and is calculated as,$${dot{E}}_{t}^{X}=({X}_{0}-{X}_{t}),{{{{{rm{and}}}}}},{dot{E}}_{t}^{Y}=({X}_{0}-{Y}_{t})$$
(4)
where ({X}_{0}) is the malaria prevalence at baseline (prior to intervention deployment with the exception of Chaccour et al.68) observed from the RCT and the model is calibrated to this endemicity. Xt is then the subsequent cross-sectional survey observed for each study, and RCTs have different numbers of surveys ranging from 1 to 4 in the published literature. The corresponding model estimate is represented by Yt. Estimates are calculated for all post-intervention timepoints in both control and intervention arms and are shown in Fig. S19A. The difference between ({dot{E}}_{t}^{X}) and ({dot{E}}_{t}^{Y}) can be used to explore how closely the model is able to predict this absolute difference observed in the trials (a value of 0 indicates exact match, high predictive ability). The model overestimates the performance of IRS only, deployed in 1995 using the pyrethroid IRS ICON CS 10% (Syngenta), but otherwise there is no difference in the models’ ability to estimate different ITN interventions or combination net and IRS interventions, be it the absence of an intervention, conventional dipped-nets, pyrethroid-only nets, pyrethroid-PBO ITNs with or without IRS (Fig. S19B). All code is available69.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Sabkha soil analysisSoil biochemical analysis shows higher concentration of salt in the middle (average of 200 ppt) which decreases to 20 ppt outside of the sabkha (Fig. 1B). Soil conductivity and total dissolved solids (TDS) were positively correlated with salinity (r = 0.8, p  More

Virtual hydrostatic model parametersVarious morphological characteristics were held constant in order to isolate and manipulate the variable of conch shape. A CT-scanned Nautilus pompilius conch was essentially morphed into ammonoid-like conch shapes, populating the Westermann morphospace22 while holding constant septal morphology, septal spacing, and shell/septal thicknesses (Fig. 9). Furthermore, body chamber proportions were determined by iteratively computing soft body volumes that yield Nautilus-like chamber liquid (~ 12% of the phragmocone volume retained)67,68. Septal spacing was measured as the angle from the ventral attachment of the current and previous septa, and the spiraling axis of the conch. Because septal spacing differs in early ontogeny (Fig. S11), only measurements from the 7th to 33rd (terminal) septum were considered. The average angle of 23.46° ± 3.32° (standard deviation) was rounded to 23° and held constant throughout the ontogeny of the hydrostatic models.Figure 9Hydrostatic models of theoretical planispiral cephalopods. These models were constructed by morphing a Nautilus pompilius conch into ammonoid shapes (see “Methods”): (a) oxycone, (b) serpenticone, (c) sphaerocone, and (d) morphospace center. The centers of buoyancy and mass are denoted by the tips of the blue (upper) and red (lower) cones. Prime symbols (′) refer to transparent, transverse views of each respective conch shape. (e) Westermann morphospace22 showing relative positions of these conch shapes. All models were rendered in MeshLab76.Full size imageShell and septal thicknesses were measured with digital calipers from a physical specimen of Nautilus pompilius (Table S13). These measurements were recorded as a ratio of inner whorl height (measured from the ventral point on the current whorl to the ventral point on the previous whorl). These ratios were used in the theoretical models to define shell and septum thicknesses (3.1% of inner whorl height for shell thickness and 2.1% of inner whorl height for septal thickness; Table S13).Hydrostatic model constructionThe near-endmember models were constructed from representative ammonoid specimens (Sphenodiscus lobatus and S. lenticularis—oxycone; Dactylioceras commune—serpenticone; Goniatites crenistria—sphaerocone). Lateral and transverse views were measured from figured specimens for the oxycone (Fig. 5 of Kennedy et al.69), serpenticone (Fig. 2 of Kutygin and Knyazev70), and sphaerocone (Figs. 17 and 20 of Korn and Ebbighausen71). These models were constructed with array algorithms similar to earlier hydrostatic models9,35,72, which were used in a piecewise manner to account for allometric changes in coiling throughout ontogeny (Table S14). These arrays replicated the adult whorl section backwards and translated, rotated, and scaled each successive one. These whorl sections were bridged together to create a single tessellated surface representing the outer interface of the shell. Shell thickness was defined by shrinking the original whorl section so that the thickness between the two was equal to 3.1% of the inner whorl height (Table S13), then using the same array to build the internal interface of the shell. The morphospace center was constructed from previously used conch measurements18 and averaging the whorl section shape in blender (Fig. S12). The corresponding Westermann morphospace parameters (Fig. S13) for each morphology are reported in Table S15.Virtual models of the septa were derived from the CT-scan of Nautilus pompilius (Fig. S14). A single septum was isolated from the adult portion of the phragmocone then smoothed to delete the siphuncular foramen. This septum was placed within the whorl section of each theoretical model and stretched in the lateral directions until it approximately fit. The “magnetize” tool in Meshmixer (Autodesk Inc.) was used to attach the septal margin to the new whorl section so that the Nautilus suture was transferred to the new whorl section. The septum was then smoothed to reconcile the first order curves with the new location of the septal margin. The respective septum for each theoretical model was then replicated with the same array instructions used to build the shell. Because each replicated object was rotated one degree (Table S14), 22 septa were deleted in between every two so that the septal spacing was equal to 23° (Fig. S11).For each theoretical model, the septa were unified with the model of the shell using Boolean operations in Netfabb (Autodesk Inc.). To perform hydrostatic calculations, virtual models must be created for each material of unique density. The virtual model of the shell constrains the shape of the soft body (within the body chamber) and chamber volumes (within the phragmocone). These internal interfaces were isolated from the model of the shell, then their faces inverted for proper, outward-facing orientations of their normals. A conservative soft body estimate was created, aligning with previously published reconstructions64,65,73. The profile shape of this soft body was scaled and maintained between each model. External interfaces of the shell and soft body were also isolated to create a model of the water displaced by each theoretical cephalopod. Each of these models are necessary for hydrostatic calculations (buoyancy and the distribution of organismal mass).Each hydrostatic model is stored in an online repository (Dataset S1; https://doi.org/10.5281/zenodo.5684906). The hydrostatic centers of each virtual model and their volumes and masses are listed in Tables S16 and S17.Hydrostatic calculationsEach theoretical model was scaled to have equal volume (near one kilogram; 0.982 kg–a result of arbitrarily scaling the sphaerocone model to 15 cm in conch diameter). An object is neutrally buoyant when the sum of organismal mass is equal to the mass of water displaced (the principle of Archimedes). The percentage of chamber liquid can be computed to satisfy this condition.$${Phi } = frac{{left( {frac{{{text{V}}_{{{text{wd}}}} {uprho }_{{{text{wd}}}} – {text{V}}_{{{text{sb}}}} {uprho }_{{{text{sb}}}} – {text{V}}_{{{text{sh}}}} {uprho }_{{{text{sh}}}} }}{{{text{V}}_{{{text{ct}}}} }}} right) – left( {{uprho }_{{{text{cl}}}} } right)}}{{left( {{uprho }_{{{text{cg}}}} – {uprho }_{{{text{cl}}}} } right)}}$$
(1)
where Vwd and ρwd are the volume and density of the water displaced, Vsb and ρsb are the volume and density of the soft body, Vsh and ρsh are the volume and density of the shell, ρcl is the density of cameral liquid, ρcg is the density of cameral gas, and Vct is the total volume of all chambers. A soft body density of 1.049 g/cm3 is used based on bulk density calculations of Nautilus-like tissues74, a seawater-filled mantle cavity, and thin calcitic mouthparts21. A shell density of 2.54 g/cm374, cameral liquid density of 1.025 g/cm375, and cameral gas density of 0.001 g/cm3 are adopted from recent hydrostatic studies.Other hydrostatic properties depend on the relative positions of the centers of buoyancy and mass. The center of buoyancy is equal to the center of volume of water displaced. This center and the centers of each virtual model of unique density were computed in the program MeshLab76. The individual centers for each organismal model (soft body, shell, cameral liquid and cameral gas) were used to compute the total center of mass, with an average weighted by material density:$$M = frac{{sum left( {L*m_{o} } right)}}{{sum m_{o} }}$$
(2)
where M is the total center of mass in a principal direction, L is the center of mass of a single object measured with respect to an arbitrary datum in each principal direction, and (m_{o}) is the mass of each object with unique density. Equation 2 was used in the x, y, and z directions to compute the 3D coordinate position of the center of mass. The centers of mass for the chamber contents (liquid and gas) were set equal to the center of volume of all chambers, a minor assumption considering the capillary retention of liquid around the septal margins in the living animals62.The hydrostatic stability index (St) is computed from the relative location of the centers of buoyancy (B) and mass (M), normalized by the cube root of volume (V) for a dimensionless metric that is independent of scale:$$S_{t} = frac{{ sqrt {left( {B_{x} – M_{x} } right)^{2} + left( {B_{y} – M_{y} } right)^{2} + left( {B_{z} – M_{z} } right)^{2} } }}{{sqrt[3]{V}}}$$
(3)
where the subscripts correspond to the x, y, and z components of each hydrostatic center.Apertural orientations were measured in blender after orienting each model so that the center of buoyancy was vertically aligned above the center of mass. Apertural angles of 0° correspond to a horizontally facing soft body, while angles of + 90° and − 90° correspond to upward- and downward-facing orientations, respectively.Thrust angles were measured from the hyponome location (ventral edge of the aperture) to the midpoint of the hydrostatic centers, with respect to the horizontal. Thrust angles of 0° infer idealized horizontal backward transmission of energy into movement, while thrust angles of + 90° and − 90° infer more efficient transmission of energy into downward and upward vertical movement, respectively.Biomimetic robot constructionTo isolate the variable of shell shape on swimming capabilities, only the external shape, and static orientation of each virtual hydrostatic model were used to build physical, 3D printed robots. That is, each model has artificially high hydrostatic stability (Tables S3) to nullify the effect of the thrust angle (the angle at which thrust energy passes through the hydrostatic centers and most efficiently transmits energy into movement; Table S4). Less stable morphotypes (e.g., serpenticones and sphaerocones) are more sensitive to the constraints imposed by this hydrostatic property.Space constraints inside each model were determined by first constructing a propulsion system and electronic components that operate the motor. The models use impeller-based water pumps (Figs. 1d and 10a) driven by a brushed DC motor. This system creates a partial vacuum by centrifugal acceleration, drawing water from a “mantle cavity” and expelling it out of a “hyponome”. This system was optimized by iteratively designing models in Blender77, then testing 3D-printed, stand-alone water pumps. After three iterations, a four-blade impeller and gently tapering hyponome (inner diameter at distal end = 6.7 mm) were chosen. The electronic components used to drive the motor consist of an Arduino Pro Micro microcontroller, a motor driver, and two batteries (Fig. 10). A 3.7 V battery operates the microcontroller, and a larger 7.4 V battery supplies power to the motor. Communication is achieved via infrared, allowing specification of the jet pulse duration, number of pulses, and the power level of the motor (using pulse-width modulation; PWM). Each of these electronic components fold into a compact cartridge capable of being plugged into 3D-printed models of each investigated shell shape (Figs. 2 and 10). Each model was designed with brackets to hold the electronics cartridge in place. The sphaerocone had the most severe space constraints, with low conch diameter to volume ratio. After determining the space required for the electronics (Fig. 10) this model was scaled to 15 cm, and all other models were scaled to have similar volumes (with subtle volume differences due to minor differences in soft body shape compared to the hydrostatic models).Figure 10Biomimetic cephalopod robot components. (a) Ventral view of the sphaerocone biomimetic robot (before covering the pump and mantle cavities) with assembled electronics cartridge to the right. (b) View of electronic components that fit into the cartridge. (c) Electronics cartridge placed in robot. These two halves are fit together with wax to create a water-tight seal. Each model component is denoted by letters in circles: A = Arduino microcontroller, B = microcontroller charger / voltage regulator, C = motor driver, D = infrared sensor, E = indicator LED, F = microcontroller battery (3.7 V), G = motor battery (7.4 V), H = brushed motor, I = impeller and water pump cavity, J = electronics cartridge. The colors of annotations correspond to components depicted in Figs. 1 and 2.Full size imageIn addition to having a propulsion system, biomimetic cephalopod robots must also be capable of neutral buoyancy, while assuming the proper orientation in the water. These robots, and their once-living counterparts, each have differing material densities and associated mass distributions for each component. To reconcile these differences, the total mass and total centers of mass for each model were manipulated by controlling the volume and 3D distribution of the 3D-printed PETG (polyethylene terephthalate glycol) thermoplastic. That is, the shape of this material holds each model component in place while correcting for these differences in hydrostatics. The PETG mass required for neutral buoyancy was found by subtracting the mass of every other model component from the mass of the water displaced by the model (i.e., electronics cartridge, bismuth counterweight, liquid, motor, batteries, electronic components, and self-healing rubber; Table S1). This model configuration also allows buoyancy to be fine-tuned in water, compensating for potential density differences between the virtual water and the actual water in the experimental settings. That is, each virtual model accounts for ~ 9 g of internal liquid, but the actual volume of this liquid can be adjusted in the physical robot with a syringe through a self-healing rubber valve (Table S1; Fig. 1).The 3D position of the total center of mass was manipulated by accounting for the local centers of mass of each material of unique density. Materials like the batteries, motor, and electronic components were each assigned bulk density values because they are made up of composite materials. While this is an approximation, their contributions to the total center of mass are low because they account for small fractions of the total model mass (Tables S1 and S2). These components, like all others, were digitally modeled in Blender77 and their volumes and centers of mass were computed in the program MeshLab76. A dense, bismuth counterweight was also modeled, and positioned to artificially stabilize each model (pulling the z component of the total center of mass downward, while maintaining the horizontal components). The virtual model of this counterweight was used to make a 3D-printed mold, allowing a high heat silicone mold to be casted. The bismuth counterweight was cast from this silicone mold and filed/sanded to the dimensions of its virtual counterpart. Hyponomes were oriented horizontally, to yield movement in this direction. To maintain the same static orientation as the virtual model (same x and y center of mass components), the PETG center of mass was computed with the following equation:$$D_{PETG} = frac{{Mmathop sum nolimits_{i = 1}^{n} m_{i} – mathop sum nolimits_{i = 1}^{n} (D_{i} m_{i} )}}{{left( {m_{PETG} } right)}}$$
(4)
where DPETG is the location of the PETG center of mass from an arbitrary datum in each principal direction. M is the total center of mass in a particular principal direction, mi is the mass of each model component, Di is the local center of mass of each model component in a particular principal direction and mPETG is the mass of the PETG required for a neutrally buoyant condition. See Tables S1 and S2 for a list of model components and measurements.Each model was 3D printed with an Ultimaker S5 3D printer using clear (natural) PETG in separate parts, allowing the internal components to be implanted (i.e., brushed DC motors and bismuth counterweights). Each model part was chemically welded together with 100% dichloromethane, with minor amounts of cyanoacrylate glue used to fill seams (e.g., the water pump lid; Fig. 10a). Each final model consists of the main body (housing the water pump, motor, and counterweight), and a “lid” with brackets that house the electronics cartridge (Figs. 2 and 10). The main body and lid were fused together before each experiment by placing wax (paraffin-beeswax blend) along a tongue and groove seam, heating it with a hairdryer, then vigorously squeezing each part together. Surplus wax extruded from the seam was removed and smoothed, producing a water-tight seal.Thrust calibrationEven though each model was designed to have equal mantle cavity and pump cavity volumes, they produced slightly different thrusts. These differences were likely due to variable degrees of friction between the impellers and the surrounding water pumps. To correct for these differences, the thrust produced by each model was measured with a Vernier Dual-Range Force Sensor (0.01 N resolution). Each robot was attached at the hyponome location, through a series of pulleys, and to the sensor with fishing line (Fig. S1; similar to the methods used for living cephalopods78). Force was recorded for 30-s intervals at a sample rate of 0.05 s. During this time, each model was recorded jetting with a 6-s pulse for 15 trials (Fig. S2A). Each trial had initial noise from setting up the model, then peaked randomly when the fishing line became taught, then stabilized after some period of oscillation. Only the stabilized portion of the thrust profile was used to record thrust at 100% voltage for each model (Fig. S2B). The true zero datum was also subtracted from each of these trials. The lowest thrust from each of the models was used as a baseline (serpenticone and oxycone). Each model was recorded again for 15 trials by lowering the motor voltage in increments of 5% until they yielded similar thrusts (0.3 N) to the original serpenticone and oxycone trials (Fig. S2C). The final power levels were then determined for each model and adjusted with pulse-width modulation (PWM) through the microcontroller: serpenticone (100%), oxycone (100%), sphaerocone (95%), and morphospace center (85%).The peak thrust measured for 1 kg extant Nautilus is around 2 N16. The time-averaged thrust during each pulse is around 23% of this value (0.46 N16). This computed value slightly overpredicts observed maximum velocities for this animal (33 cm/s instead of 25 cm/s), so the appropriate time-averaged thrust is probably slightly lower. The motor in the robots quickly reaches its maximum thrust (~ 0.3 N) once initiated then quickly declines after shutting off (Fig. S2). Therefore, the thrust produced by the robots can be treated as a conservative Nautilus-like jet thrust close to the behavior of escape jetting. One-second pulse and refill intervals are also on par with values reported for extant Nautilus16.Robot buoyancyEach of the models were made near neutrally buoyant by adjusting the allotted ~ 9 g of internal liquid with a syringe through a self-healing rubber valve. The single-pulse experiments were performed in an external pool (ranging ~ 23.5 to 26.5 °C). The three-pulse and maneuverability experiments were performed in an internal pool (the Crimson Lagoon at the University of Utah). This internal pool had slightly higher temperatures (~ 28 °C), yielding lower ambient water densities than the virtual water. These conditions required slightly less internal liquid (~ 2–5 g). These differences in internal liquid masses produced negligibly small shifts in mass distributions because they are very small proportions of total robot masses (Table S1).Perfect neutral buoyancy cannot be practically achieved, but this condition can be closely approached. Each of the biomimetic robots experience subtle upward or downward movements of the course of their 5–15 s long trials due to slightly positive or negative buoyancies. Because these differences in buoyancy influence the vertical component of movement, only the horizontal components are considered for discussion. However, a comparison of velocities computed from full, 3D movement (Eq. 5) and restricted 2D components (Eq. 6) reveals that these differences are minor (Figs. S7 and S8). These comparisons demonstrate that model buoyancy did not substantially influence kinematics other than gross trajectories (Figs. 4 and S9).3D motion trackingAfter adjusting buoyancy, each model was positioned underwater with a grabber tool. This tool was fitted with a bundle of fiber-optic cable (Fig. S4) attached to an infrared remote control. Arduino code (Dataset S2) was uploaded to the microcontroller in the robot allowing jet pulse duration, number of pulses, and power to be adjusted with this remote control. After an infrared pulse is received, the motor activates, and activity is indicated by a green LED that illuminates the model from the inside. This light is used to determine time-zero for each trial of motion tracking.After sending an infrared signal, the movement of each model was recorded with a submersible camera rig fitted with two waterproof cameras (Fig. 3). Each of the four models were monitored during a single, one-second jet for at least 9 trials each. Additionally, the laterally compressed morphotypes (serpenticone and oxycone) were monitored during three, one-second pulses for 10 trials each. The inflated morphotypes (sphaerocone and morphospace center) were not able to be monitored over longer distances because they had the tendency to rotate about the vertical axis, obscuring views of the tracking points. In addition to horizontal movement, turning efficiency (maneuverability about the vertical axis) was monitored by directing the cameras with a top-down view of each model. A 90° elbow attachment for the hyponome was fit to each model to investigate the ease or difficulty of rotation. Each model was designed to spin counter-clockwise when viewed from above so that the influence of the motor’s angular momentum was consistent between models.Footage was recorded with two GoPro Hero 8 Black cameras at 4K resolution and 24 (23.975) frames per second, with linear fields of view. Motion tracking was performed with the software DLTdv879 to record the pixel locations of each tracking point (Figs. 1c and S4). These coordinates were transformed into 3D coordinates in meters using the program easyWand580. The tracking points on each model were used for wand calibration because the distances between these sets of points were fixed. Standard deviations of the reproduced tracking point distances of less than 1 cm were considered suitable.The 3D position datasets allowed velocity, acceleration, rocking, to be computed for each experiment. Additionally angular displacement and angular velocity was of interest for the rotation experiments about the vertical axis. Velocity was computed under two scenarios: (1) using the 3D movement direction between each timestep (Eq. 5), and (2) only considering the horizontal movement direction between each time step (Eq. 6). The latter scenario was preferred to nullify the influences of model buoyancies, which were not perfectly neutral and caused some degree of vertical movement.$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} + left( {z_{i} – z_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(5)
$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(6)
where V and t are velocity and time, and the subscripts i and i −1 refer to the current and previous time steps, respectively. Coordinate components are denoted by x, y, and z at each timestep. The averaged 3D location of both tracking points was used for each model (i.e., midpoints). Note that Eq. (5) uses the 3D form of the Theorem of Pythagoras, whereas Eq. (6) uses the 2D version. Time zero for each trial was defined as the frame where the robot was illuminated by the internal LED, indicating motor activity. Acceleration was modeled by fitting a linear equation to the datapoints during the one-second pulse interval(s) using the curve fitting toolbox in MATLAB R2020A.The artificially high hydrostatic stability of each model was designed to nullify rocking during movement. This behavior was computed for each model during the one-pulse experiments with the following equation:$$theta_{dv} = cos^{ – 1} left( {frac{{left( {z_{2} – z_{1} } right)}}{{sqrt {left( {x_{2} – x_{1} } right)^{2} + left( {y_{2} – y_{1} } right)^{2} + left( {z_{2} – z_{1} } right)^{2} } }}} right) – theta_{tp}$$
(7)
where (theta_{dv}) is the angle deviated from true vertical and (theta_{tp}) is the angle of the tracking points measured from the vertical in a static setting. The subscripts 1 and 2 of the x, y, and z coordinates refer to the anterior and posterior tracking points, respectively.Maneuverability about the vertical axis was determined by computing the angle between the horizontal components of each tracking point. The net angle from the starting angle for each trial was tabulated. Angular velocity was determined by dividing the change in angle between each frame by the frame duration (1/23.975 fps).Links to example motion tracking footage, and robotic models are deposited in an online repository60,61,63 (Dataset S2; https://doi.org/10.5281/zenodo.6180801). More

Chileshe, M. N. et al. Physico-chemical characteristics and heavy metal concentrations of copper mine wastes in Zambia: Implications for pollution risk and restoration. J. For. Res. https://doi.org/10.1007/s11676-019-00921-0 (2019).Article 

Google Scholar 
Sracek, O. Formation of secondary hematite and its role in attenuation of contaminants at mine tailings: Review and comparison of sites in Zambia and Namibia. Front. Environ. Sci. 2, 1–11 (2015).ADS 
Article 

Google Scholar 
Kayika, P., Siachoono, S., Kalinda, C. & Kwenye, J. An investigation of concentrations of copper, cobalt and cadmium minerals in soils and mango fruits growing on Konkola copper mine tailings dam in Chingola, Zambia. Arch. Sci. 1, 2–5 (2017).
Google Scholar 
Nazir, R. et al. Accumulation of heavy metals (Ni, Cu, Cd, Cr, Pb, Zn, Fe) in the soil, water and plants and analysis of physico-chemical parameters of soil and water collected from Tanda Dam Kohat. J. Pharm. Sci. Res. 7, 89–97 (2015).CAS 

Google Scholar 
Surbakti, E. P., Iswantari, A., Effendi, H. & Sulistiono. Distribution of dissolved heavy metals Hg, Pb, Cd, and As in Bojonegara Coastal Waters, Banten Bay. IOP Conf. Ser. Earth Environ. Sci. 744, 012085 (2021).Article 

Google Scholar 
Van Nguyen, T. et al. Arsenic and heavy metal contamination in soils under different land use in an estuary in northern Vietnam. Int. J. Environ. Res. Public Health 13, 1091 (2016).Article 
CAS 

Google Scholar 
Yabe, J. et al. Uptake of lead, cadmium, and other metals in the liver and kidneys of cattle near a lead-zinc mine in Kabwe, Zambia. Environ. Toxicol. Chem. 30, 1892–1897 (2011).CAS 
PubMed 
Article 

Google Scholar 
Salem, M. A., Bedade, D. K., Al-ethawi, L. & Al-waleed, S. M. Heliyon Assessment of physiochemical properties and concentration of heavy metals in agricultural soils fertilized with chemical fertilizers. Heliyon 6, e05224 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Tuakuila, J. et al. Worrying exposure to trace elements in the population of Kinshasa, Democratic Republic of Congo (DRC). Int. Arch. Occup. Environ. Health 85, 927–939 (2012).CAS 
PubMed 
Article 

Google Scholar 
Setia, R. et al. Phytoavailability and human risk assessment of heavy metals in soils and food crops around Sutlej river, India. Chemosphere 263, 128321 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Burga, D. & Saunders, K. Understanding and Mitigating Lead Exposure in Kabwe: A One Health Approach (S. Afr. Inst. Policy Res, 2019).
Google Scholar 
Ikenaka, Y., Nakayama, S. M. M., Muzandu, K. & Choongo, K. Heavy metal contamination of soil and sediment in Zambia. Afr. J. Environ. Sci. Technol. https://doi.org/10.4314/ajest.v4i11.71339 (2010).Article 

Google Scholar 
Taylor, A. A. et al. Critical review of exposure and effects: Implications for setting regulatory health criteria for ingested copper. Environ. Manag. 65, 131–159 (2020).Article 

Google Scholar 
Gummow, B., Botha, C. J., Basson, A. T. & Bastianello, S. S. Copper toxicity in ruminants: Air pollution as a possible cause. Onderstepoort J. Vet. Res. 58, 33–39 (1991).CAS 
PubMed 

Google Scholar 
Cheng, S. Effects of heavy metals on plants and resistance mechanisms. Environ. Sci. Pollut. Res. 10, 256–264 (2003).CAS 
Article 

Google Scholar 
Olobatoke, R. & Mathuthu, M. Heavy metal concentration in soil in the tailing dam vicinity of an old gold mine in Johannesburg, South Africa. Can. J. Soil Sci. 96, 299–304 (2008).Article 
CAS 

Google Scholar 
Peša, I. Between waste and profit: Environmental values on the Central African Copperbelt. Extr. Ind. Soc. https://doi.org/10.1016/j.exis.2020.08.004 (2020).Article 

Google Scholar 
Trevor, M. et al. Statistical and spatial analysis of heavy metals in soils of residential areas surrounding the Nkana Copper Mine Site in Kitwe District, Zambia. Am. J. Environ. Sustain. Dev. 4, 26–37 (2019).
Google Scholar 
Nalishuwa, L. Investigation on Copper Levels in and Around Fish Farms in Kitwe, Copperbelt Province, Zambia (Sokoine University of Agriculture, 2015).
Google Scholar 
Ikenaka, Y. et al. Heavy metal contamination of soil and sediment in Zambia. Afr. J. Environ. Sci. Technol. 4, 109–128 (2014).
Google Scholar 
Sracek, O., Mihaljevič, M., Kříbek, B., Majer, V. & Veselovský, F. Geochemistry and mineralogy of Cu and Co in mine tailings at the Copperbelt, Zambia. J. Afr. Earth Sci. 57, 14–30 (2010).ADS 
CAS 
Article 

Google Scholar 
Manchisi, J. et al. Potential for bioleaching copper sulphide rougher concentrates of Nchanga Mine, Chingola, Zambia. J. S. Afr. Inst. Min. Metall. 112, 1051–1058 (2012).
Google Scholar 
Fernández-Caliani, J. C., Barba-Brioso, C., González, I. & Galán, E. Heavy metal pollution in soils around the abandoned mine sites of the Iberian Pyrite Belt (Southwest Spain). Water Air Soil Pollut. 200, 211–226 (2009).ADS 
Article 
CAS 

Google Scholar 
Prasad, R. & Chakraborty, D. Phosphorus Basics: Understanding Phosphorus Forms and Their Cycling in the Soil 1–4 (Alabama Coop. Ext. Syst, 2019).
Google Scholar 
Verma, F. et al. Appraisal of pollution of potentially toxic elements in different soils collected around the industrial area. Heliyon 7, e08122 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Hermans, S. M., Buckley, H. L., Case, B. S., Curran-cournane, F. & Taylor, M. Bacteria as emerging indicators of soil condition. Appl. Environ. Microbiol. 83, 1–13 (2017).Article 

Google Scholar 
Ndeddy Aka, R. J. & Babalola, O. O. Identification and characterization of Cr-, Cd-, and Ni-tolerant bacteria isolated from mine tailings. Bioremediat. J. 21, 1–19 (2017).Article 
CAS 

Google Scholar 
Hassan, A., Pariatamby, A., Ahmed, A., Auta, H. S. & Hamid, F. S. Enhanced bioremediation of heavy metal contaminated landfill soil using filamentous fungi consortia: A demonstration of bioaugmentation potential. Water Air Soil Pollut. 230, 1–20 (2019).Article 
CAS 

Google Scholar 
Zhou, L. et al. Restoration of rare earth mine areas: organic amendments and phytoremediation. Environ. Sci. Pollut. Res. 22, 17151–17160 (2015).CAS 
Article 

Google Scholar 
Kapungwe, E. M. Heavy metal contaminated water, soils and crops in peri urban wastewater irrigation farming in Mufulira and Kafue towns in Zambia. J. Geogr. Geol. 5, 55–72 (2013).
Google Scholar 
Sandell, E. Post-Mining Restoration in Zambia (Swedish University of Agricultural Sciences, 2020).
Google Scholar 
Kumar, V., Pandita, S. & Setia, R. A meta-analysis of potential ecological risk evaluation of heavy metals in sediments and soils. Gondwana Res. 103, 487–501 (2022).ADS 
CAS 
Article 

Google Scholar 
Kumar, V., Sihag, P., Keshavarzi, A., Pandita, S. & Rodríguez-Seijo, A. Soft computing techniques for appraisal of potentially toxic elements from Jalandhar (Punjab), India. Appl. Sci. 11, 8362 (2021).CAS 
Article 

Google Scholar 
Setia, R. et al. Assessment of metal contamination in sediments of a perennial river in India using pollution indices and multivariate statistics. Arab. J. Geosci. 14, 1–9 (2021).Article 
CAS 

Google Scholar 
Kumar, V. et al. Pollution assessment of heavy metals in soils of India and ecological risk assessment: A state-of-the-art. Chemosphere 216, 449–462 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Environmental Council of Zambia. Environment Outlook Report in Zambia (2008).Kasali, G. Clacc Capacity Strengthening in the Least Developed Countries. CLACC Working Paper (2008).Ettler, V., Mihaljevič, M., Kříbek, B., Majer, V. & Šebek, O. Tracing the spatial distribution and mobility of metal/metalloid contaminants in Oxisols in the vicinity of the Nkana copper smelter, Copperbelt province, Zambia. Geoderma 164, 73–84 (2011).ADS 
CAS 
Article 

Google Scholar 
Cook, J. M. et al. The comparability of sample digestion techniques for the determination of metals in sediments. Mar. Pollut. Bull. 34, 637–644 (1997).CAS 
Article 

Google Scholar 
Güven, D. E. & Akinci, G. Comparison of acid digestion techniques to determine heavy metals in sediment and soil samples. Gazi Univ. J. Sci. 24, 29–34 (2011).
Google Scholar 
Jha, P. et al. Predicting total organic carbon content of soils from Walkley and Black analysis. Commun. Soil Sci. Plant Anal. 45, 713–725 (2014).CAS 
Article 

Google Scholar 
Walkley, A. & Black, I. A. A critical examination of rapid method for determining organic carbon in soil. Soil Sci. 63, 251–254 (1974).ADS 
Article 

Google Scholar 
Ure, A. M. Methods of analysis for heavy metals in soils. In Heavy Metals Soils (ed. Alloway, B. J.) 58–102 (Springer, 1995).Chapter 

Google Scholar 
Staniland, S. et al. Cobalt uptake and resistance to trace metals in comamonas testosteroni isolated from a heavy-metal contaminated site in the Zambian Copperbelt. Geomicrobiol. J. 27, 656–668 (2010).CAS 
Article 

Google Scholar 
Ajmone-Marsan, F. & Biasioli, M. Trace elements in soils of urban areas. Water Air Soil Pollut. 213, 121–143 (2010).ADS 
CAS 
Article 

Google Scholar 
Adriano, D. C. Trace elements in terrestrial environments. J. Environ. Qual. 32, 374 (2003).
Google Scholar 
Adriano, D. C. Trace Elements in Terrestrial Environments: Biogeochemistry, Bioavailability and Risks of Metals (Springer, 2001).Book 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.r-project.org/ (2020).Wickham, H. Ggplot2: Elegant Graphics for Data Analysis. Springer, New York, NY, USA, (2009).Hakanson, L. Ecological risk index for aquatic pollution control. A sedimentological approach. Water Res. 14, 975–1001 (1980).Article 

Google Scholar 
Muller, G. Index of geoaccumulation in sediments of the Rhine River. Geojournal 2, 108–118. (1969).
Google Scholar 
Usero, J., A. Garcia and J. Fraidias, 2000. Andalicia Board, Environmental Counseling. 1st Edn., Seville, Editorial, pp: 164.Sikamo, J., Mwanza, A. & Mweemba, C. Copper mining in Zambia—history and future. J. S. Afr. Inst. Min. Metall. 116, 6–8 (2016).Article 
CAS 

Google Scholar 
DR Congo: copper production 2010–2020|Statista. https://www.statista.com/statistics/1276790/copper-production-in-democratic-republic-of-the-congo/.Lydall, M. I. & Auchterlonie, A. The Southern African Institute of Mining and Metallurgy 6th Southern Africa base metals conference 2011. The Democratic Republic of Congo and Zambia: A growing global ‘Hotspot’ for copper-cobalt mineral investment and explo. In The Southern African Institute of Mining and Metallurgy 25–38 (2011).Worlanyo, A. S. & Jiangfeng, L. Evaluating the environmental and economic impact of mining for post-mined land restoration and land-use: A review. J. Environ. Manag. 279, 111623 (2021).CAS 
Article 

Google Scholar 
Shengo, M. L., Kime, M. B., Mambwe, M. P. & Nyembo, T. K. A review of the beneficiation of copper-cobalt-bearing minerals in the Democratic Republic of Congo. J. Sustain. Min. 18, 226–246 (2019).Article 

Google Scholar 
Tembo, B. D., Sichilongo, K. & Cernak, J. Distribution of copper, lead, cadmium and zinc concentrations in soils around Kabwe town in Zambia. Chemosphere 63, 497–501 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Tveitnes, S. Soil productivity research programme in the high rainfall areas in Zambia. Agricultural University of Norway (1981).Esshaimi, M., El Gharmali, A., Berkhis, F., Valiente, M. & Mandi, L. Speciation of heavy metals in the soil and the mining residues, in the Zinclead Sidi Bou Othmane Abandoned mine in Marrakech area. Linnaeus Eco-Tech https://doi.org/10.15626/eco-tech.2010.102 (2017).Article 

Google Scholar 
Vítková, M. et al. Primary and secondary phases in copper-cobalt smelting slags from the Copperbelt Province, Zambia. Mineral. Mag. 74, 581–600 (2010).Article 
CAS 

Google Scholar 
Van Brusselen, D. et al. Metal mining and birth defects: A case-control study in Lubumbashi, Democratic Republic of the Congo. Lancet Planet. Health 4, e158–e167 (2020).PubMed 
Article 

Google Scholar 
Peša, I. Between waste and profit: Environmental values on the Central African Copperbelt. Extr. Ind. Soc. 8, 100793 (2021).
Google Scholar 
Muleya, F. et al. Investigating the suitability and cost-benefit of copper tailings as partial replacement of sand in concrete in Zambia: An exploratory study. J. Eng. Des. Technol. 19, 828–849 (2020).
Google Scholar 
Namweemba, M. G. Mining Induced Heavy Metal Soil and Crop Contamination in Chililabombwe on the Copperbelt of Zambia (University of Zambia, 2017).
Google Scholar 
Colombo, C., Palumbo, G., He, J.-Z., Pinton, R. & Cesco, S. Review on iron availability in soil: Interaction of Fe minerals, plants, and microbes. J. Soils Sediments 14, 538–548 (2014).CAS 
Article 

Google Scholar 
Barsova, N., Yakimenko, O., Tolpeshta, I. & Motuzova, G. Current state and dynamics of heavy metal soil pollution in Russian Federation—A review. Environ. Pollut. 249, 200–207 (2019).CAS 
PubMed 
Article 

Google Scholar 
WHO/FAO. Food additives and contaminants. Joint FAO. WHO Food Stand. Program. ALINORM 1, 1–289 (2001).
Google Scholar 
Sracek, O. et al. Mining-related contamination of surface water and sediments of the Kafue River drainage system in the Copperbelt district, Zambia: An example of a high neutralization capacity system. J. Geochem. Explor. 112, 174–188 (2012).CAS 
Article 

Google Scholar 
Hasimuna, O. J., Chibesa, M., Ellender, B. R. & Maulu, S. Variability of selected heavy metals in surface sediments and ecological risks in the Solwezi and Kifubwa Rivers, Northwestern province, Zambia. Sci. Afr. 12, e00822 (2021).
Google Scholar 
Kříbek, B. Mining and the environment in Africa. Conserv. Lett. 7, 302–311 (2011).
Google Scholar 
Crommentuijn, T., M.D.Polder & Plassche, E. J. van de. Maximum Permissible Concentrations and Negligible Concentrations for metals, taking background concentrations into account. National Institute of Public Health and the Environment Bilthoven, The Netherlands (1997).Maboeta, M. S., Oladipo, O. G. & Botha, S. M. Ecotoxicity of mine tailings: Unrehabilitated versus rehabilitated. Bull. Environ. Contam. Toxicol. 100, 702–707 (2018).CAS 
PubMed 
Article 

Google Scholar 
Festin, E. S., Tigabu, M., Chileshe, M. N., Syampungani, S. & Odén, P. C. Progresses in restoration of post-mining landscape in Africa. J. For. Res. 30, 381–396 (2019).Article 

Google Scholar 
Volk, J. & Yerokun, O. Effect of application of increasing concentrations of contaminated water on the different fractions of Cu and Co in sandy loam and clay loam soils. Agriculture 6, 64 (2016).Article 
CAS 

Google Scholar 
Pietrini, F. et al. Effect of different copper levels on growth and morpho-physiological parameters in giant reed (Arundo donax L.) in semi-hydroponic mesocosm experiment. Water (Switzerland) 11, 1837 (2019).CAS 

Google Scholar 
EPA. Ecological Soil Screening Level for Iron Interim Final 211 (US Environ. Prot. Agency – Off. Solid Waste Emerg., 2005).
Google Scholar  More

An experiment in secondary forests in the Democratic Republic of the Congo finds that calcium, an overlooked soil nutrient, is scarcer than phosphorus, and represents a potentially greater limitation on tropical forest growth.Ecology can reveal distributional patterns and dynamics in nature. One approach used is studying the elemental composition of plants, which has been linked to ecological processes such as growth, diversity or water use efficiency. More recently, elemental composition has been detected as a cofactor in governing the carbon sink capacity of plants, and therefore climate change mitigation1,2,3. This discovery has added an extra layer of urgency to the field, which now aims to better understand and predict global change. The study of nitrogen and/or phosphorus has until now received most of the attention of plant ecologists: nitrogen is the most abundant element in dry leaves after hydrogen and carbon, forms the main structure of proteins and is strongly linked to photosynthesis4. Phosphorus represents around one-tenth of nitrogen abundance in leaves and is key in energy storage and nucleic acids. However, although these represent only two of the many chemical elements that are in flux throughout ecosystems, whether others may have a dominant role in ecosystem dynamics is an open question. Writing in Nature Ecology & Evolution, Bauters et al.5 share some evidence to motivate broadening out from the dominant focus on nitrogen and phosphorus in terrestrial ecology, by revealing a crucial limiting role of calcium in the dynamics of tropical forest succession. More