Philippa Kaur

Birds and housing62 males and 8 females of house sparrows were caught with mist nets in September and October 2019 in several sites in Kraków, Poland. Before releasing them to the outdoor aviary located on the campus of the Jagiellonian University, Kraków, Poland, each bird was weighed and banded with a metal band. The aviary measured 3.5 m in width, 10.0 m in length, 2.5 m in height, and was outfitted with trees, bushes, perches, wooden shelters, a water source, and food dishes. Initially, birds were maintained with water and a mixture of seeds: wheat, barley, millet, and sunflower seeds, provided ad libitum. Additionally, they had access to sand with shells and sepia.Experimental designAfter a few weeks of acclimation to captivity, the aviary was divided into two separate parts (3.5 × 5 m): aviary no. 1 (A1) and aviary no. 2 (A2). At the same time male individuals were assigned to two crossed experimental treatments, ensuring that in each aviary birds originated from all sampled populations. The experiment comprised of two different treatments conducted simultaneously—one designed to simulate a deficiency in an environmental factor influencing colouration (the quality of available food), the other—to introduce physiological stress and facilitate trade-offs in the allocation of resources limited by the first treatment (an immune response induced by a bacteria-derived compound, S1).The dietary manipulation was achieved by feeding one group of birds with a low-quality protein food (diet reduced in exogenous amino acids, namely phenylalanine and tyrosine content, which are precursors essential for melanin synthesis; PT-reduced diet), and the other one with a wholesome diet (control diet). At the same time, two levels of immune challenge were achieved within each dietary group, by injecting half of the birds with either lipopolysaccharide (LPS) from the cell wall of Escherichia coli, or a 0.9% saline vehicle (as a control). Four females were placed in each group of males to alleviate interspecific conflicts occurring in all-male sparrow flocks, but they did not take part in the experiments. After three weeks of experiment, birds housed in A1 were moved to A2, whereas birds from A2 were moved to A1.Immune challengeBefore receiving injections, birds were first weighed and then transferred from the outdoor aviary to the laboratory. 31 house sparrows (from both dietary groups) were injected intraperitoneally with 0.026 mg LPS (serotype O55:B5, Sigma-Aldrich) diluted in 0.1 mL of 0.9% saline vehicle, so that each bird received a dose of ca. 1 mg/kg body mass, which had previously been shown to induce sickness behaviour in another passerine, the white-crowned sparrow, Zonotrichia leucophrys55. 31 control males were injected with the same volume (0.1 mL) of 0.9% saline vehicle. All individuals were injected twice throughout the experiment with an interval of three weeks between the injections. Birds were always injected at the same time in the morning and early afternoon (between 9:00 am and 12:30 pm).Diet manipulationDuring the six weeks of the experiment (S1), birds received synthetic diet ad libitum, which constituted of a mixture of protein (WPC80, free amino acids and whey protein isolate BiPRO GMP 9000 (Agropur Inc., Appleton, USA)), fats, carbohydrates, and fiber30. The ingredients were thoroughly mixed to produce small pellets (6 mm in diameter) that the sparrows consumed readily. The experimental diet had phenylalanine and tyrosine at 42% (N = 32) of their level in the control diet (N = 30)30. The food pellets were prepared by ZooLab (zoolab.pl/en/home, Sędziszów, Poland). Each bird was weighed before and after the experiment to monitor potential effects of diet on body mass of each animal. Following the experiment, during the next three consecutive days, the amount of food consumed by passerines within every 24 h (starting from 10 am each day to 10 am next day) was noted for both compartments of the aviary. Because of different numbers of individuals per aviary, an overall weight of food consumed in A1 and A2 was calculated per individual, respectively.Feathers samplingMoult of the black bib feathers was stimulated at the end of the moulting period occurring in natural conditions in early November. At day 1 of the dietary/immunological experiment (S1) a small area of the bib (around 25 mm2) was plucked from each male sparrow held in A1. At day 2 the same procedure was performed on individuals from A2. The time difference is orders of magnitude smaller than the timescale of feather growth and hence it would not affect the results in any way.Because the feather growth rate may differ during melanogenesis, with consequences for final colouration (if feathers grow at a faster rate, pigments may be deposited over a larger surface and therefore result in less intense colouration56, we measured the rate of feather development during the course of the experiment. After three weeks of the experiment, three feathers from the upper, central, and lower region of the previously plucked bib were plucked once again. The mass of the collected feathers was determined to the nearest 0.01 mg (XP26 Micro Balance, Mettler-Toledo, Greinfensee, Switzerland). The experiment was completed after six weeks after fully regrown and developed feathers from the bib and PC2 were sampled the second time (S1). Three feathers from the central part of previously plucked bib region were collected to perform transmission electron microscopy (TEM) imaging, whereas the feathers obtained from the rest of the regrown bib area were subjected to electron paramagnetic resonance (EPR) spectroscopy and feather microstructure analyses (greater spatial density of melanized barbs or barbules may affect colouration17.Feathers measurementsReflectance measurementsAn USB4000 spectrophotometer (range 300–700 nm) with the PX-2 Pulsed Xenon Lamp (Ocean Optics, Dunedin, FL, USA) and a bifurcated probe with 7 × 400 μm optical fibres, equipped with a permanently attached 3 mm long black collar, was used to quantify the brightness of the bib feathers collected at the end of the experiment. The measurements were taken with 90 ms integration time and the probe held at 90° to a feather’s surface. Calibration measurements of a Spectralon white standard (Ocean Optics. Largo, FL, USA) were taken every 15 min during measurements. The order in which the samples were measured was randomized in terms of belonging to the experimental group. From each sample (N = 62), seven feathers were chosen and stacked in one pile on a piece of black paper. Ten reflectance measurements were taken on each pile, avoiding distal, brighter parts of the feathers. The obtained spectra were averaged and smoothed in the package ‘pavo’57. Brightness was calculated as a sum of the reflectance values over all wavelengths of a spectrum, and its lower values were interpreted as those indicative of a more melanin-rich feathers (i.e., absorbing more light).Feather developmentEach feather (3 per individual; N = 62 individuals) was laid on a white card and covered by a microscope slide to flatten the naturally curved feathers. Digital photographs were taken using camera (Canon EOS 7D) and imported to ImageJ v1.52a Software (National Institutes of Health, USA). The lengths of fully developed and undeveloped (still in sheath) parts of each feather were measured. To estimate the degree of a feather’s development, the length of the developed part of the vane was divided by its total length (quill with rachis plus the developed vane, Fig. 4A).Figure 4House sparrow feathers sampled from bib after three weeks of the experiment. Feathers during development (A), a TEM cross-sections of feather sampled from bib after the experiment (B).Full size imageFeather densityBarb density measurements were performed on the sampled regrown black bib feathers (N = 2–3 for each individual; N = 62 individuals), but because of their sparser structure we calculated the number of non-down (i.e., rigid) barbs on both sides of the vane, and divided this number by two (to obtain an average single-sided number of barbs) and then by the length of the rachis.Melanosome density (TEM)Feathers sampled from the bib of male sparrows (N = 62) were fixed for transmission electron microscopy (TEM) analysis in a mixture of 0.25 M sodium hydroxide and 0.1% Tween for 20 to 30 min on a bench-top shaker. Next, the feathers were treated with formic acid and ethanol in the ratio of 2:3 for 2.5 h and dehydrated twice for 20 min in 100% ethanol. Samples were embedded in a mixture of the PolyBed 812 resin (20 ml), DDSA (9 ml), NMA (12 ml) and DMP-30 (0.82 ml). Resin infiltration was gradual from 15% resin content in ethanol through 50%, 70% to 100% without alcohol. Each step lasted for 24 h. Then, the feathers were placed in silicone embedding moulds (Agar Scientific) and transferred to an oven. The polymerization proceeded at the temperature of 60 °C for 16 h. The epoxy resin blocks were then trimmed to get rid of excess resin. The surface of each block was prepared by its trimming, starting from the end of the feather, to approximately 5 mm using a glass knife. Next, ultrathin sections (70 nm) were cut with a diamond knife (DIATOME A. G., Berno, Switzerland) on a microtome (UC7, Leica, Wetzlar, Germany) and collected on single slot grids coated with a formvar film. The sections were then contrasted in uranyl acetate and lead citrate for 3 min. They were viewed and photographed with a transmission electron microscope (TEM) JEOL 2100HT (Jeol Ltd, Tokyo, Japan) for the purpose of investigating the number and density of the embedded pigment granules. For each individual three photographs of the cross-sections from a similar feather region were selected. Melanosome density was measured as the number of melanin granules observed in the barb cross-section divided by its area. Images were analysed using Adobe Photoshop (cross-sections area) and ImageJ (number of melanosomes, Fig. 4B).Melanin content: electron paramagnetic resonance (EPR) spectroscopyQuality and quantity of melanin pigments58 in individual feather samples obtained from the bib of house sparrows (N = 57) were characterized using a Varian E3 spectrometer (Varian, Sunnyvale, LA, USA) equipped with a rectangular resonance (TE 102) cavity. Five milligrams of feathers per individual were placed inside the Wilmad finger quartz dewar WG-816-Q (Rototec-Spintec GmbH, Griesheim, Germany). Prior to inserting the vessel into the resonance cavity of the EPR spectrometer, feathers were pressed down the quartz finger to a height of approximately 0.5 cm to ensure comparable volumes of each sample. Measurements were performed at room temperature, at X-band (9.26–9.27 GHz frequency), using the following parameters: magnetic field range 3240–3340 Gs, microwave power 1 mW, modulation frequency 100 kHz, modulation amplitude and time constant—5 Gs and 0.3 s for quantitative analysis, 1 Gs and 0.1 s for qualitative analysis. An EPR signal was recorded as its first derivative, averaged from three consecutive scans, lasting 160 s each (giving a total of 480 s of scan time per EPR spectrum). Then, the following parameters were measured: peak-to-peak amplitude, area under the microwave absorption curve (the integral intensity of the recorded signal) and linewidth of the EPR absorption curve (ΔH;59).Statistical analysesStatistical analysis was performed in R (version 4.0.2,60) using a two-way ANOVA test, with bird’s diet (control vs. PT-reduced) and applied immune challenges (LPS vs. saline-injections) as the independent variables. The following parameters were used as the dependent variables: feathers reflectance (brightness), feather growth rate, feather density (number of barbs per mm), and melanisation level (expressed as the EPR spectrum amplitude measured in arbitrary units [a.u.]). The density of melanosomes was analysed by fitting a linear mixed-effects model. In this model, melanosome density was used as the dependent variable, with diet, immunological challenge, and slice ID as independent variables, and individual ID as a random-effect term. Additionally, to assess the reliability of measurements, the intraclass correlation coefficient (i.e., technical repeatability) was calculated. The models’ residuals were checked for normality and homoscedasticity. Mean food consumption per individual was analysed by the Friedman test. Body mass before and after the experiment was analysed by fitting a linear mixed-effect model. Body mass was used as the dependent variable, whereas diet, immunological challenge, and time as the independent variables, and individual ID as a random-effect term. The model included the following interaction terms: time × diet, time × injection, and diet × injection, and was reduced by removing the non-significant interactions. Results are reported with appropriate statistical tests and estimates (accompanied by standard errors) signifying relevant factor contrasts (relative to the reference group, which in all analyses was diet: control; injection: LPS, body mass: before experiment).
Ethical noteAll applicable national and institutional guidelines for the care and use of animals were followed. The research was performed under permit no. 25/2019 (with a supplementary permit no. 78/2020) from the 2nd Local Institutional Animal Care and Use Committee in Kraków. More

UTrack atmospheric moisture tracking modelThe UTrack atmospheric moisture tracking model is a novel Lagrangian model that tracks parcels of moisture forward in three-dimensional space9. UTrack is the first moisture tracking model to employ ERA5 reanalysis data8. The basic principle of the model is that for each mm of evaporation, a certain number of “moisture parcels” is released and subsequently tracked through time and space. At each time step, the moisture budget of the parcels is updated based on evaporation and precipitation at the respective time and location, meaning that for each location of evaporation, a detailed image of the “footprint” of evaporation can be created. All types of evapotranspiration are included, and here is simply called evaporation.For each mm of evaporation, 100 parcels are released 50 hPa above the surface height at random spatial locations within each 0.25° grid cell of input evaporation data. The trajectories of the parcels are based on interpolated three-dimensional ERA5 wind speed and wind direction data, which also have a horizontal resolution of 0.25° and consist of 25 pressure layers in the atmospheric column. The spatial coordinates of each parcel are updated at each time step of 0.1 h. Also, at each time step, there is a certain probability that a parcel is redistributed randomly along the atmospheric column such that, on average, every parcel is redistributed every 24 h (see methods section Moisture recycling dataset: validation and uncertainties below for further details). The relative probability of the new position in the atmospheric column is scaled with the vertical moisture profile. Parcels are tracked for 30 days or until 99% of their moisture has precipitated.To allocate a certain fraction of any moisture parcel to precipitation events at the current time and location, ERA5 hourly total precipitation (P) and total precipitable water (TPW) are interpolated to the simulation time step of 0.1 h. The amount of moisture that precipitates at a certain time step equals the amount of precipitation at that time step over the total precipitable water in the atmospheric water column (P/TPW). Specifically, precipitation A in mm per time step at location x, y at time t that originated as evaporation from a particular source is described as:$${A}_{x,y,t}={P}_{x,y,t}frac{{W}_{{{{{{{{rm{parcel,t}}}}}}}}}{E}_{{{{{{{{rm{source,t}}}}}}}}}}{{{{{{mathrm{TP}}}}}}{{{{{{mathrm{W}}}}}}}_{x,y,t}}$$
(1)
with P being precipitation in mm at time step t, Wparcel,t (mm) the amount of moisture in the parcel of interest, Esource,t the fraction of moisture present in the parcel at time t that has evaporated from the source, and TPW (mm) the precipitable water in the atmospheric water column. The moisture content of parcels is updated each time step using evaporation and precipitation at its current location:$${W}_{{{{{{{{rm{parcel,t}}}}}}}}}={W}_{{{{{{{{rm{parcel,t-1}}}}}}}}}+({E}_{{{{{{{{rm{x,y,t}}}}}}}}}-{P}_{{{{{{{{rm{x,y,t}}}}}}}}})frac{{W}_{{{{{{{{rm{parcel,t-1}}}}}}}}}}{{{{{{mathrm{TP}}}}}}{{{{{{mathrm{W}}}}}}}_{{{{{{{{rm{x,y,t}}}}}}}}}}$$
(2)
The moisture (fraction) that has evaporated from the source is updated as follows:$${E}_{{{{{{{{rm{source,t}}}}}}}}}=frac{{E}_{{{{{{{{rm{source,t-1}}}}}}}}}{W}_{{{{{{{{rm{parcel,t-1}}}}}}}}}{A}_{x,y,t}}{{W}_{{{{{{{{rm{parcel,t}}}}}}}}}}$$
(3)
The moisture flow mij from evaporation in cell i to precipitation in cell j is aggregated on a monthly basis (mm/month), where [x, y] ∈ j becomes:$${m}_{ij}=mathop{sum }limits_{t=0}^{{{{{{{{rm{month}}}}}}}}}{A}_{j,t}frac{{E}_{i,t}}{{W}_{i,t}}$$
(4)
with Wi,t being the tracked amount of moisture from the source cell i at time t. These simulations were performed for all evaporation on Earth during 2008–2017. The results were then aggregated on a mean-monthly basis to produce monthly means, and stored at 0.5 degree resolution. This dataset can be downloaded from ref. 53. For details on how to process the data, we refer to the accompanying paper by ref. 3.Moisture recycling dataset: validation and uncertaintiesAs with all moisture recycling simulations, the ones used in this study rely on a number of assumptions that may affect the moisture recycling rates. All offline moisture recycling models use atmospheric model output to simulate the path of evaporation through the atmosphere to the location where it precipitates. Therefore, there are two sources of uncertainty that affect the moisture recycling estimates: (1) the quality of the atmospheric forcing data and (2) the assumptions in the moisture tracking model.Tuinenburg and Staal (2020)9 explored these sources of uncertainty for a number of locations globally. The effects of a decrease in the quality of the atmospheric forcing data were most important in the vertical resolution of the atmospheric data: the forcing data should have enough vertical levels to resolve any vertical shear in atmospheric moisture transport. If the forcing data has a low vertical resolution, the moisture tracking model is forced with the mean atmospheric flow over a number of layers. In many regions, there are surface moisture flows that are in a different direction than the moisture flow aloft, resulting in a very small vertically integrated transport, which would distort the simulation of atmospheric moisture transport. Compared to the vertical resolution of the forcing data, the horizontal and temporal resolutions were less important in order to keep errors as small as possible. Because of the importance of this high vertical resolution, it was recommended9 to use the ERA5 dataset8 as its forcing dataset, as this currently is the atmospheric reanalysis dataset with the highest vertical resolution.In addition, the change of ERA-interim to ERA5 resulted in a much better land-surface scheme with monthly varying vegetation and better bare soil evaporation. Also, many more observations are assimilated, which results in a better precipitation product compared to ERA-interim. Following this, the tracking of atmospheric moisture using ERA5 allows for a better quality of the atmospheric moisture cycle than before. But, of course, also the already high horizontal resolution of 0.5∘ × 0. 5∘ has the limitation that very localized moisture recycling features like orography and locally varying land use cannot be resolved. Out of these reasons, the uncertainty in the evaporation estimates is a lot larger than that in the precipitation estimates, because of the lack of global evaporation measurements and the difficulty in measuring evaporation in general54,55.There are also uncertainties due to the assumptions in the moisture tracking model that can be split into a category of simulation assumptions and physical assumptions. The simulation assumptions include model formulation (Eulerian vs. Lagrangian model set-ups), time step lengths, number of parcels released, and types of interpolation. Of these simulation assumptions, the most important aspects were the model formulation, with Lagrangian models better able to resolve complex terrain and atmospheric flows. For the other model assumptions (see methods section UTrack atmospheric moisture tracking model), it was chosen to simulate with the highest level of precision before any more information (e.g., more parcels) would no longer affect evaporation footprints and moisture recycling statistics (see ref. 9 for further details). Even though the ERA5 dataset is known to have some precipitation biases in the tropics, the results of UTrack (forced by ERA5) have recently been validated across the tropics by independent measurements of deuterium excess, a measure of a stable isotope that depends on terrestrial precipitation recycling56. UTrack estimates and isotope-based estimates of terrestrial moisture recycling corresponded, especially in tropical rainforests (Kendall’s (overline{tau }=0.52)56), which are found to be moisture recycling hubs on a global scale.Network constructionMotivated by the network-like structure of the data, we here employ a network perspective to study moisture flows. Hence, nodes in such a network are grid cells on a regular spherical grid and edges represent the moisture transported. However, interpreting the dataset directly as a weighted network is neither computationally feasible nor does a weighted network allow for identifying motifs, the building blocks of complex networks17. We, therefore, aim for an approach utilizing an unweighted network.As shown in Fig. S1, moisture recycling strengths are heterogeneously distributed over multiple powers of magnitude. Thus, it is not appropriate to just withdraw the moisture transport volume and include all moisture transport connections within the dataset as equal and unweighted links. Instead, we attempt to highlight the strongest moisture pathways and, thus, the backbone of the Earth’s moisture recycling network. To, on the one hand, include as much moisture volume as possible but also keep the absolute volume of moisture transport represented per edge as similar as possible, we decided to include edges in a data-adaptive way: we step-wise include links starting from the strongest and stop this procedure as the total moisture transport volume exceeds the variable threshold ρ. The resulting edges then represent the backbone of the global moisture recycling network. In the main text, we have shown the results for a network where all edges together represent ρ = 25% of the total moisture transport. Here and in the SI figures, we add a sensitivity analysis for ρ = 20% and ρ = 30% and find that the results are stable for this broader range of total moisture volume thresholds.Network measures and motifsThe topology of an unweighted network is typically encoded in an adjacency matrix A with elements aij indicating if there exists an edge from node i to node j (aij = 1) or not (aij = 0). The degree k of a node i describes the number of adjacent edges pointing towards or away from node i. Hence, the in-degree is defined by25$${k}_{{{{{{mathrm{in}}}}}}}^{i}=mathop{sum }limits_{i=1}^{N}{a}_{ji}$$
(5)
and out-degree is defined by25$${k}_{{{{{{mathrm{out}}}}}}}^{i}=mathop{sum }limits_{i=1}^{N}{a}_{ij}.$$
(6)
To further analyze the topology of a network and, in particular, the local connectivity patterns, we study the presence of three motifs—the feed-forward loop, the neighboring loop, and the zero loop.The feed-forward loop (FFL) consists of three nodes, A, B, and C, where nodes A and C are directly connected via a detour over node B (intermediary node). Therefore, we have two different pathways that focus on node C. Hence, this motif can be referred to as a directed lens, due to its focused flow from two nodes on one singular and its purely directed linkage. This network motif has been studied in the context of tipping elements and has been proven to facilitate tipping cascades by lowering critical thresholds19. The zero loop (ZL) is made up of a bidirectional connection of two nodes. In contrast to the FFL, where node A does not receive feedback from node C, here, both nodes are dependent on each other without a preferred direction of network flow. This facilitates tipping to a much lesser degree than the FFL motif19. The neighboring loop (NBr) is an extension of the ZL. In this case, there is an additional node connected to one of the nodes of a zero loop. Hence, there is a two-step directionality in the motif, but in contrast to the FFL, this motif is characterized by reciprocity.We count the number of motifs a certain node is involved in the network. The number of FFLs is counted as the number that a certain node is a so-called “target” node. The target node is the node, on which the triangular structure of the motif is converging to, i.e., the node that has been referred to as node C above. The ZL is a symmetric motif for the two involved nodes. Therefore, the number of ZLs of a certain node in the network is counted directly as the number of bidirectional interactions of the inspected node. Lastly, the number of NBrs of a certain node is the number of being in the center of a neighboring loop. With this procedure, each node is characterized by its number of FFLs, ZLs, and NBrs (cf. ref. 19).Motif strength and their spatially aggregated differenceTo assess the presence of motifs and, in particular, their relative frequency, we first determine the numbers of FFLs, ZLs, and NBrs per node. Subsequently, we normalize these counts by the respective maximum to obtain the motif strength, which is shown for each network motif in Fig. S5. In Fig. S5a–c, we display the motifs for the global network, and in Fig. S5d–f for the land-to-land network.To specifically characterize the focus regions by means of the network topology, we evaluate which motifs dominate in which region. Consequently, we compute the difference of the motif strengths shown in Fig. S5 and reveal the patterns shown in Fig. 2. For spatially aggregated motif strength differences (Fig. 2c, d), we then compute the average of the respective values inside the highlighted boxes.Sensitivity to link threshold ρ
The network analysis featured in the main text uses those moisture recycling edges that together represent ρ = 25% of all atmospheric moisture recycling on Earth. As we aimed to focus on the strongest moisture flows, we chose a threshold of ρ = 25% aggregating the strongest moisture transport pathways. This allows us to reveal the regions of strongest moisture connections, which are located in and close to the tropics, as we expected. Overall, the aim of this thresholding procedure is to utilize a network approach with unweighted edges but also take into account the large spread of moisture recycling strengths. To test the robustness of the results to the threshold value, we here show the same figures as above in the main text but with different thresholds ρ. Note that the error bars in Fig. 2 are based on the analysis featured in this part (the resulting differences using thresholds of ρ = 20% and ρ = 30%).Figures S6 and S7 show the in- and out-degree of the all-to-all and land-to-land network using a threshold of ρ = 20% (Fig. S6) and ρ = 30% (Fig. S7). Note that the color bar has been adjusted as the number of links differs substantially between the networks. The main difference between Figs. S6 and S7 is the greater emphasis on moisture recycling in the mid-latitudes in Fig. S7. This is a direct consequence of considering more, and thus also some weaker, links. Acknowledging this difference, we stress that especially the land-to-land patterns (Figs. S6c, d, S7c, d) are consistent. In particular, the four focus regions, as defined in the main text, stand out as the main global land-to-land moisture recycling hubs. To support this visual analysis of the in- and out-degree pattern, we furthermore compute the motif strengths for both network configurations for quantitative validation of the results.In line with the main text, we compare the FFL and ZL strength (see Fig. 2a–d). Not only the spatial patterns in our sensitivity analysis agree remarkably well with the results in the main text above, but also the focus regions remain basically the same (cf. Fig. S8 for ρ = 20% and Fig. S9 for ρ = 30% with Fig. 2). The only slight change is the shift towards a directed lens (spatially aggregated FFL and ZL strength difference) for the Amazon basin in the all-to-all network for increasing ρ (cf. Fig. S8c vs Fig. S9c vs Fig. 2c). We attribute the overproportional increase of the number of FFLs to those that include at least one oceanic grid cell to this noticeable shift. This underscores our characterization of the Amazon basin as a directed lens.The spatially aggregated FFL and NBr difference (Figs. S10, S11) is structurally the same as above, where we computed the FFL and ZL difference (see Figs. S8, S9). The spatial patterns and the aggregated values are robust against shifts of ρ. However, for the Amazon basin (AB), the number of FFLs increases overproportionally in the all-to-all network when we include more links in our analysis. In other words, the spatially aggregated FFL-strength for AB increases for higher thresholds ρ (cf. Figs. S10c, S11c and Fig. 2g).Sensitivity to the size of the focus regionsAnother aspect affecting the results is the spatial extent chosen as a focus region (i.e., the rectangles in Fig. 2). Varying the size of these rectangles affects the spatially aggregated measures. For all focus regions besides the Amazon Basin (AB), the values are not significantly affected by changing the rectangle size, as the values close to the focus regions are either coherently negative, as for the Congo Rainforest (CR) and the Indonesian Archipelago (IA), or close to zero (South Asia: SA). The AB is characterized by positive values (tendency to lensing), whereas the more southern parts along the Andes are marked by more negative (corridor/washing machine) values.Hence, we assess the stability of the results by using the spatial region covered by the Amazon rainforest (the extent of the Amazon rainforest is based on ref. 6) and compare them to the ones obtained by using the rectangle. The results featured in Fig. S12 indicate that only considering the rainforest-covered parts of the AB leads to similar or even more positive (lensing) values, confirming our conclusions that the Amazon rainforest region functions differently from the other focus regions.Notes on mapsThis paper makes use of perceptually uniform color maps developed by ref. 57. The underlying world maps have been created by cartopy58. More

Cavalleri, A., Masumoto, M., Minaei, K., Mound, L. & Ulitzka, M. R. ThripsWiki – providing information on the World’s thrips. https://thrips.info/wiki/Main_Page (2022).Kirk, W. D. J., de Kogel, W. J., Koschier, E. H. & Teulon, D. A. J. Semiochemicals for thrips and their use in pest management. Annu. Rev. Entomol. 66, 101–119. https://doi.org/10.1146/annurev-ento-022020-081531 (2021).Article 
CAS 
PubMed 

Google Scholar 
Mota-Sanchez, D. & Wise, C. J. The Arthropod Pesticide Resistance Database. Available at http://www.pesticideresistance.org (2022).von Kèler, S. Entomologisches Wörterbuch. mit besonderer Berücksichtigung der morphologischen Terminologie. 3rd edn. (Akademie, Berlin, 1963).
Google Scholar 
Pascini, T. V. & Martins, G. F. The insect spermatheca: an overview. Zoology Jena (Germany) 121, 56–71. https://doi.org/10.1016/j.zool.2016.12.001 (2017).Article 

Google Scholar 
Bode, W. Der Ovipositor und die weiblichen Geschlechtswege der Thripiden (Thysanoptera, Terebrantia). Z. Morph. Tiere (Zeitschrift für Morphologie der Tiere) 81, 1–53; https://doi.org/10.1007/BF00290072 (1975).Moritz, G. Zur Morphologie und Anatomie des Fransenflüglers Aeolothrips intermedius Bagnall, 1934 (Aeolothripidae, Thysanoptera, Insecta) III. Mitteilung: Das Abdomen. Investigation on the Morphology and Anatomy in Aeolothrips intermedius Bagnall, 1934 (Aeolothripidae, Thysanoptera, Insecta) 3. The Abdomen. Zoologische Jahrbücher/ Abteilung für Anatomie und Ontogenie der Tiere 108, 293–340 (1982).Moritz, G. Die Ontogenese der Thysanoptera (lnsecta) unter besonderer Berücksichtigung des Fransenflüglers Hercinothrips femoralis (O. M. REUTER, 1891) (Thysanoptera, Thripidae, Panchaetothripinae) VI. Imago – Abdomen. The Ontogenesis of Thysanoptera (Insecta) with Special Reference to the Panchaetothripine Hercinothrips femoralis (O. M. REUTER, 1891) (Thysanoptera, Thripidae, Panchaetothripinae) VI. Imago – Abdomen. Zoologische Jahrbücher/ Abteilung für Anatomie und Ontogenie der Tiere 119, 157–217 (1989).Heming, B. S. Postembryonic development of the female reproductive system in Frankliniella fusca (Thripidae) and Haplothrips verbasci (Phlaeothripidae) (Thysanoptera). Misc. Publ. Entomol. Soc. Am. 7, 197–234 (1970).
Google Scholar 
Heming, B. S. History of the germ line in male and female thrips. In Thrips Biology and Management. International conference on Thysanoptera: Towards Understanding Thrips Management, Vermont, edited by B. L. Parker, M. Skinner & T. Lewis (Springer, Berlin, 1995), pp. 505–535.Dallai, R., Del Bene, G. & Lupetti, P. Fine structure of spermatheca and accessory gland of Frankliniella occidentalis (Pergande) (Thysanoptera: Thripidae). Int. J. Insect Morphol. Embryol. 25, 317–330; https://doi.org/10.1016/0020-7322(95)00018-6 (1996).Jordan, K. Anatomie und Biologie der Physapoda. Zeitschrift für wissenschaftliche Zoologie 47, 541–620 (1888).
Google Scholar 
Bournier, A. L`appareil gènital femelle de Caudothrips buffai Karny et sa pompe spermatige (Thysan.). Bulletin de la Sociètè Entomologique de France 67, 203–207 (1962).Dhileepan, K. & Ananthakrishnan, T. Impact of sex-limited and alary polymorphism on spermathecal diversity and reproductive behaviour in some mycophagous Tubulifera. Proc. Indian Natl. Acad. Sci. 4, 329–336 (1987).
Google Scholar 
Bhatti, J. S. The spermatheca as a useful character for species differentiation in Coleothrips Haliday (Insecta: Terebrantia: Aeolothripidae). J. Pure Appl. Zool. 1, 111–116 (1988).
Google Scholar 
Klocke, F. Beiträge zur Anatomie und Histologie der Thysanopteren. Zeitschrift für wissenschaftliche Zoologie 128, 1 (1926).
Google Scholar 
Priesner, H. Die Thysanopteren Europas (F. Wagner Verlag, Wien, 1926–1928).Bode, W. Spermienstruktur und Spermatohistogenese bei Thrips validus Uzel (Insecta, Thysanoptera). Zoologische Jahrbücher/ Abteilung für Anatomie und Ontogenie der Tiere 109, 301–318 (1983).
Google Scholar 
Moritz, G. Thripse. Fransenflügler, Thysanoptera. 1st ed. (Westarp Wissenschaften, Hohenwarsleben, 2006).Bournier, A. Contribution à l’étude de la parthénogénèse des thysanoptères et de sa cytologie. Archives de Zoologie expérimentale et générale 93, 219–318 (1956).
Google Scholar 
Heming, B. S. Postembryonic development of the male reproductive system in Frankliniella fusca (Thripidae) and Haplothrips verbasci (Phlaeothripidae) (Thysanoptera). Misc. Publ. Entomol. Soc. Am. 7, 235–272 (1970).
Google Scholar 
Krueger, S., Jilge, M., Mound, L. & Moritz, G. B. Reproductive behavior of Echinothrips americanus (Thysanoptera: Thripidae). J. Insect Sci. 17; https://doi.org/10.1093/jisesa/iex043 (2017).Krueger, S. & Moritz, G. Sperm ultrastructure in arrhenotokous and thelytokous Thysanoptera. Arthropod Struct. Dev. 64, 101084; https://doi.org/10.1016/j.asd.2021.101084 (2021).Lewis, T. Thrips. Their biology, ecology and economic importance (Academic Press, London, 1973).Jacobs, W. & Seidel, F. Systematische Zoologie, Insekten (Fischer, 1975).
Google Scholar 
Gotoh, A., Ito, F. & Billen, J. Vestigial spermatheca morphology in honeybee workers, Apis cerana and Apis mellifera, from Japan. Apidologie 44, 133–143; https://doi.org/10.1007/s13592-012-0165-6 (2013).Gotoh, A., Billen, J., Hashim, R. & Ito, F. Degeneration patterns of the worker spermatheca during morphogenesis in ants (Hymenoptera: Formicidae). Evol. Dev. 18, 96–104; https://doi.org/10.1111/ede.12182 (2016).Schoeters, E. & Billen, J. The importance of the spermathecal duct in bumblebees. J. Insect Physiol. 46, 1303–1312; https://doi.org/10.1016/S0022-1910(00)00052-4 (2000).Gobin, B., Ito, F., Peeters, C. & Billen, J. Queen-worker differences in spermatheca reservoir of phylogenetically basal ants. Z. Zellforsch (Zeitschrift für Zellforschung und Mikroskopische Anatomie) 326, 169–178; https://doi.org/10.1007/s00441-006-0232-2 (2006).Gobin, B., Ito, F., Billen, J. & Peeters, C. Degeneration of sperm reservoir and the loss of mating ability in worker ants. Die Naturwissenschaften 95, 1041–1048; https://doi.org/10.1007/s00114-008-0420-x (2008).Gotoh, A., Billen, J., Hashim, R., Ito, F. Comparison of spermatheca morphology between reproductive and non-reproductive females in social wasps. Arthropod. Struct. Dev. 37, 199–209; https://doi.org/10.1016/j.asd.2007.11.001 (2008).Gotoh, A., Billen, J., Tsuji, K., Sasaki, T. & Ito, F. Histological study of the spermatheca in three thelytokous parthenogenetic ant species, Pristomyrmex punctatus, Pyramica membranifera and Monomorium triviale (Hymenoptera: Formicidae). Acta Zool. 93, 200–207; https://doi.org/10.1111/j.1463-6395.2010.00498.x (2012).Buffa, P. Studi intorno al ciclo partenogenetico dell´ Heliothrips haemorrhoidales (Boúche). REDIA 7, 71–109 (1911).
Google Scholar 
Bene, G. D., Cavallo, V., Lupetti, P. & Dallai, R. Ultrastructure of the accessory gland in the parthenogenetic thrips Heliothrips haemorrhoidalis (Bouché) (Thysanoptera. Thripidae). Int. J. Insect Morphol. Embryol. 27, 255–261; https://doi.org/10.1016/S0020-7322(98)00018-X (1998).Nakao, S. & Yabu, S. Ethological and chemical discrimination between thelytokous and arrhenotokous Thrips nigropilosus Uzel, with discussion of taxonomy. Jpn. J. Appl. Entomol. Zool. 42, 77–83; https://doi.org/10.1303/jjaez.42.77 (1998).Article 
CAS 

Google Scholar 
Arakaki, N., Miyoshi, T. & Noda, H. Wolbachia-mediated parthenogenesis in the predatory thrips Franklinothrips vespiformis (Thysanoptera: Insecta). Proc. R. Soc. B Biol. Sci. 268, 1011–1016; https://doi.org/10.1098/rspb.2001.1628 (2001).Kumm, S. & Moritz, G. First detection of Wolbachia in arrhenotokous populations of thrips species (Thysanoptera: Thripidae and Phlaeothripidae) and its role in reproduction. Environ. Entomol. 37, 1422–1428. https://doi.org/10.1603/0046-225X-37.6.1422 (2008).Article 
PubMed 

Google Scholar 
Moritz, G. The biology of thrips is not the biology of their adults: a developmental view. In Thrips and Tospovirus: Proceedings of the 7th International Symposium on Thysanoptera, edited by L. A. Mound & R. Marullo (Australian National Insect Collection CSIRO, Canberra, 2002), pp. 259–267.Moritz, G., Schäfer, E., Kumm, S., Steller, A. & Tschuch, G. D. Alien-Thrips: Suocerathrips linguis – Biologie und Verhalten. Mitteilungen der Deutschen Gesellschaft für allgemeine und angewandte Entomologie 14, 177–181 (2004).
Google Scholar 
Gehlsen, U. Ernährungssystem, Verhalten und Wehrsekret des subsozialen Phlaeothripinen Suocerathrips linguis MOUND & MARULLO, 1994 (Insecta, Thysanoptera, Tubulifera). PhD-thesis (Martin-Luther University Halle-Wittenberg, Germany, 2009).Kumm, S. Reproduction, progenesis and embryogenesis of thrips (Thysanoptera, Insecta). Dissertation. Martin-Luther-Universität Halle-Wittenberg, 2002.Krueger, S., Mound, L. A. & Moritz, G. B. Offspring sex ratio and development are determined by copulation activity in Echinothrips americanus MORGAN 1913 (Thysanoptera: Thripidae). J. Appl. Entomol. 140, 462–473; https://doi.org/10.1111/jen.12280 (2016).Oetting, R. D., Beshear, R. J., Liu, T.-X., Braman S. K., & Baker, J. R. Biology and identification of thrips on greenhouse ornamentals. Univ. Ga. Res. Bull. 414, 20 (1993).Mound, L. A., Nielsen, M.-C. & Hastings, A. Thysanoptera Aotearoa. Thrips of New Zealand. Available at https://keys.lucidcentral.org/keys/v3/nz_thrips/index.html.Schindelin, J. et al. Fiji an open-source platform for biological-image analysis. Nat. Methods 9, 676–682. https://doi.org/10.1038/nmeth.2019 (2012).Article 
CAS 
PubMed 

Google Scholar 
Uzel, H. Monographie der Ordnung Thysanoptera (Nabu Oress, Charleston SC, United States, 1895).Marzo L. De. Dettagli anatomici dei genitali interni in Melanthrips fuscus (Sulzer) e altri tisanotteri. Entomologica 36, 109–119. https://doi.org/10.15162/0425-1016/747 (2002).Melis, A. Nuove osservazioni anatomo-istologiche sui diversi stati postembrionali del Liothrips oleae Costa. REDIA 21, 263–334 (1934).
Google Scholar 
van der Kooi, C. J. & Schwander, T. On the fate of sexual traits under asexuality. Biol. Rev. 89, 805–819. https://doi.org/10.1111/brv.12078 (2014).Article 
PubMed 

Google Scholar 
Sloan, N. S. & Simmons, L. W. The evolution of female genitalia. J. Evol. Biol. 32, 882–899. https://doi.org/10.1111/jeb.13503 (2019).Article 
PubMed 

Google Scholar 
Buffa, P. Tisanotteri esotici esistenti nel Museo Civico di Storia Naturale di Genova. REDIA 5, 157–172 (1909).
Google Scholar 
Osborn, H. Note on a New Species of Phloeothrips, with description. Proc. Iowa Acad. Sci. 3, 228 (1895).
Google Scholar 
Mound, L. A. & Marullo, R. New thrips on mother-in-law`s tongue. Entomol. Mon. Mag. 130, 95–98 (1994).
Google Scholar 
Zur Strassen, R. Die terebranten Thysanopteren Europas und des Mittelmeer-Gebietes (Goecke und Evers, Keltern, 2003).Davies, R. G. The postembryonic development of the female reproductive system in Limothrips cerealium Haliday (Thysanoptera: Thripidae). Proc. Zool. Soc. Lond. 136, 411–437. https://doi.org/10.1111/j.1469-7998.1961.tb05883.x (1961).Article 

Google Scholar 
Gerber, G. H. Evolution of the methods of spermatophore formation in pterygotan insects. Can. Entomol. 102, 358–362 (1970).Article 

Google Scholar 
Dallai, R., Afzelius, B. A., Lanzavecchia, S. & Bellon, P. L. Bizarre flagellum of thrips spermatozoa (Thysanoptera, Insecta). J. Morphol. 209, 343–347. https://doi.org/10.1002/jmor.1052090309 (1991).Article 
CAS 
PubMed 

Google Scholar 
Paccagnini, E., Mencarelli, C., Mercati, D., Afzelius, B. A. & Dallai, R. Ultrastructural analysis of the aberrant axoneme morphogenesis in thrips (Thysanoptera, Insecta). Cell Mot. Cytoskel. 64, 645–661. https://doi.org/10.1002/cm.20212 (2007).Article 

Google Scholar 
Paccagnini, E., Lupetti, P., Afzelius, B. A. & Dallai, R. New findings on sperm ultrastructure in thrips (Thysanoptera, Insecta). Arthropod. Struct. Dev. 38, 70–83. https://doi.org/10.1016/j.asd.2008.07.004 (2009).Article 
PubMed 

Google Scholar 
Paccagnini, E., Mercati, D., Giusti, F., Conti, B. & Dallai, R. The spermatogenesis and the sperm structure of Terebrantia (Thysanoptera, Insecta). Tissue & cell 42, 247–258. https://doi.org/10.1016/j.tice.2010.04.008 (2010).Article 

Google Scholar 
Pitnick, S., Wolfner, M. F. & Dorus, S. Post-ejaculatory modifications to sperm (PEMS). Biol. Rev. Camb. Philos. Soc. 95, 365–392. https://doi.org/10.1111/brv.12569 (2020).Article 
PubMed 

Google Scholar 
Karr, T. L., Swanson, W. J. & Snook, R. R. The evolutionary significance of variation in sperm–egg interactions. In Sperm biology. An evolutionary perspective, edited by T. R. Birkhead. 1st ed. (Academic Press/Elsevier, Amsterdam, 2009), pp. 305–365.Friedländer, M., Jeshtadi, A. & Reynolds, S. E. The structural mechanism of trypsin-induced intrinsic motility in Manduca sexta spermatozoa in vitro. J. Insect Physiol. 47, 245–255. https://doi.org/10.1016/s0022-1910(00)00109-8 (2001).Article 
PubMed 

Google Scholar 
Hughes, M. & Davey, K. G. The activity of spermatozoa of Periplaneta. J. Insect Physiol. 15, 1607–1616. https://doi.org/10.1016/0022-1910(69)90181-4 (1969).Article 

Google Scholar 
Longo, G. et al. Ultrastructural changes in sperm of Eyprepocnemis plorans (Charpentier) (Orthoptera: Acrididae) during storage of gametes in female genital tract. Inverteb. Reprod. Dev. 24, 1–6. https://doi.org/10.1080/07924259.1993.9672325 (1993).Article 

Google Scholar 
Makielski, S. K. The structure and maturation of the spermatozoa of Sciara coprophila. J. Morphol. 118, 11–41. https://doi.org/10.1002/jmor.1051180103 (1966).Article 
CAS 
PubMed 

Google Scholar 
Giuffrida, A. & Rosati, F. Changes in sperm tail of Eyprepocnemis plorans (Insects, Orthoptera) as a result of in vitro incubation in spermathecal extract. Inverteb. Reprod. Dev. 24, 47–52. https://doi.org/10.1080/07924259.1993.9672330 (1993).Article 

Google Scholar 
Giuffrida, A., Focarelli, R., Lampariello, R., Thole, H. & Rosati, F. Purification and properties of a 35 kDa glycoprotein from spermathecal extract of Eyprepocnemis plorans (Insecta, Orthoptera) with axonemal cytoskeleton disassembly activity. Insect Biochem. Mol. Biol. 26, 347–354. https://doi.org/10.1016/0965-1748(95)00095-X (1996).Article 
CAS 

Google Scholar 
Arakaki, N., Noda, H. & Yamagishi, K. Wolbachia-induced parthenogenesis in the egg parasitoid Telenomus nawai. Entomol. Exp. Appl. 96, 177–184. https://doi.org/10.1046/j.1570-7458.2000.00693.x (2000).Article 

Google Scholar 
Pannebakker, B. A. et al. Sexual functionality of Leptopilina clavipes (Hymenoptera: Figitidae) after reversing Wolbachia-induced parthenogenesis. J. Evol. Biol. 18, 1019–1028. https://doi.org/10.1111/j.1420-9101.2005.00898.x (2005).Article 
CAS 
PubMed 

Google Scholar 
Stouthamer, R., Russell, J. E., Vavre, F. & Nunney, L. Intragenomic conflict in populations infected by Parthenogenesis Inducing Wolbachia ends with irreversible loss of sexual reproduction. BMC Evol. Biol. 10, 229. https://doi.org/10.1186/1471-2148-10-229 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
Sänger, K. & Helfert, B. Comparative studies on number and position of the micropyles and the shape of the eggs of Saga pedo, S. natoliae and S. ephippigera (Orthoptera: Tettigoniidae). Entomologia 19, 49–56. https://doi.org/10.1127/ENTOM.GEN/19/1994/049 (1994).Article 

Google Scholar 
Gottlieb, Y. & Zchori-Fein, E. Irreversible thelytokous reproduction in Muscidifurax uniraptor. Entomol. Exp. Appl. 100, 271–278. https://doi.org/10.1046/j.1570-7458.2001.00874.x (2001).Article 

Google Scholar 
Stouthamer, R., Breeuwer, J. A. J. & Hurst, G. D. D. WOLBACHIA PIPIENTIS. Microbial manipulator of arthropod reproduction. Annu. Rev. Microbiol. 53, 71–102. https://doi.org/10.1146/annurev.micro.53.1.71 (1999).Article 
CAS 
PubMed 

Google Scholar 
Schwander, T., Crespi, B. J., Gries, R. & Gries, G. Neutral and selection-driven decay of sexual traits in asexual stick insects. Proc. Biol. Sci. 280, 20130823. https://doi.org/10.1098/rspb.2013.0823 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Werren, J. H., Baldo, L. & Clark, M. E. Wolbachia: master manipulators of invertebrate biology. Nat. Rev. Micro. 6, 741–751. https://doi.org/10.1038/nrmicro1969 (2008).Article 
CAS 

Google Scholar  More

Yang, X., Quam, M. B., Zhang, T. & Sang, S. Global burden for dengue and the evolving pattern in the past 30 years. J. Travel Med. 28, taab146 (2021).Simmons, C. P., Farrar, J. J., van Vinh Chau, N. & Wills, B. Dengue. N. Engl. J. Med. 366, 1423–1432 (2012).Article 
CAS 
PubMed 

Google Scholar 
Brady, O. J. et al. Refining the global spatial limits of dengue virus transmission by evidence-based consensus. PLoS Negl. Trop. Dis. 6, e1760 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Morin, C. W., Comrie, A. C. & Ernst, K. Climate and dengue transmission: Evidence and implications. Environ. Health Perspect. 121, 1264–1272 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Brady, O. J. et al. Global temperature constraints on Aedes aegypti and Ae. albopictus persistence and competence for dengue virus transmission. Parasites Vectors 7, 338 (2014).Betanzos-Reyes, Á. F. et al. Association of dengue fever with Aedes spp. abundance and climatological effects. Salud Pública de México 60, 12 (2017).WHO. Dengue and severe dengue. (2022).Gubler, D. J. Dengue and dengue hemorrhagic fever. Clin. Microbiol. Rev. 11, 480–496 (1998).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Delatte, H. et al. Aedes albopictus, vecteur des virus du chikungunya et de la dengue à La Réunion : biologie et contrôle. Parasite 15, 3–13 (2008).Article 
CAS 
PubMed 

Google Scholar 
Kles, V., Michault, A., Rodhain, F., Mevel, F. & Chastel, C. A serological survey regarding Flaviviridae infections on the island of Reunion (1971–1989). Bull. Soc. Pathol. Exot. 1990(87), 71–76 (1994).
Google Scholar 
Pierre, V. et al. Epidémie de dengue 1 à la Réunion en 2004. Journal de Veille Sanitaire (2005).Vincent, M. et al. From the threat to the large outbreak: dengue on Reunion Island, 2015 to 2018. Eurosurveillance 24, (2019).Cellule Santé Publique France en Région, ARS. Situation de la dengue à La Réunion au 15 décembre 2020. https://www.lareunion.ars.sante.fr/avec-le-retour-de-lete-agissons-des-maintenant-contre-la-dengue (2020).Agence Régionale de Santé. Communiqué de presse: dengue à La Réunion. Situation au 28 juillet 2021. https://www.lareunion.ars.sante.fr/system/files/2021-07/2021-07-28-Dengue-Situation à La Réunion_0.pdf (2021).Hafsia, S. et al. Overview of dengue outbreaks in the southwestern Indian Ocean and analysis of factors involved in the shift toward endemicity in Reunion Island: A systematic review. PLoS Negl. Trop. Dis. 16, e0010547 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Reiter, P., Fontenille, D. & Paupy, C. Aedes albopictus as an epidemic vector of chikungunya virus: Another emerging problem?. Lancet. Infect. Dis 6, 463–464 (2006).Article 
PubMed 

Google Scholar 
Njenga, M. K. et al. Tracking epidemic Chikungunya virus into the Indian Ocean from East Africa. J. Gen. Virol. 89, 2754–2760 (2008).Article 
CAS 

Google Scholar 
Soumahoro, M.-K. et al. The Chikungunya epidemic on La Réunion Island in 2005–2006: a cost-of-illness study. PLoS Negl. Trop. Dis. 5, e1197 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Larrieu, S., Balleydier, E., Renault, P., Baville, M. & Filleul, L. [Epidemiological surveillance du chikungunya on Reunion Island from 2005 to 2011]. Médecine tropicale : Revue du Corps de Santé colonial 72 Spec No, 38–42 (2012).Soghigian, J. et al. Genetic evidence for the origin of Aedes aegypti, the yellow fever mosquito, in the southwestern Indian Ocean. Mol. Ecol. 29, 3593–3606 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Kotsakiozi, P. et al. Population structure of a vector of human diseases: Aedes aegypti in its ancestral range Africa. Ecol. Evol. 8, 7835–7848 (2018).Article 
PubMed 

Google Scholar 
MacGregor, M. E. Aedes (Stegomyia) mascarensis, MacGregor: A new Mosquito from Mauritius. Bull. Entomol. Res. 14, 409–412 (1924).Article 

Google Scholar 
Salvan, M. & Mouchet, J. Aedes albopictus et Aedes aegypti à l’Ile de La Réunion. Ann. Soc. Belg. Med. Trop. 74, 323–326 (1994).CAS 
PubMed 

Google Scholar 
Bagny, L., Delatte, H., Quilici, S. & Fontenille, D. Progressive Decrease in Aedes aegypti distribution in Reunion Island since the 1900s. J. Med. Entomol. 46, 1541–1545 (2009).Article 
PubMed 

Google Scholar 
Le Vassal, J. J. paludisme à l’Ile de La Réunion. Per Gli Stud Della Maria 8, 18–27 (1907).
Google Scholar 
Delatte, H. et al. Geographic distribution and developmental sites of Aedes albopictus (Diptera: Culicidae) during a chikungunya epidemic event. Vector-Borne Zoon. Dis. 8, 25–34 (2008).Article 
CAS 

Google Scholar 
Hamon, J. Etudes biologique et systématique des Culicinae de l’Ile de La Réunion. Mem. Inst. Scient. Madagascar 4, 521–541 (1953).
Google Scholar 
Bouyer, J. & Lefrançois, T. Boosting the sterile insect technique to control mosquitoes. Trends Parasitol. 30, 271–273 (2014).Article 
PubMed 

Google Scholar 
Haramboure, M. et al. Modelling the control of Aedes albopictus mosquitoes based on sterile males release techniques in a tropical environment. Ecol. Model. 424, 109002 (2020).Article 

Google Scholar 
Bouyer, J., Yamada, H., Pereira, R., Bourtzis, K. & Vreysen, M. J. B. Phased conditional approach for mosquito management using sterile insect technique. Trends Parasitol. 36, 325–336 (2020).Article 
PubMed 

Google Scholar 
Bouyer, J. & Vreysen, M. J. B. Yes, irradiated sterile male mosquitoes can be sexually competitive!. Trends Parasitol. 36, 877–880 (2020).Article 
CAS 
PubMed 

Google Scholar 
Organization, W. H. Guidelines for laboratory and field testing of mosquito larvicides. WHO/CDS/WHOPES/GCDPP/2005.13 (2005).World Health Organization and Special Programme for Research and Training in Tropical Diseases and World Health Organization. Department of Control of Neglected Tropical Diseases and World Health Organization. Epidemic and Pandemic Alert. Dengue: Guidelines for diagnosis, treatment, prevention and control. (World Health Organization, 2009).Yap, H. H. Preliminary report on the color preference for oviposition by Aedes albopictus (Skuse) in the field. Southeast Asian J. Trop. Med. Public Health 6, 1–2 (1975).
Google Scholar 
Yap, H. H., Lee, C. Y., Chong, N. L., Foo, A. E. S. & Lim, M. P. Oviposition site preference of Aedes albopictus in the laboratory. J. Am. Mosquito Control Assoc. Mosquito News 11, 128–132 (1995).CAS 
PubMed 

Google Scholar 
Marin, G., Mahiba, B., Arivoli, S. & Tennyson, S. Does colour of ovitrap influence the ovipositional preference of Aedes aegypti Linnaeus 1762 (Diptera: Culicidae). Int. J. Mosq. Res 7, 11–15 (2020).CAS 

Google Scholar 
Claudel, I. et al. To bait or not to bait? Optimizing the use of adult mosquito traps for monitoring arbovirus vector populations in La Réunion Island. (2022). https://doi.org/10.21203/rs.3.rs-1798972/v1.Cleveland, W. S. Visualizing data. (Hobart press, 1993).Lamigueiro, Ó. P. Displaying time series, spatial, and space-time data with R. (Chapman; Hall/CRC, 2018).Yoshioka, M. et al. Diet and density dependent competition affect larval performance and oviposition site selection in the mosquito species Aedes albopictus (Diptera: Culicidae). Parasites Vectors 5, (2012).Wong, J., Stoddard, S. T., Astete, H., Morrison, A. C. & Scott, T. W. Oviposition site selection by the dengue vector Aedes aegypti and its implications for dengue control. PLoS Negl. Trop. Dis. 5, e1015 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Hawley, W. A. The biology of aedes albopictus. J. Am. Mosquito Control Assoc. Suppl 1, 1–39 (1988).Delatte, H., Gimonneau, G., Triboire, A. & Fontenille, D. Influence of temperature on immature development, survival, longevity, fecundity, and gonotrophic cycles of Aedes albopictus, vector of Chikungunya and dengue in the Indian Ocean. J. Med. Entomol. 46, 33–41 (2009).Article 
CAS 
PubMed 

Google Scholar 
Daugherty, M. P., Alto, B. W. & Juliano, S. A. Invertebrate carcasses as a resource for competing Aedes albopictus and Aedes aegypti (Diptera: Culicidae). J. Med. Entomol. 37, 364–372 (2000).Article 
CAS 
PubMed 

Google Scholar 
Papaj, D. R. & Rausher, M. D. Individual variation in host location by phytophagous insects. Herbivorous Insects: Host seeking behavior and mechanisms 77–127 (1983).Valladares, G. & Lawton, J. H. Host-plant selection in the holly leaf-miner: Does mother know best?. J. Anim. Ecol. 60, 227 (1991).Article 

Google Scholar 
Ellis, A. M. Incorporating density dependence into the oviposition preference-offspring performance hypothesis. J. Anim. Ecol. 77, 247–256 (2008).Article 
PubMed 

Google Scholar 
Juliano, S. A., OMeara, G. F., Morrill, J. R. & Cutwa, M. M. Desiccation and thermal tolerance of eggs and the coexistence of competing mosquitoes. Oecologia 130, 458–469 (2002).Costanzo, K. S., Kesavaraju, B. & Juliano, S. A. Condition-specific competion in container mosquitoes: The role of non-competing life-history stages. Ecology 86, 3289–3295 (2005).Article 
PubMed 

Google Scholar 
Sanchez, M. & Probst, J.-M. Distribution and conservation status of the Manapany day gecko, Phelsuma inexpectata MERTENS, 1966, an endemic threatened reptile from Réunion Island (Squamata: Gekkonidae). Cahiers scientifiques de l’océan Indien occidental 2, (2011).Braks, M. A. H., Honório, N. A., Lounibos, L. P., De-Oliveira, R. L. & Juliano, S. A. Interspecific competition between two invasive species of container mosquitoes, Aedes aegypti and Aedes albopictus (Diptera: Culicidae), in Brazil. Ann. Entomol. Soc. Am. 97, 130–139 (2004).Article 

Google Scholar 
Moore, C. G. & Fisher, B. R. Competition in mosquitoes.1 Density and species ratio effects on growth, mortality, fecundity, and production of growth retardant2. Ann. Entomol. Soc. Am. 62, 1325–1331 (1969).Madeira, N. G., Macharelli, C. A. & Carvalho, L. R. Variation of the Oviposition Preferences of Aedes aegypti in Function of Substratum and Humidity. Mem. Inst. Oswaldo Cruz 97, 415–420 (2002).Article 
CAS 
PubMed 

Google Scholar 
Bellini, R. et al. Use of the sterile insect technique against Aedes albopictus in Italy: first results of a pilot trial. in Area-wide control of insect pests 505–515 (Springer, 2007).Boussès, P., Dehecq, J. S., Brengues, C. & Fontenille, D. Inventaire actualisé des moustiques (Diptera : Culicidae) de l’île de La Réunion, océan Indien. Bulletin de la Société de pathologie exotique 106, 113–125 (2013).Article 
PubMed 

Google Scholar 
Zeileis, A., Kleiber, C. & Jackman, S. Regression models for count data in R. J. Stat. Softw. 27, (2008).Sileshi, G. Selecting the right statistical model for analysis of insect count data by using information theoretic measures. Bull. Entomol. Res. 96, 479–488 (2006).CAS 
PubMed 

Google Scholar 
Guthery, F. S., Burnham, K. P. & Anderson, D. R. Model selection and multimodel inference: A practical information-theoretic approach. J. Wildl. Manag. 67, 655 (2003).Article 

Google Scholar 
Hurvich, C. M. & Tsai, C.-L. Model selection for extended quasi-likelihood models in small samples. Biometrics 51, 1077–1084 (1995).Article 
CAS 
PubMed 
MATH 

Google Scholar 
Burnham, K. P. & Anderson, D. R. Model selection and multimodel inference: a practical information-theoretic approach. 496 (Springer-Verlag, 2002).Manly, B. F. J. Randomization, bootstrap and Monte Carlo methods in biology. 399 (CRC Press / Chapman & Hall, 2006). https://doi.org/10.1201/9781315273075.Holm, S. A simple sequentially rejective multiple test procedure. Scand. J. Stat. 65–70 (1979).R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2022).Lesnoff, M. & Lancelot, R. aods3: analysis of overdispersed data using S3 methods. (2018).Barton, K. MuMIn: Multi-Model Inference. (2022).Xie, Y., Dervieux, C. & Riederer, E. R Markdown Cookbook. (Chapman; Hall/CRC, 2020). https://doi.org/10.1201/9781003097471. More

NEWS AND VIEWS
02 November 2022

Predicting the risk of extinction from climate change requires an understanding of the interactions between species. An analysis of how changes in rainfall affect competition between plant species offers a way of tackling this challenge. More

Study site and future climate treatmentOur study site is located at the Bad Lauchstädt Field Research Station, Bad Lauchstädt, Germany (51° 22060 N, 11° 50060 E), which belongs to the Helmholtz Centre for Environmental Research–UFZ. Long-term mean annual precipitation in the area is 489 mm and the mean annual temperature is 8.9 °C (ref. 32). During 2018 and 2019, Europe experienced a record-setting drought that was especially severe in 2018 (refs. 33,34); the mean annual precipitation at our study site in 2018 and 2019 was 254 mm and 353 mm, respectively, whereas 2017 was a more normal year, with a mean annual precipitation of 403 mm. Mean annual temperatures were above average: 2017, 10.5 °C; 2018, 10.8 °C; 2019, 11.2 °C (data from the weather station at the Bad Lauchstädt field station). The soils in the study area are fertile Haplic Chernozem type32,35.Our eDiValo experiment was conducted in the GCEF, which was designed to investigate climate change effects under different land-use scenarios32. We used 10 ‘extensively’ used pastures of the GCEF in our experiment; that is, 384-m2 (16 × 24 m) areas of grassland (hereafter called ‘pastures’) that were grazed by a flock of 20 sheep 2–3 times each year. Grazing was implemented as short-time high-intensity grazing events, each lasting 24 h (ref. 32). This type of high-intensity but short-term grazing is considered better in maintaining species richness as it gives plants more time to recover between grazing events36. It is also a recommended management type for nature conservation areas in Germany37. Vegetation in the pastures was species-rich grassland vegetation that is typical of drier regions of central Germany32,38. The whole GCEF was fenced to exclude native large mammalian herbivores (for example, deer); however, European hare (Lepus europaeus), wood mice (Apodemus sylvaticus) and voles (Microtus arvalis) are common at the site.Our experimental design was originally intended to test the dependence of light competition on nutrient and herbivory under current and future climatic scenarios. Although we included both climate treatments in our data, climate was never significant for richness and Shannon diversity, either alone or in interaction with other factors, and our focus was therefore on the other treatments. Five of the above random pastures received future climatic treatment which was based on different dynamic regional climate models for Germany, all predicting an increased mean temperature by approximately 2 °C year-round, strongly decreased summer precipitation and slightly increased spring and autumn precipitation (https://www.regionaler-klimaatlas.de/) (ref. 32). Passive night-time (after sunset and before sunrise) warming through the use of roller blinds attached to the GCEF roof and eastern and western wall structures was used to increase the air temperature. In each spring (1 March–31 May) and autumn (1 September–30 November), future climate plots received 110% of the ambient rainfall and in the summer (1 June–31 August), they received 80% of the ambient rainfall. The precipitation treatment was adjusted weekly and compensated for a possible night-time reduction in rainfall due to temperature treatment. A detailed description of the future climate treatment is provided in a previous report32.Fertilization, herbivore exclusion and light additionWe first tested whether adding light can offset the negative effect of fertilization on plant diversity. In May 2017, we established a full-factorial experiment of fertilization and light addition. Within each 10 pastures (5 in ambient climatic conditions, 5 in future climatic conditions), we established 4 plots of 1.4 × 1.4 m, separated by a 1-m buffer zone (hereafter called ‘blocks’), in total 40 plots and 10 blocks. At the time the experiment was established, vegetation in the whole experimental area (that is, in a block of 4 plots and the surrounding 1-m area) was trimmed to a height of 5 cm to make conditions uniform and the whole area was temporarily fenced to let the experiment establish and fertilization effects develop. The temporary fence was removed in August when the herbivore exclusion treatment was started. Therefore, there was no grazing by sheep in the experimental plots in the summer of 2017. Two randomly chosen plots received fertilizer treatment and two were controls. For the former (fertilizer-treatment plots), slow-release granular NPK fertilizer (a mixture of Haifa Multicote 2 M 40-0-0 40% N; Triple Super Phosphate (TSP) 45% P205; and potassium sulfate fertilizer 50% K2O, 45% SO3) was added twice per growing season, in a total of 10 g N, 10 g P and 10 g K per m² (see ref. 3 for a similar protocol that is used in grasslands worldwide). In 2017, the first fertilization was done at the beginning of June right after establishing the experiment and the second fertilization was done at the beginning of July. In the subsequent years, the first fertilization was done at the beginning of the growing season (late March–April) and the second fertilization was done in June. In 2019, two previously unfertilized plots were accidentally fertilized and were thereafter treated as fertilized plots. To manipulate light, 1.4 × 1.4-m plots were further divided into two subplots, 0.7 m × 1.4 m each, and one of these was randomly assigned to the light-addition treatment, resulting in 80 subplots (Fig. 1). We installed two 120-cm-long and 3.5-cm-wide recently developed LED lamps (C65, Valoya) parallel to each other and at a 28-cm distance from each other to each light-addition subplot. To increase light for the small understory plants that are the most likely to suffer from competition for light, we installed the lamps 10 cm above the smallest plants. The lamps were gradually uplifted over the course of the growing season to follow the growth of the smallest plants. As our light-addition treatment was intended to mimic natural sunlight (that is, making a gap in a dense vegetation and allowing the sunshine in), we chose the spectrum of the lamps to include all wavelengths of sunlight, including small amounts of ultraviolet and infrared. Each lamp added roughly 350–400 µmol and did not alter the air or aboveground soil surface temperature (Fig. 1b), which is an improvement on previous studies12. Each year, we added light during the active growing season: the lamps were switched on early in the spring (March–April), when temperatures were clearly above zero, and switched off and removed when temperatures dropped close to zero in November–December and aboveground plant parts had died and formed litter. Each day, the lamps were set to switch on two hours after sunrise, and to switch off two hours before sunset, and when the temperature exceeded 28 °C to prevent overheating. We did not install unpowered lamps to unlighted plots because our modern, narrow LED lamps caused minimal disturbance (see below) and no heating (Fig. 1b), and because unpowered lamps would have added an artefact in that they create shade that does not occur when the lamps are on in lighted plots.At the end of August 2017, after running the fertilization–light-addition experiment for one growing season, we expanded the experiment by implementing the herbivore exclusion treatment in a full-factorial combination with the other treatments. Two of the previously established 1.4 m × 1.4-m plots, one with and one without the fertilization treatment, were randomly allotted to the herbivore (sheep) exclusion treatment and fenced with rectangular metal fences of 1.8 m × 1.8 m, 82 cm height and 10 cm mesh size. At the same time, the temporary fence established in May 2017 was removed from around the whole experimental area, allowing the grazing of sheep in unfenced plots. The fences did not exclude mice, voles and hares. For the time of each grazing event, lamps in grazed subplots were removed and switched off in the ungrazed subplots. Uplifting the lamps from grazed plots did not cause disturbance because vegetation in grazed plots was always short and did not reach above the lamps. Inside exclosures, lamps were always kept in place during the growing season, and plants could freely grow around and above them.Plant community and trait samplingIn July 2017, we established 50 cm × 50-cm permanent quadrats in every subplot for plant community sampling. We visually estimated the per cent areal cover for all species occurring in the quadrats, and litter cover, from the beginning of June to mid-June 2019, when the vegetation was at its peak biomass. The 2017 sampling happened later, in mid-July, because vegetation in all plots and surrounding areas was trimmed to a height of 5 cm at the time of the establishment of the experiment at the end of May, and it took later for vegetation to reach its peak biomass. In 2018, the effects of drought were devastating, and most plants had senesced or died before the planned sampling date; we therefore omitted the year 2018. At the beginning of each growing season—that is, when the lamps were installed and switched on—there was very little live biomass in the plots, and the maximum height of existing plants was approximately 5 cm (in all plots). During the peak biomass the maximum plant height was up to approximately 1 m; however, it varied greatly between the treatments and was especially low in grazed plots. All vegetation surveys were done by the same trained and experienced person with a minimum estimate threshold of 0.1%. We used plant cover data to calculate species richness and Shannon diversity.In May–June 2020, we measured plant height (centimetres), SLA (leaf area in square millimetres per milligram of dry mass), foliar C:N (based on the per cent C and N in plant leaves) and LWC (leaf water content as 1,000 − LDMC (the ratio of leaf dry mass to saturated fresh mass), expressed as milligrams per gram39) for most species occurring in the experimental plots, and complemented the trait data from the TRY Plant Trait Database40,41,42 (v.5.0; https://www.try-db.org/TryWeb/Home.php) and for one species one trait value from another source9. The trait data were collected from seven to ten individuals per species from the study site or close areas; the collection and handling followed standard protocols39. We chose these traits because they are widely documented to be associated with responsiveness to soil nutrients, herbivory and light9,26,27,43,44,45,46. We used all traits as, although they partially reflect similar ecological adaptations (for example, leaf economics spectrum43), they could also potentially reflect independent and distinctive processes, and differently mediate the responses of species to our treatments. For example, SLA and LWC in our dataset correlated weakly (r2 = 0.16), but were to a greater extent uncorrelated (Extended Data Table 6), and could function differently, for example, in light capture and drought tolerance26,39. In 2017, our trait data covered on average 97.7–98.6% of the total cover in the plots, the value slightly differing depending on the trait as we did not have all traits for all species. Our own trait collections covered on average 96.6–97.6% and TRY data covered on average 0.9–2% of the total cover. In 2019, the whole trait data covered on average 99.5% of the total cover in the plots, again slightly depending on the trait. Our own trait collections covered on average 94.2–96.5% and TRY data covered on average 2.7–5.3% of the total cover.Abiotic environmental measurementsWe measured several soil and other environmental properties from the experimental plots. Light availability (photosynthetically active radiation; PAR) in unlighted and lighted (under lamps) subplots was measured using LI-190R and LI-250A meters (LI-COR), approximately 7–10 cm under the lamps and 15–20 cm above ground level. We measured light availability from the same distance to the ground in unlighted plots. Measurements of light availability were done in mid-July 2020 on three consecutive cloudless days around noon. Note that in grazed plots, light levels between lighted and unlighted plots are more similar than inside exclosures (Fig. 1), because herbivores keep the vegetation short, and natural sunlight can therefore reach under the lamps where the light measurements were taken. Air temperature and humidity were recorded from unlighted and lighted (under lamps) subplots using loggers (HOBO MX2301A, Onset Computer Cooperation) that were installed approximately 7 cm under the lamps and to the same height from the ground in unlighted plots, and were replicated under different combinations of fertilization, herbivore exclusion and light addition in ambient climatic conditions three times (n = 3). The logger data were collected in May 2019 before the effects of drought were visible.Statistical analysisWe analysed our data in two steps. First, to test whether competition for light mediates the effect of fertilization on diversity, we analysed the effects of fertilization and light and their interaction on species richness and Shannon diversity using data from 2017, when the herbivore exclusion treatment had not yet been implemented. We also analysed the effects of treatment on total vegetation cover and litter cover. We fit LME models in which diversity (species richness and Shannon diversity), total cover and litter cover, each in their own model, were explained by fertilization, light addition and their interaction (fixed variables). All treatments were categorical variables with two levels (treated and untreated). In each model, subplot was nested within plot, which was nested within block (nested random variable). We simplified the models using the anova() function for model comparison in the nlme and lme4 packages in R (ref. 47) (on the basis of log likelihood ratio tests; P ≥ 0.05; Extended Data Table 2). This was done to uncover the significance of the main effects and interaction terms, to avoid overparametrization47,48 and to provide model-derived parameter estimates for the figures (Extended Data Table 5). However, we also provide full model results that are qualitatively similar to the results of simplified models (Extended Data Tables 3 and 4); therefore, model choice did not affect our conclusions. Climate treatment was included in all original models but was never significant for richness and diversity, and was not considered further. Total cover and litter results for 2017 are reported in Extended Data Figs. 1a,b and 3a). As there was heterogeneity in the variance structure between treatments, we used the varIdent() function in the nlme package in R to allow each treatment combination to have a different variance. Model fit was inspected using model diagnostic plots in the package nlme. In the full design with climate included, the number of replicates per treatment combination was ten.Second, to include herbivore exclusion to the experimental design and to test whether competition for light mediates the effect of herbivore exclusion on diversity, and whether competition for light, herbivory and fertilization interact, we analysed the effects of herbivore exclusion, fertilization, light and their interactions on species richness and Shannon diversity using data from 2019. All treatments were categorical variables with two levels (treated and untreated). We also analysed the effects of treatment on total vegetation cover and litter cover. We fit similar models to those described above, except that herbivore exclusion was an additional fixed factor in the models. We simplified the models, used the varIdent() function to account for heteroscedasticity and checked the model fit using model diagnostic plots, as above. Climate treatment was included in all original models but was significant for litter cover only, and was not considered further. In the full design with climate included, the number of replicates per treatment combination was five.To further assess which plant traits increased the probability of species benefiting from the addition of light, we first created a binary response variable: those species that increased from unlighted to lighted plots (that is, had a higher value in a lighted than an unlighted plot) were given a value of 1 and those that did not were given a value of 0. This response variable takes into account rare species that emerged or persisted in the lighted plots but were absent in the unlighted plots (that is, species gains and losses) and changes in small, subordinate species (those that are likely to benefit from light addition) with small but consistently trait-dependent changes in response to light. It is also in line with our species richness analyses, as species gains and losses ultimately determine richness responses. We did not use different indexes (for example, lnRR or RII) because these could not handle multiple zero values and species losses or gains (that is, species having zero cover in either unlighted or lighted subplots). Second, we fit GLME models with a binomial error structure (family = “binomial”, link = “logit”) in which a probability of a species increasing from unlighted to lighted plots was explained by categorical experimental treatments (fertilization, herbivore exclusion and their interactions), traits (SLA, height, LWC, foliar C:N), and interactions between the treatments and traits. Each trait was analysed in its own model as some of the traits were correlated (Extended Data Table 6), and to avoid overly complex models and overparametrization47,48. We included all species for which we had traits in the models. As we calculated the increase in cover from unlighted to lighted plots, our smallest experimental unit in trait analyses was a plot (not a subplot, unlike in other analyses). As there were several species in the same plots, we nested species within plots, and plots within blocks. We similarly simplified the models to include only significant variables (on the basis of χ2 tests; P ≥ 0.05). We did not include a crossed random effect for species in the models because the full models with a more complex random structure did not converge; however, when we refitted the simplified models with a crossed random effect for species, we found that the models converged (with scaled data) and that the significance of the effects remained qualitatively the same. Climate was included in all original models but was never significant. In addition, C:N and height did not predict the responsiveness of species to light in either year (P ≥ 0.13 for both); results are therefore not shown. In the full design with climate included, the number of replicates per treatment combination was five; however, the number of observations was greater (see Fig. 4 and Extended Data Fig. 4). To make sure that our results for SLA and LWC were not influenced by whether they were analysed in separate models or in the same model, or by the order in which they were in the models, we also performed analyses in which both SLA and LWC were included (in both orders). Results remained qualitatively similar and are not discussed further.Furthermore, to check whether our trait results were driven primarily by species gains and losses or changes in abundance, we ran additional trait analyses for which we calculated the change in cover between lighted and unlighted subplots (cover in lighted subplot − cover in unlighted subplot), and analysed the ‘change’ with otherwise similar trait models to those described above, except that we used Gaussian error structure. With this index, which gives a disproportionate importance to the abundant species, we found that traits were poor predictors of changes in cover between lighted and unlighted plots (all interactions were non-significant, P  > 0.05, except for a marginally significant C:N × fertilization interaction in 2017 that was no longer visible in 2019; results not shown; codes and data available in the Dryad repository). We also analysed presence–absence-based species losses and gains. In these models, each species was given a value of 1 when it was present in the lighted subplot but absent from the unlighted subplot; otherwise, these models were similar to the binomial trait models described above. These models produced, to a large extent, similar results to our models using the probability of increase in response to light as a response variable (results not shown; codes and data available in the Dryad repository). These additional analyses and results support using the probability of increase in response to light as our response variable, rather than abundance-based metrics, as it includes both gains and losses and abundance aspects, and is therefore a general test that is well suited to assessing species gains and extinctions and changes in subordinate species.All statistical analyses were performed using R v. 4.0.0 (ref. 49). We used the nlme package (v.3.1.147) for LME models50, the lme4 package (v.1.1.23) for GLME models51, and the car package52 for P values (v.3.07).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Schimper, A. F. W. Plant Geography upon a Physiological Basis (Clarendon Press, 1903).Alexander, J. M., Diez, J. M. & Levine, J. M. Novel competitors shape species’ responses to climate change. Nature 525, 515–518 (2015).Article 
ADS 
CAS 
PubMed 

Google Scholar 
HilleRisLambers, J., Harsch, M. A., Ettinger, A. K., Ford, K. R. & Theobald, E. J. How will biotic interactions influence climate change-induced range shifts? Ann. N. Y. Acad. Sci. 1297, 112–125 (2013).PubMed 

Google Scholar 
Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366 (2000).Article 

Google Scholar 
Loarie, S. R., Weiss, S. B., Hamilton, H., Branciforte, R. & Kraft, N. J. B. The geography of climate change: implications for conservation biogeography. Divers. Distrib. 16, 476–487 (2010).Article 

Google Scholar 
Callaway, R. M. et al. Positive interactions among alpine plants increase with stress. Nature 417, 844–848 (2002).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Dybzinski, R. & Tilman, D. Resource use patterns predict long‐term outcomes of plant competition for nutrients and light. Am. Nat. 170, 305–318 (2007).Article 
PubMed 

Google Scholar 
Hautier, Y., Niklaus, P. A. & Hector, A. Competition for light causes plant biodiversity loss after eutrophication. Science 324, 636–638 (2009).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Levine, J. M. & HilleRisLambers, J. The importance of niches for the maintenance of species diversity. Nature 461, 254–257 (2009).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Kraft, N. J. B., Godoy, O. & Levine, J. M. Plant functional traits and the multidimensional nature of species coexistence. Proc. Natl Acad. Sci. USA 112, 797–802 (2015).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Knapp, A. K. et al. Rainfall variability, carbon cycling, and plant species diversity in a mesic grassland. Science 298, 2202–2205 (2002).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Sandel, B. et al. Contrasting trait responses in plant communities to experimental and geographic variation in precipitation. New Phytol. 188, 565–575 (2010).Article 
PubMed 

Google Scholar 
Esch, E. H., Ashbacher, A. C., Kopp, C. W. & Cleland, E. E. Competition reverses the response of shrub seedling mortality and growth along a soil moisture gradient. J. Ecol. 106, 2096–2108 (2018).Article 

Google Scholar 
Alon, M. & Sternberg, M. Effects of extreme drought on primary production, species composition and species diversity of a Mediterranean annual plant community. J. Veg. Sci. 30, 1045–1061 (2019).Article 

Google Scholar 
Chesson, P. Updates on mechanisms of maintenance of species diversity. J. Ecol. 106, 1773–1794 (2018).Article 

Google Scholar 
Barabás, G., D’Andrea, R. & Stump, S. M. Chesson’s coexistence theory. Ecol. Monogr. 88, 277–303 (2018).Article 

Google Scholar 
Ellner, S. P., Snyder, R. E., Adler, P. B. & Hooker, G. An expanded modern coexistence theory for empirical applications. Ecol. Lett. 22, 3–18 (2019).Article 
ADS 
PubMed 

Google Scholar 
Adler, P., Hillerislambers, J. & Levine, J. A niche for neutrality. Ecol. Lett. 10, 95–104 (2007).Article 
PubMed 

Google Scholar 
Germain, R. M., Mayfield, M. M. & Gilbert, B. The ‘filtering’ metaphor revisited: competition and environment jointly structure invasibility and coexistence. Biol. Lett. 14, 20180460 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Pau, S. et al. Predicting phenology by integrating ecology, evolution and climate science. Glob. Change Biol. 17, 3633–3643 (2011).Article 
ADS 

Google Scholar 
Fargione, J. & Tilman, D. Niche differences in phenology and rooting depth promote coexistence with a dominant C4 bunchgrass. Oecologia 143, 598–606 (2005).Article 
ADS 
PubMed 

Google Scholar 
Godoy, O., Kraft, N. J. B. & Levine, J. M. Phylogenetic relatedness and the determinants of competitive outcomes. Ecol. Lett. 17, 836–844 (2014).Article 
PubMed 

Google Scholar 
Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167–171 (2016).Article 
ADS 
PubMed 

Google Scholar 
Kunstler, G. et al. Plant functional traits have globally consistent effects on competition. Nature 529, 204–207 (2016).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Diffenbaugh, N. S., Swain, D. L. & Touma, D. Anthropogenic warming has increased drought risk in California. Proc. Natl Acad. Sci. USA 112, 3931–3936 (2015).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Swain, D. L., Langenbrunner, B., Neelin, J. D. & Hall, A. Increasing precipitation volatility in twenty-first-century California. Nat. Clim. Change 8, 427–433 (2018).Article 
ADS 

Google Scholar 
Chesson, P. Geometry, heterogeneity and competition in variable environments. Phil. Trans. R. Soc. Lond. B 330, 165–173 (1990).Article 
ADS 

Google Scholar 
Aronson, J., Kigel, J., Shmida, A. & Klein, J. Adaptive phenology of desert and Mediterranean populations of annual plants grown with and without water stress. Oecologia 89, 17–26 (1992).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Santa Barbara County Public Works water resources hydrology: historical rainfall data: daily and monthly rainfall. County of Santa Barbara http://www.countyofsb.org/pwd/water/downloads/hydro/421dailys.pdf (2019).Kandlikar, G. S., Kleinhesselink, A. R. & Kraft, N. J. B. Functional traits predict species responses to environmental variation in a California grassland annual plant community. J. Ecol. 110, 833–844 (2022).Article 
CAS 

Google Scholar 
Cleland, E. E. et al. Sensitivity of grassland plant community composition to spatial vs. temporal variation in precipitation. Ecology 94, 1687–1696 (2013).Article 
PubMed 

Google Scholar 
Usinowicz, J. et al. Temporal coexistence mechanisms contribute to the latitudinal gradient in forest diversity. Nature 550, 105–108 (2017).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Kandlikar, G. S., Johnson, C. A., Yan, X., Kraft, N. J. B. & Levine, J. M. Winning and losing with microbes: how microbially mediated fitness differences influence plant diversity. Ecol. Lett. 22, 1178–1191 (2019).PubMed 

Google Scholar 
Kleinhesselink, A. R., Kraft, N. J. B., Pacala, S. W. & Levine, J. M. Detecting and interpreting higher order interactions in ecological communities. Ecol. Lett. 25, 1604–1617 (2022).Article 
PubMed 

Google Scholar 
Saavedra, S. et al. A structural approach for understanding multispecies coexistence. Ecol. Monogr. 87, 470–486 (2017).Article 

Google Scholar 
Levine, J. I., Levine, J. M., Gibbs, T. & Pacala, S. W. Competition for water and species coexistence in phenologically structured annual plant communities. Ecol. Lett. 25, 1110–1125 (2022).Article 
PubMed 

Google Scholar 
Farrior, C. E. et al. Resource limitation in a competitive context determines complex plant responses to experimental resource additions. Ecology 94, 2505–2517 (2013).Article 
PubMed 

Google Scholar 
Harrison, S., Grace, J. B., Davies, K. F., Safford, H. D. & Viers, J. H. Invasion in a diversity hotspot: exotic cover and native richness in the Californian serpentine flora. Ecology 87, 695–703 (2006).Article 
PubMed 

Google Scholar 
Pérez-Harguindeguy, N. et al. New handbook for standardised measurement of plant functional traits worldwide. Aust. J. Bot. 61, 167–234 (2013).Article 

Google Scholar 
Godoy, O. & Levine, J. M. Phenology effects on invasion success: insights from coupling field experiments to coexistence theory. Ecology 95, 726–736 (2014).Article 
PubMed 

Google Scholar  More

All the results were reported relative to the control, unless specifically stated to the contrary or for clarity.Growth of cucumber plants in response to different biofertilizersSoilThere was no significant difference in cucumber growth before microbial fertilizer was applied. However, some microbial fertilizers significantly increased cucumber height and stem diameter when they were applied within 4 weeks from when the seedlings were planted (Fig. 1a,b,e,f). In the second week, SHZ and SMF increased plant height by 11.2 and 9.5%, respectively. In the third week, S267, SBS, SBH, SM and SHZ increased plant height by 12.0, 13.8, 15.0, 20.5 and 26.9%, respectively (Fig. 1a). In the fourth and fifth weeks, some treatments significantly increased cucumber height. In the second and third weeks, S267 significantly increased stem diameter by 21.2 and 16.8% (Fig. 1b).Figure 1Effect of different biofertilizer treatments on the growth of cucumber seedlings produced in soil or substrate in a greenhouse. S267 = Trichoderma Strain 267 added to soil; SBH = Bacillus subtilis and T. harzianum biofertilizers added to soil; SBS = B. subtilis biofertilizer added to the soil; SM = Compound biofertilizer added to soil; SHZ = T. harzianum biofertilizer added to soil; SCK = Untreated soil. US267 = T.267 biofertilizer added to substrate; USBH = B. subtilis and T. harzianum biofertilizers added to substrate; USBS = B. subtilis biofertilizer added to substrate; USM = Compound biofertilizer added to substrate; USHZ = T. harzianum biofertilizer added to substrate; USCK = Untreated substrate.Full size imageOver the subsequent 5 weeks, some microbial fertilizer treatments decreased cucumber height and stem diameter (Fig. 1g,h).SubstrateThere were no significant differences in cucumber growth before microbial fertilizer microbial fertilizer was applied (Fig. 1c,d,g,h). However, within 4 weeks of applying the microbial fertilizer, each biofertilizer treatment applied significantly increased cucumber height (Fig. 1c). US267 and USHZ significantly increased cucumber height by 39.8–75.4% and 56.1–86.1%, respectively. US267, USM and USHZ significantly increased the stem diameter by 76.8–108.9%, 71.1–97.6% and 80.4–122.4%, respectively (Fig. 1d).Over the subsequent 5 weeks, US267, USM and USHZ treatments continued to significantly increase cucumber height and stem diameter (Fig. 1g,h).Changes in the taxonomic composition of soil-borne fungal pathogensSoilBiofertilizers application significantly reduced the taxonomic composition of soil-borne fungal pathogens at different times during the cucumber growth period (Tables 1 and 2). Fusarium spp. were significantly reduced (T, 63.8% reduction, P  More