Philippa Kaur

Mobile phone data can be used to inform different aspects of COVID-19 response (Table 1). At the population level, quantifying changes in human mobility or clustering can help evaluate the impact of an NPI and identify hotspots where additional or different interventions may need to be applied. At the individual level, mobile phone data may be used to understand patterns of individual contacts and enhance contact tracing.
Table 1 Summary of types, metrics, and proposed applications of mobile phone data.
Full size table

Evaluating current interventions and monitoring their release
The most widely used application of mobile phone data in public health to date is the use of telecom geolocation data to track population movements11,12. Mobile phone operators routinely collect Call Detail Records (CDRs) that contain a timestamp and GPS location with a unique identifier for all subscribers. These data thus are typically readily available and offer high coverage to estimate mobility patterns of individuals using their mobile devices. We note that similar time-resolved GPS location data may be passively collected through certain applications, though typically for only a subset of subscribers that may introduce further bias.
CDRs can be used to generate a number of metrics for characterizing large, population-level mobility patterns. Origin-Destination (OD) matrices reflect the number of times a trip is made between two locations (of varying spatial resolution) in a certain period. These matrices can be analyzed over time to detect temporal trends (i.e., holidays, seasonality, weekday vs weekend) and regular hotspots of attraction. These spatial and temporal flows of individuals between locations, including the magnitude and frequency of these movements, can be used to understand the risk of importation from areas with ongoing outbreaks to areas without sustained transmission where there is a risk of reintroduction and resurgence. Aggregate flows can also be used to retrace the likely introduction and spread of an outbreak in new areas and to inform future projections of disease risk or burden across space and decision making around the design and implementation of travel restrictions or increased surveillance.
Aggregate mobility patterns may also be critical pieces of evidence when evaluating the effectiveness of various NPIs. Most NPIs are reliant on modifying physical behavior. Monitoring the volume, frequency, and average distance of flow during interventions can be used to directly quantify the adoption and effect of these interventions, and identify areas of high potential risk to target with different interventions. There are already identified associations between reductions in population-level mobility within and between different locations and COVID-19 incidence6,10,29, though further exploration of which population-level metrics are most closely related to changes in disease risk and whether these associations are sustained throughout an outbreak is needed30. These associations would ideally be interrogated to identify individual behaviors associated with mobility measures that are also associated with individual risk of COVID-19.
The effect on NPIs can also be monitored through subscriber density metrics that combine the recorded GPS location and timestamp of CDRs to capture the real-time population density and identify potential hotspots. When using finer-scale GPS location data, these density metrics may quantify the likelihood or frequency that users came into proximal contact. A third metric derived from CDR or GPS location data, the radius of gyration, quantifies the range over which a single person may travel in a specified time period. Importantly, the data required for these applications are non-identifiable; they cannot be used to identify any given individual’s interactions, but provide population-level insight into the average clustering and movement of individuals. These metrics, along with traditional OD matrix flows, were recently employed in Italy as a way to evaluate the impact of its national lockdown31. Traffic flow between provinces and probability of colocation were reduced initially in the northern provinces, where the COVID-19 outbreak was first observed, a clear signal of reactive social distancing. As the epidemic progressed, and especially once the national lockdown was enforced, the entire country saw a reduction in traffic between provinces; however, the probability of colocation remained highly dependent on province and was likely attributed to the number of cases reported in each province. Interestingly, the average distance traveled by individuals was significantly reduced across all provinces after the initial outbreak was confirmed.
The use of Bluetooth data (records of proximal interactions between Bluetooth-enabled devices) to quantify physical clustering or real-time density of subscribers at small spatial scales (e.g., zip codes) and fine temporal resolution has been explored for the purposes of contact tracing (see below). The use of these data has been considered less for population-level analyses, though it offers another source of information on behavioral changes under different NPIs. When activated, mobile phones will emit a Bluetooth beacon that is detected by other activated phones. When two Bluetooth-enabled devices are within range, the date, time, distance and duration of interaction can be recorded. The frequency or number of these interactions (analyzed anonymously to form, broadly, measures of clustering or proximal interaction rates over time) may be important given the role of sustained interaction or overcrowding of individuals32,33,34 and contact structure in SARS-CoV-2 transmission35. Furthermore, Bluetooth data in combination with GPS data or a network of Bluetooth sensors can be used to quantify the amount of time people spend at home or other identified locations when lockdown measures are in place to determine if policies are effective.
These data and measures of population-level mobility or clustering patterns would be exceedingly difficult to collect on a similar scale without mobile phone data. These data are often continuously collected, in near real-time, allowing for continued analysis as an outbreak unfolds. Importantly, though, a baseline understanding of contact or clustering patterns prior to any interventions is necessary to inform estimates of intervention impact.
Facilitating contact tracing
Opt-in applications (apps)36,37,38,39,40,41,42 that rely on digital approaches to enumerate and contact individuals who may have been in proximity with someone infected with COVID-19 have been proposed to increase efficiency and decrease the very large burden of manual contact tracing programs43,44,45. By enabling rapid tracing of perhaps higher proportions of affected individuals, these apps can reduce the amount of time that a potentially infected person would have to infect others, particularly in asymptomatic or pre-symptomatic phases of infection46. Most contact tracing apps collect Bluetooth and/or GPS location data to create trails of contacts over a moving time window (14-28 days). Unlike the data needed to understand population-level, aggregated behaviors described above, these data must be linked to single individuals and capture pairwise interactions with other identifiable individuals. Once a case has been identified, they are added to a list of infected users that is queried by the other phones in the network. If the infected user is detected in the trail of contacts, then the user and their contacts are alerted, either by the app or by a public health official, to initiate isolation and quarantine.
This contact tracing process occurs either in a centralized manner, where user information is sent to a remote computer where matching occurs, or in a decentralized manner, where the matching process occurs on the user’s phone. In order for these approaches to feed directly into public health decision making, a direct line between the developers, public health response teams, and users needs to be put in place. This will also be key to mitigating any privacy concerns, which should be dealt with in a transparent and direct manner. Although there has been little discussion to date, routinely collected, individually-identifiable Bluetooth or fine-scale GPS location data may also be used to infer and quantify high-resolution proximity network structures which may further inform contact tracing efforts, but will also raise additional privacy concerns47,48.
Frameworks to process and analyze mobile phone data
Luckily, computing resources and methods to analyze and extract these data will not likely be the limiting factor in these instances. Groups such as Flowminder and Telenor Research Group have worked for multiple years to develop more streamlined processes to analyze these data, particularly aggregate mobility data, that are able to directly interface with mobile phone operators. Flowminder has produced a suite of CDR aggregates, such as counts of active subscribers per region or counts of travelers, that can then be used to calculate indicators of mobility, such as crowdedness, population mixing, locations of interest, and intra-/inter-regional travel49. The code to extract these metrics is publicly available at50. Telenor Research Group works directly with mobile phone operators to provide researchers with spatially aggregated CDR/mobility data51. Facebook’s Data For Good program provides aggregated mobility data to researchers that come from their subscribers, and companies like Cuebiq provided mobility data for a number of COVID-19 studies that summarize the distance users travel or the proportion of users that stay at home52. These existing frameworks – not only the analyses, but also the privacy considerations and data sharing agreements – will provide standardized methods that facilitate integrating mobility data into intervention assessments.
Data privacy
Various forms of identifiable personal information are generated when using mobile phones, including names, identification numbers, fine spatial and temporal data on where the device was used, other users’ identification numbers who may have been detected by Bluetooth, and personal details that might be entered into an app. In light of the growing number of digital privacy concerns and regulations, one must carefully consider the exact form and use of mobile phone data being collected against the legal and ethical need to protect users’ data security and confidentiality. While maintaining user confidentiality is often seen as a hindrance to the use of mobile phone data, in that it limits the use of individual-level data and typically requires aggregation to coarse spatial and temporal resolutions, there are a number of existing frameworks that can help provide guidance for the effective, privacy-conscious use of mobile phone data53.
Exactly which model of data privacy will best suit the use of mobile phone data for COVID-19 response will depend on the exact form and proposed use of the data. As discussed above, there already exist many data processing and analysis frameworks to provide anonymized indicators of population mobility. These standard procedures, though, could result in aggregated data with insufficient spatial and temporal resolution to be effective for monitoring the spread of SARS-CoV-2. Privacy regulations, such as the European Union’s General Data Protection Regulation (GDPR)54, offer exceptions for the use of non-anonymous data that may be needed for other response efforts. For example, opt-in applications for contact tracing may seek consent of the data subject to collect and analyze identifiable data, though the ability to scale opt-in approaches to a wide enough population and to maintain user compliance and participation remains unclear. GDPR and other regulations also provide an exception for anonymization of data to be used in public service, but the regulatory hurdles to gain this exception can be substantial and would require clear use policies and applications for these data. The use of mobile phone data, particularly forms such as those proposed through contact tracing applications, must be weighed against the possible infringements of privacy and civil liberties versus the potential public health benefit. More

Land use mapping and accuracy assessment
According to the land use planning map of Zhuhai city, the characteristics of the city, the status of human activities and land use, and the types of natural ecosystems, we identified and categorized land use into 10 types: woodland, grassland, rainfed cropland, paddy fields, aquaculture areas, reservoirs and pit ponds, tidal flats, rivers and shallow water, built-up land and unutilized land (Supplemental Materials S1: Land use types and descriptions). The ecological land types include woodland, grassland, reservoirs and pit ponds, tidal flats, and rivers and shallow water. Rainfed cropland, paddy fields and aquaculture areas were not included as ecological land types because they are agricultural land mainly used for agricultural production. These land use types are greatly disturbed by humans, their ecological functions are very fragile, and they are affected by economic interests and have low ecological value. Unutilized land provides few ecological benefits and may be converted into built-up land in the short term; thus, its ecological benefits are unsustainable.
After the preprocessing and splicing of multiperiod satellite RS images, we completed object-based multiscale automatic segmentation and land use classification of the images using eCognition Developer software. Specifically, the Estimation of Scale Parameters (ESP) tool was first used to obtain the local variance parameter, which reflects the internal homogeneity of the segmentation object; then, the rate of change (ROC) of the local variance (LV) parameter was calculated37,38. When the ROC reaches its peak, the corresponding segmentation scale can be used as the optimal segmentation scale37. At the optimal segmentation scale, classification is based on the object unit using the nearest neighbor method of eCognition Developer. The nearest neighbor method is a commonly used supervised classification method that is simple and easy to understand, and it is suitable for multiclassification problems39.
Finally, based on the preliminary results data of the four stages automatically classified by eCognition Developer, obvious errors and omissions in the data of the preliminary results were revised and improved through manual visual interpretation. The final revised result data were used for the subsequent analysis of the land pattern and its changes.
This study first drew land use maps for four years: 1991, 2000, 2010, and 2018. We extracted no less than 200 regions of interest (ROIs) in each study year and compared high-resolution Google Earth images to perform a land mapping accuracy assessment. To ensure that the accuracy of each land type was reliably estimated, we confirmed that each land type had at least 10 ROIs when laying out the ROI area. Table 1 shows the land use classification accuracy for the 1991–2018 period. The overall accuracy of the land mapping for 1991, 2000, 2010, and 2018 was 93.4%, 94.1%, 91.1%, and 94.5%, respectively, and the Kappa coefficients were 0.925, 0.933, 0.890, and 0.938, respectively, meeting the research requirements.
Table 1 Classification accuracy of land use types in Zhuhai city.
Full size table

Spatial patterns and dynamics of ecological land
From 1991 to 2018, the ecological land in Zhuhai was dominated by woodland and rivers and shallow water, and the overall area of ecological land continuously decreased (Fig. 1). In 1991, the total area of ecological land was 849.4 km2, accounting for 53.7% of Zhuhai’s urban area. In 2018, the area was reduced to 574.6 km2, accounting for only 36.3% of Zhuhai’s urban area.
Figure 1

The net change in ecological land in Zhuhai city, 1991–2018. The area of woodland is the largest, followed by the area of rivers and shallow water. The proportions of woodland and grassland in the total area of ecological land increased by 7.6% and 1.3%, respectively. Rivers and shallow water and tidal flats showed downward trends, decreasing by 8.7% and 1.8%, respectively. Reservoirs and pit ponds increased slightly and showed dynamic changes.

Full size image

In 28 years, the amount of ecological land decreased by 32.3%, of which woodland decreased by 24.2% (129.6 km2), tidal flats decreased by 67.2% (19.3 km2), and rivers and shallow water decreased by as much as 51.8% (132.3 km2). The reduction in rivers and shallow water represented the bulk of the reduction in ecological land area (48.1%). In contrast, the area of reservoirs and pit ponds grew slightly while maintaining a steady state, increasing by 1.1 km2. Compared with 1991, the grassland area grew slightly, increasing by 5.3 km2, mainly due to the construction of golf courses and parks. Clearly, there is an order of magnitude difference between the increase and decrease in ecological land.
From the temporal perspective (Fig. 2), the change in ecological land mainly occurred in the 1991–2000 period. During this period, the reduction in ecological land was the largest (212.3 km2), mainly distributed in the contiguous area of woodland and built-up land in the central and western areas of the Doumen District and in the coastal areas of the Jinwan District and Xiangzhou District. At the same time, there was a small increase in ecological land, mainly due to the restoration and regulation of tidal flats and reservoirs and pit ponds.
Figure 2

Ecological land gains and losses in Zhuhai city, 1991–2018. (a,c,e) show an increase in ecological land; (b,d,f) show a decrease in ecological land. The decrease in ecological land is obviously higher than the increase, and there is an increase in the degree of patch fragmentation. The reduced patches are mostly marginal woodland and river and shallow water areas. The boundaries of the map come from the Zhuhai Natural Resources Bureau. The drawing of the map was completed with the support of ArcGIS 10.7 software.

Full size image

Since 2000, ecological environmental protection and construction work have gradually been taken more seriously, and the State Council of China promulgated the “National Ecological Environmental Protection Program”. Local governments at all levels have gradually strengthened their awareness of ecological environmental protection. The occupation of ecological land by urban development has rapidly decreased, while the area of new ecological land formed by ecological protection and ecological restoration has gradually and steadily increased. From 2000 to 2010, the ecological land in Zhuhai decreased by 130.1 km2 and increased by 53.6 km2, with a net reduction of 76.5 km2. From 2010 to 2018, the decrease and increase in ecological land were similar, and the net reduction in area was only 18.6 km2; thus, the spatial distribution and quantity of ecological land in Zhuhai city was approximately stable (Fig. 3).
Figure 3

Losses and gains in ecological land area in Zhuhai city, 1991–2018. Green indicates an increase in ecological land, and red indicates a decrease in ecological land. From 1991 to 2000, the net reduction in ecological land was 177.9 km2. From 2000 to 2010, the net reduction in ecological land was 76.5 km2. From 2010 to 2018, the net reduction in ecological land was 18.6 km2.

Full size image

In the 28-year monitoring period of this paper, the reduction in ecological land in the first 10 years (1991–2000) was 0.99 times that in the subsequent 18 years (2000–2018). The total amount of ecological land added in the subsequent 18 years (2000–2018) was 3.5 times that of the first 10 years (1991–2000).
Landscape characteristics
At the landscape level (Table 2), the edge density (ED) of ecological land in the study area is significantly lower than that of nonecological land. The ED exhibited a pattern of first increasing, then decreasing, and subsequently slightly increasing (with values of 33.6 in 1991, 37.7 in 2000, 31.8 in 2010, and 34.7 in 2018). The patch density (PD), landscape shape index (LSI), and largest patch index (LPI) had the same trend as that of the ED. These changes indicate that over time, the landscape of ecological land began to experience an increase in fragmentation and a decrease in regularity and continuity; then, the landscape was reintegrated into a more regular and continuous pattern.
Table 2 Changes in landscape-level indexes in Zhuhai city, 1991–2018.
Full size table

In addition, from 1991 to 2018, the contagion index (CONTAG) of all land in Zhuhai city fluctuated slightly at approximately 55%, and the degree of landscape pattern aggregation did not change much. However, the CONTAG of ecological land was approximately 70%, which was significantly higher than that of nonecological land; this result indicates that the CONTAG and connectivity of ecological land were higher than those of nonecological land. Shannon’s diversity index (SHDI) and Shannon’s evenness index (SHEI) did not change much in the time series, indicating that the landscape diversity of Zhuhai city has basically been stable over the past 28 years. However, compared with 1991, the SHDI and SHEI decreased slightly, indicating that the ecological landscape diversity and uniformity decreased in the study area, while the landscape heterogeneity increased.
At the class level (Table 3), the PD and the area-weighted mean contiguity index (CONTIG_AM) of woodland remained basically unchanged, the LSI increased from 19.99 to 21.7, and the LPI decreased from 9.6 to 3.9. These changes were caused by the following processes: the expansion of built-up land, the preferential occupation of marginal forestland by built-up land, the reduction in the dominance of the landscape type, and the increasing complexity of the original geometry. However, woodland mainly exists in a continuous form, and these encroachment behaviors have little effect on the number, spatial connectivity or proximity of woodland patches.
Table 3 Changes in class-level indexes in Zhuhai city, 1991–2018.
Full size table

The PD and LSI of grassland showed downward trends, while the LPI and CONTIG_AM showed upward trends. This result is closely related to the increase in grassland in the study area. The increased grassland caused the number of patches to increase slightly, improving the superiority of the landscape. The construction of artificial grassland is more regular in the shape of grass patches, and the connectivity is enhanced between landscape units.
In addition, the PD, LSI and LPI of tidal flats showed downward trends, indicating that the development and utilization of tidal flat reclamation were strengthened, the number decreased, and the shape tended to be regular. The landscape characteristics of reservoirs and pit ponds and rivers and shallow water were basically the same: the LSI showed an upward trend, indicating that the patches were seriously disturbed by human activities, the large patches experienced continuous fragmentation, and the landscape type shapes were complicated. In contrast, the LPI showed a downward trend, indicating that activities such as sea filling led to a continuous decrease in sea area.
Ecological quality evaluation
Ecological quality is used to characterize the conditions of the ecosystem; the ecosystem is disturbed by human activities and land use change, and the ability to provide services is also affected40. The value of ecosystem services is an important comprehensive indicator reflecting ecological quality, and the ecological service value of ecological land is higher than that of nonecological land41. Based on the ecosystem service value coefficient proposed by Xie et al.28, we normalized the coefficient value to 0–1 and used the equivalent area and the average equivalent area, which were used to evaluate the ecological service quality of ecological land.
The transformation matrix of ecological land and nonecological land shows the following (Table 4): the probability of ecological land being transformed into nonecological land in the periods 1991–2000, 2000–2010 and 2010–2018 was 25.0%, 19.4% and 14.3%, respectively. The contributions of ecological land to nonecological land were 23.3%, 13.2% and 8.4%, respectively. The transformation of ecological land to nonecological land showed a weakening trend after 2000, and the ecological quality showed improvement.
Table 4 Probability of ecological land being transformed into nonecological land in Zhuhai city, 1991–2018.
Full size table

From 1991 to 2018, the equivalent area of ecological land continued to decrease, but the downward trend gradually stabilized after 2000 (Fig. 4). In 1991, the equivalent area of regional ecological land was 849.4 km2, and in 2000, it was 673.2 km2, indicating a significant decrease in the equivalent area, with a reduction of 20.7%. In 2010, the equivalent area of ecological land further dropped to 600.2 km2, a reduction of 10.8%, although the decrease was significantly smaller than that in the previous period. In 2018, the equivalent area was 574.6 km2, representing a reduction of only 4.3%.
Figure 4

Dynamic changes in ecological land quality in Zhuhai city, 1991–2018. From 1991 to 2018, the equivalent area of ecological land in Zhuhai city showed a downward trend, with a decrease of 274.8 km2, i.e., 32.3%. The average equivalent area index showed an upward trend, with an increase of 0.11, i.e., 9.3%.

Full size image

As shown in Fig. 4, the average equivalent area of ecological land showed a continuous upward trend. Specifically, the average equivalent area was 1.14 in 1991, 1.22 in 2000, 1.24 in 2010, and 1.25 in 2018. This result shows that although the ecological land area decreased, the quality of the ecological land gradually improved. In reality, this pattern was manifested as follows: the area of grasslands and reservoirs and pit ponds gradually increased, the degree of landscape fragmentation weakened, and the landscape dominance became more obvious. In addition, these land types have relatively high ecosystem service values among all land types.
Changes in the center of gravity of ecological land
From 1991 to 2018, the center of gravity of ecological land shifted to the northeast, and the center of gravity of built-up land shifted to the southwest (Fig. 5).
Figure 5

Changes in the center of gravity of ecological land and built-up land in Zhuhai city, 1991–2018. From 1991 to 2018, the center of gravity of ecological land in Zhuhai moved to the northeast by 1346 m. The center of gravity of built-up land moved in the opposite direction, moving 7254 m to the southwest. The boundaries of the map come from the Zhuhai Natural Resources Bureau, and the base map in the main map is the China Online Community Basemap in ArcGIS. The drawing of this map was completed with the support of ArcGIS 10.7 software.

Full size image

From 1991 to 2000, the center of gravity of ecological land moved 404 m to the east and 409 m to the north, and the overall movement was 578 m to the northeast. From 2000 to 2010, the center of gravity of ecological land moved 24 m to the east and 355 m to the north, and the overall movement trend was northward. From 2010 to 2018, the center of gravity of ecological land moved 273 m to the east and 236 m to the north, and the overall movement was 473 m to the northeast. In these three periods, the center of gravity of built-up land moved to the southwest by 2871 m, 3983 m and 424 m. The urban expansion and internal construction mainly experienced a rapid and then slow evolution from the northeast to the southwest.
From the spatial distribution of all ecological land types, the center of gravity of woodland moved to the southeast (0.68 km) from 1991 to 2018. This movement occurred because human construction activities such as deforestation, urban expansion, and infrastructure construction were prominent in the western and northern parts of Zhuhai during the 1991–2000 period. The movement of the center of gravity of grassland to the east and south was highly related to the construction of golf courses, such as the Zhuxiandong Golf Club in the Xiangzhou District, the Dananshan Cuihu Golf Course in Jinding Town, a golf club in the Jinwan District, and Zhuhai Stadium in the Xiangzhou District. The center of gravity of reservoirs and pit ponds moved southward (2.9 km); the center of gravity of tidal flats moved eastward (5.8 km); and the center of gravity of rivers and shallow water moved northward (3.5 km). These changes were closely related to the reclamation engineering carried out by Zhuhai city in recent years.
Modeling the ecological land change process
Changes in urban ecological land are mainly due to the expansion of the outer edge of cities and the oppression of urban internal land development. Therefore, we selected four indicators of natural geography and regional development that might reflect changes in urban expansion and urban construction: elevation, slope, distance from built-up land, and growth rate of built-up land.
With the support of SPSS software, the equation of the transformation probability of ecological land to nonecological land in Zhuhai can be obtained through the binary logistic regression analysis module. Specifically, this equation is expressed as follows (see Supplemental Materials S2: Parameter of the driving factors for modeling):

$$P = 1 – frac{1}{{{1 + }e^{{{ – }left( {{0}{text{.069}} times {text{A } + text{ 0}}{.033} times {text{B } + text{ 0}}{.473} times {text{C } – text{ 1}}{.079} times {text{D } – text{ 0}}{.963} times {text{E } – text{ 0}}{.853}} right)}} }}$$
(1)

where A is the slope; B is the elevation; C is the distance from built-up land; and D and E are the built-up land growth rates of categories 4 and 5, respectively. The squared maximum likelihood of the numerical values (− 2 log-likelihood) of the model was 18,155.4, and the value of the χ2(5) comprehensive test statistic was 7871.2 (p  More

Study design
The study was a single-center, randomized, double-blind, controlled study, stratified by sex in four parallel groups with a 1:1:1:1 allocation ratio: the Test 1, Control 1, Test 3 and Control 3 groups, receiving one (Test 1 and Control 1) or three (Test 3 and Control 3) bottles per day of the Test or the Control product. The study period was split into three subperiods (Fig. 1): a 2-week washout period (day 14 to day 0), a 4-week period of Test or Control product consumption (day 0 to day 28) and a 4-week follow-up period (day 28 to day 56). Dietary restrictions were imposed throughout the entire study period (from day 14 to day 56), with prohibition of the consumption of other fermented dairy products, probiotics, vitamins and mineral supplements, to limit potential interference with the evaluation of the Test product effects. Each subject attended five visits to a clinical unit (Harrison Clinical Research, Munich, Germany): inclusion visit (V1-day 14), randomization visit (V2-day 0), two evaluation visits (V3-day 14, V4-day 28), and an end-of-study evaluation visit (V5-day 56). Blood and stool samples were collected for assessments of eligibility and of the safety evaluation criteria at V1, 2, 3 and 4 (blood) and V2, 3, 4, and 5 (stool). Each visit had to take place within 2 days of the scheduled visit date (± 2 days) to ensure a consistent adequacy between the times of clinical and biological measures and the duration of each corresponding period of product intake or follow-up between subjects. This study was performed in accordance with the principles of the Declaration of Helsinki, the French Huriet law, and ICH-GCP recommendations, and was approved by the ethics committee of the Bavarian Medical Association, Munich, Germany. All volunteers provided written informed consent. This trial was registered on the ClinicalTrials.gov, with the registration number NCT01108419 (date of registration April 22, 2010). The study was funded by Danone Research (France).
Figure 1

Clinical study design.

Full size image

Subject selection
Subjects were screened between March and April 2010, and the study lasted from March 29th 2010 (first subject included) to June 25th 2010 (last subject completed). The following eligibility criteria were assessed at subject inclusion (V1). The inclusion criteria were: male or female volunteers providing written informed consent, aged from 18 to 55 years, with a body mass index (BMI) of 18.5 to 30.0 kg/m2, free-living and considered to be in good health on the basis of a clinical examination, with a normal defecation pattern and either menopausal or with an approved method of contraception if female. Non-inclusion criteria were: any allergy, hypersensitivity to any component of the study product, including lactose, systemic or topical treatment (at the time of inclusion or in the previous 4 weeks) likely to interfere with the evaluation of the study parameters (antibiotics, intestinal or respiratory antiseptics, antirheumatic agents, anti-inflammatory drugs [except for aspirin or equivalent at doses preventing from platelet aggregation or blood clotting] and steroids prescribed for chronic inflammatory diseases), any symptoms of respiratory or gastrointestinal common infectious diseases, a history of chronic metabolic or gastrointestinal disease, abdominal pain or any other severe progressive or chronic disease (cardiac, respiratory, etc.), immunodeficiency, eating disorders or a medicated diet, pregnancy or breast-feeding. The following eligibility criteria were also assessed at the randomization visit (V2): compliance with the dietary and medication restriction (as defined in the non-inclusion criteria) between V1 and V2, negative pregnancy test and parameters within the normal range in the blood samples collected at V1, and absence of common infectious disease symptoms.
Product intervention
The Test product was a fermented dairy drink containing Lactobacillus paracasei CNCM I-1518, Lactobacillus paracasei CNCM I-3689 and Lactobacillus rhamnosus CNCM I-3690 strains, with 107 to 109 colony-forming units (CFU)/g of product, and four yogurt strains (Lactobacillus bulgaricus CNCM I-2787, Streptococcus thermophilus CNCM I-2773, Streptococcus thermophilus CNCM I-2835, Streptococcus thermophilus CNCM I-2778). Counts were measured for each of the bacterial strains present in the Test product, at the start and end of the authorized storage period (shelf life). Means and ranges of strains counts from the batches of product used in the study are provided in Supplementary Table S1. The Control product was a non-fermented dairy drink, acidified with lactic acid and containing pectin as a stabilizer. Both the Test and Control products were sweetened and multi-fruit flavored. Both products were similar in terms of their appearance, packaging, nutritional content (isocaloric) and taste, to ensure the maintenance of double-blinding (both the participants and key study personnel, including the outcome assessors) until the database was locked and the request by the statistician for unblinding (the only staff not blinded being those involved in the preparation of the study products). Products were manufactured in a pilot plant approved by the national health authorities for the production of dairy products for human consumption. They were supplied by Danone Research, France and stored at + 4 ± 2 °C, with a shelf life of 37 days. Analyses were performed to guarantee the absence of microbiological contaminants in all products. Subjects were randomly assigned to the Test or Control group according to a randomization list established before the start of the study by an external statistician. The randomization list contained balanced blocks, stratified by sex, with the allocation of an incremental number linked to product number given by an IWRS system, and was kept confidential at the sponsor’s premises in order to ensure allocation concealment. The subjects were then asked to ingest either one (100 g) or three (3 × 100 g) bottles of the Test or Control product daily, in accordance with their randomization group, for the entire 4-week product-consumption period (28 days). Subjects with three doses per day were recommended to consume no more than two doses at the same time. Compliance was evaluated by the investigator on the basis of the daily reporting of product consumption by each participant in a personal diary and a count of unused bottles.
Outcomes
The primary aim of the study was to compare product safety between the Test 1 and Control 1 groups over the 4-week period of product consumption. The safety evaluation was based on the following parameters: adverse events, physical examination, hematology, metabolism profile, markers of hepatic, kidney and thyroid function, inflammatory markers, bowel habits and frequency of digestive symptoms. Additional information about safety parameters is provided in Supporting Information.
As secondary criteria, safety parameters were also analyzed for the Test 3 and Control 3 groups, over the period of product consumption (V2 to V4), and for both 1 and 3 product doses during other periods: the follow-up period (V4–V5) and the whole experimental period (V2–V5). Stool samples were also subjected to testing to detect and quantify the strains present in the Test product and to analyze the microbiota, for both doses and different study periods (see details and methods below).
Procedure
At each visit, from V1 to V5, subjects underwent a physical examination and vital signs were recorded. Subjects completed a personal diary throughout the 10-week study period, which was collected and examined at each visit by the investigator. This diary included daily reports of study product consumption, the intake of unauthorized products, concomitant medication, symptoms, frequency and consistency of stool and a weekly scoring from the Frequency of Digestive Symptoms questionnaire. The physical activity and smoking habits of the subjects were recorded at each visit. Blood samples were collected for analyses after overnight fasting every two weeks from V1 to V4. The measure of calprotectin concentration, the detection and quantification of strains from the Test product, and the evaluation of the microbiota profile were performed on stool samples collected at each visit from V2 to V5. The study was performed in accordance with the protocol and the statistical analysis plan with no major change during the course of the trial.
Safety monitoring committee
A safety and monitoring committee (SMC), composed of three independent experts in internal medicine, hepato-gastro-enterology and pharmacology, performed an unblinded review of the subject withdrawals, the protocol deviations, the statistical analyses of study parameters and the individual data in the event of abnormal values for safety results. The statistical results were presented after the database lock by the study scientist and statistician to the SMC during two meetings. The SMC then presented its conclusions concerning the safety of the daily ingestion of the Test product at the two doses evaluated.
Stool collection, DNA extraction
We collected fecal samples from 90 subjects at four time points (Test 1 (N = 22), Test 3 (N = 23), Control 1 (N = 21), Control 3 (N = 24)) in RNAlater solution (Ambion, Courtaboeuf, France). Fecal DNA was extracted by mechanical lysis (FastprepFP120; ThermoSavant, Illkirch, France) followed by phenol/chloroform-based extraction, as previously described39. The DNA preparation was subjected to quality control by spectrophotometry on a NanoDrop 2000c spectrophotometer (Thermo Fisher). The DNA was analyzed by quantitative polymerase chain reaction (qPCR), 16S rRNA gene sequencing and whole-genome sequencing.
Quantitative PCR
Three strains, Lactobacillus paracasei subsp. paracasei CNCM I-1518, Lactobacillus paracasei subsp. paracasei CNCM I-3689 and Lactobacillus rhamnosus CNCM I-3690, were quantified by qPCR, as previously described39, with specific primers (Supplementary Table S2). Values were reported as median and interquartile range.
16S RNA gene sequencing, processing and analysis
16S RNA gene sequencing was performed as previously described18. Amplification was performed with the V3-V4 primers for the 16S rRNA (forward: CCTACGGGNGGCWGCAG, reverse: GACTACHVGGGTATCTAATCC). The samples were loaded into flow cells in an Illumina MiSeq 300PE Sequencing Platform, in accordance with the manufacturer’s instructions. Analyses were performed with QIIME (v. 19). The sequences were filtered for quality and a mean of 99,437 ± 36,973 reads per sample were retained. Reads were clustered into operational taxonomic units (OTUs; 97% identity threshold) with VSEARCH, and representative sequences for each OTU were aligned and taxonomically assigned with the SILVA database (v. 119). Alpha-diversity was assessed at genus level. Beta diversity was assessed with Bray–Curtis dissimilarity, Jensen-Shannon divergence, and weighted and unweighted UniFrac on genera and OTUs.
Metagenomic shotgun sequencing and preprocessing
Following standard DNA quality control and quantification, sequencing libraries were prepared with the Nextera XT DNA sample preparation kit in accordance with the manufacturer’s instructions. An overview of the bioinformatic pipeline used in this study is provided in Supplementary Fig. S1. We generated a mean of 35 million (± 8 million) paired-end reads per sample. Read cleaning, filtering and mapping were performed with NGLess version 0.740. An augmented catalog was built from the Integrated Gene Catalog (IGC)41 enriched with genes from the sequencing and de novo assembly of these 107 metagenomes and the seven bacterial genomes present in the Test product (Supplementary Fig. S2). Mapping and count matrix generation were also performed with NGLess. The taxonomic profile was extracted from the count matrix with the Metagenomic Species Pan-Genomes database42. For functional characterization, the catalog was annotated with functional data from the Kyoto encyclopedia of genes and genomes (KEGG, https://www.genome.jp/kegg/)43.
Functional contribution
Metagenomic gene count matrices were aggregated at KEGG orthologous (KO) levels, for the whole gene set and for genes from L. rhamnosus and L. paracasei from the Test product only. We estimated the contribution of the Test product to each KO, by dividing each KO relative abundance level for the Test product by the corresponding value for the whole gene set. A pseudocount of one was added. Corresponding KO relative abundances for the 31 universally distributed marker genes from Ciccarelli et al.44 were also obtained, to estimate the minimal functional contribution of each Test product gene. All KOs for the Test product with a contribution strictly higher than the minimal contribution, constituting a significant functional contribution of the Test product to the gut metagenome, were extracted for downstream analysis. KEGG BRITE and module annotations were used to explore this functional contribution, focusing on enzymes and transporters. We then assessed the extent to which this significant functional contribution set was shared by the other metagenomic species pan-genomes (MSPs).
Statistical analysis
Clinical parameters
No data on adverse events were available to assess the sample size required. The decision to include 24 subjects per group was thus made on the basis of previously published safety studies45,46. For assessment of the safety of consuming the Test product, in comparison to the Control product, adverse events were recorded (MedDRA version 13) and used to evaluate the number of subjects with at least one adverse event, and the total number of adverse events overall, and by relationship to the study product, intensity, seriousness, action taken, and subject outcome. Additional physical examination data, blood parameters, calprotectin concentration in feces, and questionnaires about bowel movements, stool consistency and the frequency of digestive symptoms were collected throughout the period of product consumption and were analyzed as raw data or in terms of clinical significance relative to the baseline value. No formal statistical tests has been performed to assess the safety and study conclusions were based on nominal statistics as described hereafter, on individual data and on overall agreement of the SMC. For quantitative variables, Cohen’s d was calculated for the change from baseline after 4-week product consumption in Test and Control groups as follows: Cohen’s d = (Average raw change from baseline in Test group − Average raw change from baseline in Control group)/Pooled standard deviation at baseline. Cohen’s d values around 0.50 are considered to be of medium magnitude, and those around or above 0.80 are considered to be large47,48. In this study, an absolute Cohen’s d value above 0.5 was considered to be large enough to detect a potential difference between the Test and Control groups. For qualitative binary parameters, the relative risk (RR) and its 95% confidence interval (CI) were calculated by the normal approximation method. Safety analyses were performed on all randomized subjects who had consumed the Test or Control product at least once, i.e. the full analysis set (FAS) population. Statistical analyses were performed with the Statistical Analysis Systems statistical software package version 9.1.3 (Windows XP Professional; SAS Institute, Cary, NC, USA).
Gut microbiota
We used non-parametric tests to analyze qPCR data, alpha and beta-diversity, gene and species richness within individuals, between groups, at baseline and over time. Differential analyses were performed with DESeq2 (version 1.14.1)49 and ZIBR50. For all tests, the alpha risk was set at 0.05 after FDR adjustment by the Benjamini–Hochberg procedure. Network analysis was performed with the SPIEC-EASI R package (version 1.0.751). All statistical analyses were performed, and graphs were plotted with R software (version 3.6.0). Details of the analyses and parameters are provided in Supporting Information. More

Coupling of SM and VPD confounds ecosystem dryness stress
The difficulty to disentangle the respective effects of SM and VPD stems from the fact that SM and VPD are strongly coupled through land–atmosphere interactions7,20. In addition, field experiments that manipulate atmospheric humidity and temperature at the ecosystem scale are lacking21. Given the strong SM-VPD coupling (Fig. 1c), e.g., on the yearly scale, both lower SM and higher VPD are associated with lower ecosystem gross primary production (GPP), indicated by SIF (Fig. 1a, b). This underlies the use of either SM or VPD alone as proxy for dryness stress on ecosystem production in many current models. Note a global spatially contiguous SIF data set was mainly used in this study, which was generated by using the machine-learning algorithm to train SIF observations from Orbiting Carbon Observatory-2 (OCO-2)22. We display the yearly scale because it is typically used to represent the condition of strong SM-VPD coupling globally11, and the study time period mainly spans from 2001 to 2016. However, as SM and VPD are strongly coupled, it is possible that the correlation between SM and SIF is a byproduct of the correlation between VPD and SIF, or vice versa. As a consequence of SM-VPD coupling, the correlations of yearly SM and VPD with SIF is very similar globally (Fig. 1d). Consequently, the correlation between SM and VPD constitutes a confounding factor that is often overlooked when assessing the role of SM and VPD in determining the impact of dryness stress on ecosystem production. There are still low correlations between SIF and SM or VPD in the northern high latitudes or tropical regions, which suggests possible temperature or radiation effects and requires further investigation.
Fig. 1: Strong coupling of soil moisture and vapor pressure deficit confounds ecosystem dryness stress.

a–c Spatial distribution of Pearson’s correlation coefficient between solar-induced chlorophyll fluorescence (SIF) and soil moisture (SM) (r(SIF, SM)), SIF and vapor pressure deficit (VPD) (r(SIF, VPD)), and SM and VPD (r(SM, VPD)), at the yearly scale. Regions with sparse vegetation and regions without valid data are masked in gray. d Relationship between yearly r(SIF, VPD) and yearly r(SIF,SM) across land vegetated areas. Color shows the relative density of data points, with higher density in black and lower density in yellow.

Full size image

Decoupling of SM and VPD globally
At yearly scale, there is a strong negative correlation between SM and VPD, indicating that low SM is always accompanied by high VPD (Fig. 1c), which is consistent with previous findings7,20. From yearly to monthly, weekly, and daily scale, the correlations between SM and VPD are generally decreasing (Fig. 2d), but remain large across extensive areas, such as central South America, Sub-Saharan Africa, India, and Southeast Asia (Fig. 2a and Supplementary Fig. 1). However, when binning the data into 10 bins according to percentiles of either SM or VPD per pixel, we find that the correlation coefficient between SM and VPD in each bin becomes approximately zero (Fig. 2b–d and Supplementary Figs. 2 and 3). This shows that SM and VPD are generally decoupled at daily scale in both SM and VPD bins.
Fig. 2: Decoupling of soil moisture and vapor pressure deficit.

a–c Spatial distribution of Pearson’s correlation coefficient between soil moisture (SM) and vapor pressure deficit (VPD) at daily scale, averaged over daily SM bins, and averaged over daily VPD. Regions with sparse vegetation and regions without valid data are masked in gray. d Violin plots of correlations between SM and VPD from yearly to daily bins across land vegetated areas. White dots indicate the median values, gray boxes cover the interquartile range, and thin gray lines reach the 5th and 95th percentiles.

Full size image

Disentangling the relative role of SM of VPD
We now disentangle the respective effects of SM and VPD in limiting ecosystem production by exploiting the fact that SM and VPD are decoupled in binned daily SM or VPD data (Fig. 2). SM and VPD are also largely decoupled in 4-day bins, which is the temporal resolution of the mainly used SIF data set (Supplementary Figs. 4 and 5). The analysis is guided by the assumption that if SM dominates dryness stress, low SM will limit ecosystem production regardless of VPD variations (Supplementary Fig. 6a, c). In the same way, if VPD dominates dryness stress, high VPD will limit ecosystem production regardless of SM variations (Supplementary Fig. 6b, d).
To illustrate this further, we select an example pixel located in Mali (West Africa). Without decoupling SM and VPD, it is difficult to conclude whether the decrease in SIF is caused by low SM, high VPD, or both in conjunction (Fig. 3a, b). However, when looking at the variation of SIF across VPD gradients in SM bins (without SM-VPD coupling), high VPD does not reduce SIF but even increase SIF a bit under moderate SM conditions (Fig. 3c). In contrast, low SM reduces SIF noticeably in VPD bins (Fig. 3d). This shows that high VPD does not limit SIF in the absence of the SM-VPD coupling at the example pixel, whereas low SM can still limit SIF. In other words, the apparent VPD limitation on SIF is largely the byproduct of SM-VPD coupling. The respective effects of SM and VPD on SIF is also illustrated in Fig. 3e. The changes in SIF from low VPD to high VPD without SM-VPD coupling (termed ΔSIF(VPD|SM)) can quantify the VPD stress on SIF. Likewise, changes in SIF from high SM to low SM without SM-VPD coupling (termed ΔSIF(SM|VPD)) quantify the SM stress on SIF. The effect of SM and VPD on SIF is estimated using two approaches: (i) SIF in the maximum VPD bin minus SIF in the minimum VPD bin or SIF in the minimum SM bin minus SIF in the maximum SM bin; (ii) using linear regression to derive changes in SIF caused by high VPD or low SM. The two approaches lead to similar results (Methods and Supplementary Fig. 16). As shown in Fig. 3f, the SM effect is strong at the example location (ΔSIF(SM|VPD) = −0.17 mW m−2 nm−2 sr−1), in contrast to the VPD effect (ΔSIF(VPD|SM) = −0.03 mW m−2 nm−2 sr−1). Thus, the comparison of (ΔSIF(SM|VPD) and ΔSIF(VPD|SM) enables the disentangling of their relative role in governing dryness stress.
Fig. 3: Disentangling soil moisture and vapor pressure deficit limitation effects.

a Daily solar-induced chlorophyll fluorescence (SIF) versus daily vapor pressure deficit (VPD). b Daily SIF versus daily soil moisture (SM). c Daily SIF versus daily VPD, binned by SM. d Daily SIF versus daily SM, binned by VPD. c, d circles denote the averaged SIF within each bin of VPD and SM. e Average SIF in each percentile bin of SM and VPD. The cyan arrows indicate the VPD limitations on SIF without SM-VPD coupling (ΔSIF(VPD|SM)), and the orange arrows indicate the SM limitations on SIF without SM-VPD coupling (ΔSIF(SM|VPD)). For better readability, only four arrows are shown. f Distribution of ΔSIF(VPD|SM) and ΔSIF(SM|VPD). Circles denote the ΔSIF(VPD|SM) and ΔSIF(SM|VPD) in each bin. Squares denote the corresponding mean. The example pixel is located in Mali, West Africa at 14.25°N, −4.75°E. See Methods for more details.

Full size image

Next, we examine the respective SM and VPD effects on SIF globally. To ensure comparability in space, the SIF time series at each pixel are normalized by the average SIF exceeding the 90th percentile. Temperature and radiation can also limit ecosystem production, therefore, we have filtered out days when other meteorological drivers were likely to be more important than SM or VPD in limiting ecosystem carbon and water fluxes throughout the analyses, following previous studies12,23. We find that ΔSIF(SM|VPD) is negative across most vegetated land areas, robustly indicating the limiting role of low SM to SIF (Fig. 4a, b) and consistent with plant physiological understanding and previous studies4,7. The units refer to the fractions relative to average SIF exceeding the 90th percentile in each grid cell. Large ΔSIF(SM|VPD) are identified in mid-latitudes, including southern North America, central Eurasia, southern Africa, and Australia. In contrast, ΔSIF(VPD|SM) is small and close to 0 across large areas, but it was larger than ΔSIF(SM|VPD) in tropical Africa surrounding the equator (Fig. 4c, d). Globally, a change from the wettest SM to the driest SM under constant VPD reduces SIF by up to 14.9% on average, whereas a change in VPD from lowest to highest quantiles under constant SM has little effect on SIF (−3.8%) on average. Locally, the areas where the strength of SM effects on SIF (|ΔSIF(SM|VPD)|) exceeds that of VPD effects (|ΔSIF(VPD|SM)|) are widespread, which is also visible along the latitudinal gradient (Fig. 4e, f). In total, |ΔSIF(SM|VPD)| is larger than |ΔSIF(VPD|SM)| across 71.3% of land vegetated areas with valid data, by contrast, VPD is more important than SM in 26.7% of corresponding areas. Furthermore, our findings suggest that many previous estimates of the role of VPD on ecosystem production are likely exaggerated16,24 as they did not account for the strong SM-VPD coupling as a confounding factor. In boreal and tropical regions, both SM and VPD have little effect on SIF, which is controlled by radiation and temperature7,25. The spatial patterns of ΔSIF(SM|VPD)—ΔSIF(VPD|SM) are robust to the choice of the particular forcing data set (Supplementary Figs. 7–11). However, when using the GOME-2 SIF and SCIAMACHY SIF with the local overpass time at 9:30 am and 10:00 am, the VPD effects are weaker than that in CSIF (reducing SIF by 0.1% and 0.02% on average globally), including most of Africa (excluding the Sahara) as well as large areas of central South America, southern Asia, and Australia (Supplementary Figs. 9–11). This raise a caveat that using SIF retrieved in the morning would underestimate the VPD effects. To further test the robustness of our result, we standardized the SIF by photosynthetically active radiation (PAR) to remove possible radiation effects26, limited the data to a narrow temperature range to remove possible temperature effects and aggregated data to a coarser time resolution or using 20 percentile bins, yielding similar results (Supplementary Figs. 12–15). Thus, we demonstrate that SM is the dominant factor in driving the response of ecosystem production to dryness at the ecosystem scale across most land vegetated areas, except for tropical and boreal areas.
Fig. 4: Effect of soil moisture and vapor pressure deficit on ecosystem production globally.

a, c, e Spatial distribution of the changes in solar-induced chlorophyll fluorescence (SIF) caused by low soil moisture (SM) (ΔSIF(SM|VPD)) and high vapor pressure deficit (VPD) (ΔSIF(VPD|SM)), and their differences in absolute values (i.e., |ΔSIF(SM|VPD)|−|ΔSIF(VPD|SM)|). b, d, f Zonal means of SM and VPD effects on SIF and their differences in absolute values. The units refer to the fractions relative to average SIF exceeding the 90th percentile in each grid cell. Black lines indicate the mean values, and gray shaded bands show the standard deviation. Regions with sparse vegetation and regions without valid data are masked in white.

Full size image

Different from a recent global assessment of SM stress on ecosystem production that estimates the relation between SM stress and background climate from a small sample of flux sites18, our results build on data with global coverage and hence provide spatially explicit information of SM stress. Further converting the SIF decrease to the actual carbon loss would largely help quantify changes in terrestrial carbon fluxes under drought. Furthermore, our conclusions contradict many laboratory experiments that show strong VPD effects on stomatal conductance at the leaf scale27,28. This again indicates that the stomatal sensitivity to VPD do not definitely determine the same VPD response of plant water and carbon fluxes at the ecosystem scale29,30, but some ecosystem scale measurements reveal that stomatal sensitivity to VPD can matter in some cases11,12. Key processes driving the weak plant photosynthesis response to VPD at the ecosystem scale need to be addressed in future work, such as the role of ecosystem water use efficiency, water storage and hydraulic strategies29.
Dependence of SM stress on climate and vegetation gradients
We find that SM limitation effects (ΔSIF(SM|VPD) are largest in semi-arid ecosystems (Fig. 5a), including shrubland, grassland, and savannah ecosystems. These are the ecosystems that are the main drivers of the interannual variability in global terrestrial CO2 flux31,32. In contrast, VPD effects are much weaker in these regions (Fig. 4c). This suggests that SM could be more important than VPD in driving interannual variability of global terrestrial carbon uptake. As SM stress is strongest in drylands, the projected expansion of drylands33 is likely to increase the influence of SM on the future global carbon cycle. In addition, we find that regions with lower tree fraction exhibit a larger response to SM stress globally (Fig. 5b). This is in line with recent findings34, and further verifies the robustness of our results. Our findings also highlights the differential dryness response of ecosystems along a tree cover gradient.
Fig. 5: Dependence of soil moisture dryness stress on climate and vegetation gradients.

Violin plots of soil moisture (SM) limitation effects (ΔSIF(SM|VPD)) across a aridity gradients and b tree cover gradients. c Violin plots of the sensitivity of solar-induced chlorophyll fluorescence (SIF) to SM (i.e., (frac{{delta SIF}}{{delta SM}}|_{VPD})) within different plant functional types: SHR(S), shrubland (south of 45° N); GRA, grassland; CRO, cropland; WSA(S), woody savanna (south of 45° N); SAV, savanna. White dots indicate the median values, gray boxes cover the interquartile range, and thin gray lines reach the 5th and 95th percentiles.

Full size image

The representation of dryness stress on plant photosynthetic CO2 assimilation can differ largely between TEMs and is considered one of the largest uncertainties in predicting future land carbon uptake and climate35,36,37. Their representations in TEMs often uses an empirical function that only varies by plant functional type (PFT)38, which have generally not been validated against observational empirical data. Therefore, we explored the observed standardized sensitivity of SIF to SM. We find that the sensitivity of ecosystem production to changes in SM can vary largely even in the same PFT with strong observed dryness effects (Fig. 5c). This is consistent with recent findings that the grassland’s sensitivity to dryness can vary greatly39. The differences of dryness response in the same PFT are, e.g., related to plant species, plant height and plant hydraulic processes, such as plasticity variations in xylem and mesophyll conductance, embolism resistance, or water storage40. At present, evaluating and incorporating more plant hydraulic processes into the next generation of terrestrial ecosystems is on the way41. Our results of dryness effects on ecosystem production thus enables an evaluation of further TEM evolution.
In summary, we provide global results of SM and VPD stress on SIF and demonstrate that SM, rather than VPD, is the dominant driver leading to drought limitation on vegetation productivity at the ecosystem level across most vegetated land areas. VPD stress on ecosystem production is almost lost across large areas without SM-VPD coupling. We thus make the case for revisiting the role of VPD in previous studies that neglected the strong SM-VPD coupling. Furthermore, models that do not correctly disentangle the respective VPD and SM limitations cannot adequately predict the dryness stress on ecosystems and associated rough risks to human well-being. The next challenge is to incorporate the observations to constrain the representation of dryness stress on plants in models, which would also reduce uncertainties in the projection of terrestrial CO2 fluxes and associated climate projections. More

1.
Maddrell, S. & Gardiner, B. Excretion of alkaloids by Malpighian tubules of insects. J. Exp. Biol. 64, 267–281 (1976).
CAS  PubMed  Google Scholar 
2.
Després, L., David, J.-P. & Gallet, C. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol. Evol. 22, 298–307 (2007).
PubMed  Article  Google Scholar 

3.
Goulson, D., Nicholls, E., Botías, C. & Rotheray, E. L. Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science 347, 1255957 (2015).
PubMed  Article  CAS  Google Scholar 

4.
Richardson, L. L. et al. Secondary metabolites in floral nectar reduce parasite infections in bumblebees. Proc. R. Soc. B Biol. Sci. 282, 20142471 (2015).
Article  Google Scholar 

5.
Manson, J. S., Otterstatter, M. C. & Thomson, J. D. Consumption of a nectar alkaloid reduces pathogen load in bumble bees. Oecologia 162, 81–89 (2010).
ADS  PubMed  Article  Google Scholar 

6.
Alaux, C. et al. Interactions between Nosema microspores and a neonicotinoid weaken honeybees (Apis mellifera). Environ. Microbiol. 12, 774–782 (2010).
PubMed  PubMed Central  Article  Google Scholar 

7.
Vidau, C. et al. Exposure to sublethal doses of fipronil and thiacloprid highly increases mortality of honeybees previously infected by Nosema ceranae. PLoS ONE 6, e21550 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
McMillan, L. E., Miller, D. W. & Adamo, S. A. Eating when ill is risky: immune defense impairs food detoxification in the caterpillar Manduca sexta. J. Exp. Biol. 221, jeb173336 (2018).
PubMed  Article  Google Scholar 

9.
King, R. L. & Taylor, A. B. Malpighamœba locustae, n. sp. (Amoebidae), a protozoan parasitic in the Malpighian tubes of grasshoppers. Trans. Am. Microsc. Soc. 55, 6–10 (1936).
Article  Google Scholar 

10.
Taylor, A. B. & King, R. L. Further studies on the parasitic amebae found in grasshoppers. Trans. Am. Microsc. Soc. 56, 172–176 (1937).
Article  Google Scholar 

11.
Bailey, L. Honey bee pathology. Annu. Rev. Entomol. 13, 191–212 (1968).
Article  Google Scholar 

12.
Harry, O. G. & Finlayson, L. H. Histopathology of secondary infections of Malpighamoeba locustae (Protozoa, Amoebidae) in the desert locust, Schistocerca gregaria (Orthoptera, Acrididae). J. Invertebr. Pathol. 25, 25–33 (1975).
Article  Google Scholar 

13.
Harry, O. G. & Finlayson, L. H. The life-cycle, ultrastructure and mode of feeding of the locust amoeba Malpighamoeba locustae. Parasitology 72, 127 (1976).
Article  Google Scholar 

14.
Liu, T. P. Scanning electron microscope observations on the pathological changes of Malpighian tubules in the worker honeybee, Apis mellifera, infected by Malpighamoeba mellificae. J. Invertebr. Pathol. 46, 125–132 (1985).
Article  Google Scholar 

15.
Wright, S. H. & Dantzler, W. H. Molecular and cellular physiology of renal organic cation and anion transport. Physiol. Rev. 84, 987–1049 (2004).
CAS  PubMed  Article  Google Scholar 

16.
Gaertner, L. S., Murray, C. L. & Morris, C. E. Transepithelial transport of nicotine and vinblastine in isolated Malpighian tubules of the tobacco hornworm (Manduca sexta) suggests a P-glycoprotein-like mechanism. J. Exp. Biol. 201, 2637–2645 (1998).
CAS  PubMed  Google Scholar 

17.
Rheault, M. R., Plaumann, J. S. & O’Donnell, M. J. Tetraethylammonium and nicotine transport by the Malpighian tubules of insects. J. Insect Physiol. 52, 487–498 (2006).
CAS  PubMed  Article  Google Scholar 

18.
Leader, J. P. & O’Donnell, M. J. Transepithelial transport of fluorescent p-glycoprotein and MRP2 substrates by insect Malpighian tubules: confocal microscopic analysis of secreted fluid droplets. J. Exp. Biol. 208, 4363–4376 (2005).
CAS  PubMed  Article  Google Scholar 

19.
Rossi, M., De Battisti, D. & Niven, J. E. Transepithelial transport of P-glycoprotein substrate by the Malpighian tubules of the desert locust. PLoS ONE 14, e0223569 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Dermauw, W. & Van Leeuwen, T. The ABC gene family in arthropods: comparative genomics and role in insecticide transport and resistance. Insect Biochem. Mol. Biol. 45, 89–110 (2014).
CAS  PubMed  Article  Google Scholar 

21.
Eytan, G. D., Regev, R., Oren, G., Hurwitz, C. D. & Assaraf, Y. G. Efficiency of P-glycoprotein–mediated exclusion of rhodamine dyes from multidrug-resistant cells is determined by their passive transmembrane movement rate. Eur. J. Biochem. 248, 104–112 (1997).
CAS  PubMed  Article  Google Scholar 

22.
Murray, C. L. A P-glycoprotein-like mechanism in the nicotine-resistant insect, Manduca sexta (University of Ottawa, Ottawa, 1996).
Google Scholar 

23.
O’Donnell, M. Insect excretory mechanisms. Adv. Insect Physiol. 35, 1–122 (2008).
Article  Google Scholar 

24.
Berridge, M. J. The physiology of excretion in the cotton stainer, Dysdercus fasciatus, Signoret. IV. Hormonal control of excretion. J. Exp. Biol. 44, 553–566 (1966).
CAS  PubMed  Google Scholar 

25.
Ramsay, J. A. Active transport of water by the Malpighian tubules of the stick insect, Dixippus Morosus (Orthoptera, Phasmidae). J. Exp. Biol. 31, 104–113 (1954).
CAS  Google Scholar 

26.
Maddrell, S. Active transport of water by insect Malpighian tubules. J. Exp. Biol. 207, 894–896 (2004).
PubMed  Article  Google Scholar 

27.
Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
R Core Team. R: a language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2019). https://www.R-project.org/.

29.
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Article  Google Scholar 

30.
Burnham, K. P. & Anderson, D. R. A practical information-theoretic approach. in Model Selection Multimodel Inference 2nd edn (Springer, New York, 2002).

31.
Maddrell, S. H. P. & O’Donnell, M. J. Insect Malpighian tubules: V-ATPase action in ion and fluid transport. J. Exp. Biol. 172, 417–429 (1992).
CAS  PubMed  Google Scholar 

32.
Wieczorek, H., Beyenbach, K. W., Huss, M. & Vitavska, O. Vacuolar-type proton pumps in insect epithelia. J. Exp. Biol. 212, 1611–1619 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Garrett, M. A., Bradley, T. J., Meredith, J. E. & Phillips, J. E. Ultrastructure of the Malpighian tubules of Schistocerca gregaria. J. Morphol. 195, 313–325 (1988).
PubMed  Article  Google Scholar 

34.
Ugwu, M. C., Oli, A., Esimone, C. O. & Agu, R. U. Organic cation rhodamines for screening organic cation transporters in early stages of drug development. J. Pharmacol. Toxicol. Methods 82, 9–19 (2016).
CAS  PubMed  Article  Google Scholar 

35.
Maddrell, S. H. P., Gardiner, B. O. C., Pilcher, D. E. M. & Reynolds, S. E. Active transport by insect Malpighian tubules of acidic dyes and of acylamides. J. Exp. Biol. 61, 357–377 (1974).
CAS  PubMed  Google Scholar 

36.
Hinks, C. F. & Ewen, A. B. Pathological effects of the parasite Malameba locustae in males of the migratory grasshopper Melanoplus sanguinipes and its interaction with the insecticide, cypermethrin. Entomol. Exp. Appl. 42, 39–44 (1986).
CAS  Article  Google Scholar 

37.
Sreeramulu, K., Liu, R. & Sharom, F. J. Interaction of insecticides with mammalian P-glycoprotein and their effect on its transport function. Biochim. Biophys. Acta BBA Biomembr. 1768, 1750–1757 (2007).
CAS  Article  Google Scholar 

38.
Bernays, E. A. & Chapman, R. F. Plant chemistry and acridoid feeding behaviour. Biochem. Asp. Plant Anim. Coevol. 99, 41 (1978).
Google Scholar 

39.
Habig, W. H., Pabst, M. J. & Jakoby, W. B. Glutathione S-transferases the first enzymatic step in mercapturic acid formation. J. Biol. Chem. 249, 7130–7139 (1974).
CAS  PubMed  PubMed Central  Google Scholar 

40.
Stahlschmidt, Z. R., Acker, M., Kovalko, I. & Adamo, S. A. The double-edged sword of immune defence and damage control: do food availability and immune challenge alter the balance?. Funct. Ecol. 29, 1445–1452 (2015).
Article  Google Scholar 

41.
Jeschke, V., Gershenzon, J. & Vassão, D. G. A mode of action of glucosinolate-derived isothiocyanates: detoxification depletes glutathione and cysteine levels with ramifications on protein metabolism in Spodoptera littoralis. Insect Biochem. Mol. Biol. 71, 37–48 (2016).
CAS  PubMed  Article  Google Scholar 

42.
Phillips, J. E. Rectal absorption in the desert locust, Schistocerca gregaria Forskal. I. Water. J. Exp. Biol. 41, 15–38 (1964).
CAS  PubMed  Google Scholar 

43.
Phillips, J. Comparative physiology of insect renal function. Am. J. Physiol. Regul. Integr. Comp. Physiol. 241, R241–R257 (1981).
CAS  Article  Google Scholar 

44.
Proux, J. Lack of responsiveness of Malpighian tubules to the AVP-like insect diuretic hormone on migratory locusts infected with the protozoan Malameba locustae. J. Invertebr. Pathol. 58, 353–361 (1991).
CAS  Article  Google Scholar 

45.
Phillips, J. E. Rectal absorption in the desert locust, Schistocerca gregaria Forskal. II. Sodium, potassium and chloride. J. Exp. Biol. 41, 39–67 (1964).
CAS  PubMed  Google Scholar 

46.
Misof, B. et al. Phylogenomics resolves the timing and pattern of insect evolution. Science 346, 763–767 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

47.
Venter, I. G. Egg development in the brown locust, Locustana pardalina (Walker), with special reference to the effect of infestation by Malameba locustae. South Afr. J. Agric. Sci. 9, 429–434 (1966).
Google Scholar  More

Study areas
Time series of population survey data were used from nonmigratory elk populations in three different locations along the West Coast of the USA (Fig. 4). Five of the populations were in the Prairie Creek drainage (Davison), the Lower Redwood Creek drainage (Levee Soc), the Stone Lagoon area, the Gold Bluffs region, and the Bald Hills region of Redwood National and State Parks (RNSP), Humboldt County, California (41.2132° N, 124.0046° W). These populations occupy an area of about 380 km2. The climate in this region was mild, with cool summers and rainy winters. Annual precipitation was usually between 120 and 180 cm and most of the precipitation fell between October and April. Snow was rare since average winter temperatures rarely dropped below freezing and ranged from 3 to 5 °C. Average summer temperatures ranged from 10 to 27 °C, depending on the distance inland. Elk in RNSP were not legally hunted, and displayed strong social bonding between females, juveniles, and sub-adult males7.
Figure 4

Map of study areas in Arid Lands Ecology (ALE) Reserve, southern part of Redwood National and State Parks, and Tomales Point Elk Reserve in Point Reyes National Seashore. This map was created in ArcMap (Version 10.6; https://desktop.arcgis.com/en/arcmap/).

Full size image

An elk population in the Point Reyes National Seashore inhabited part of the Point Reyes Peninsula in Marin County, California (38.0723° N, 122.8817° W). The elk were restricted to an area of 10.52 km2 on the northern tip of the peninsula by a 3-m-tall fence. The climate of this study area was Mediterranean, with an average annual precipitation of 87 cm27. Most of the precipitation fell from autumn to early spring. Temperatures averaged about 7 °C in winter and 13 °C in summer27,35.
Another elk population was in the Arid Lands Ecology (ALE) Reserve and occupied a 300 km2 area within the U.S. Department of Energy’s Hanford Site, Washington (46.68778° N, 119.6292° W). The climate in this area was semi-arid with dry, hot summers and wet, moderately-cold winters. Average summer temperatures were around 20 °C and average winter temperatures were around 5 °C with an average annual precipitation of 16 cm, half of which fell in the winter as rain36.
Population surveys
In RNSP, females, juveniles, and subadult males were often in the same group and tended to use open meadow habitat more frequently than adult males37,38. These behavioral patterns likely explain why females, juveniles, and subadult males were sighted more frequently than adult males7. Moreover, in size-dimorphic ungulates such as elk, recruitment was strongly correlated with female abundance and weakly correlated with male abundance7,13,39. In RNSP, the abundance of groups of females, juveniles, and subadult males drove the dynamics of the group and of adult males7. Therefore, for the RNSP populations, we used herd counts where a herd was comprised of females, juveniles, and subadult males. We also used herd counts for the Point Reyes and ALE Reserve populations to remain consistent.
Systematic herd surveys of elk were conducted during January from 1997 to 2019 in RNSP. Surveys in the Davison meadows, the Levee Soc area, the Stone Lagoon area, the Gold Bluffs region, and the Bald Hills region were conducted by driving specified routes 4 to 10 times on different days throughout the month of January. The time series for these five herds ranged from 19 to 23 years of data. The elk were counted and classified by age and sex as adult males, subadult males, females, and juveniles. Females could not be visually differentiated into adult and subadult age categories37. The highest count of females, juveniles, and subadult males from the surveys conducted each year was used as an index of abundance of each herd since the detection probabilities were high both on an absolute basis ( > 0.8) and relative to variation in detection probabilities (CVsighting = 0.05)7,40. For the Bald Hills herd, which is the only herd in RNSP where harvests occurred, we added hunter harvests to the highest count of each year to account for this source of mortality. These harvests occurred only when elk from the Bald Hills herd left RNSP.
Elk population surveys were conducted at the Point Reyes National Seashore from 1982 to 2008. Weekly surveys were conducted after the mating season. Surveys were conducted on foot or horseback of female elk that were ear-tagged or had a collar containing radio telemetry32,35. Individuals counted were classified as females, juveniles, subadult males, and adult males. Data were not available for the years 1984 to 1989 and 1993, so the time series included 20 years of data. We used the highest count of females, juveniles, and subadult males in each year in our analyses. This herd was also not hunted.
Elk population surveys in the ALE Reserve were conducted in winters after hunting and before parturition. From 1982 to 2000, biologists used aerial telemetry studies, in which they located all collared elk during each survey and classified them by sex and age. We used the total counts of females, juveniles and subadult males. For years in which multiple surveys were conducted, we used the highest count in each year as an index of abundance for that year25,41. We omitted population survey data collected in 1982 from our analysis because individuals were not classified by sex and age in this year. Consequently, the time series included 18 years of data. For all years of data used, we added hunter harvests to the highest count of each year to account for this source of mortality. The count in 2000 was much lower than in the previous year, likely due in part to a large wildfire which occurred in the summer of 2000, which probably had an immediate effect of reducing available elk forage in the reserve and caused elk to spend more time outside of the ALE Reserve42,43. In addition, the highest recorded number of elk (about 291) were harvested that year43.
Ricker growth models
We fit linearized Ricker growth models simultaneously to the seven time series to estimate population growth parameters as well as temporal variation in r and β. We estimated K as the x-intercept of the Ricker growth model (i.e., when r = 0). Notably, preliminary analyses showed that not accounting for observer error did not bias our results (see Supplementary Information).
We used a Bayesian Markov Chain Monte Carlo (MCMC) algorithm with 3 chains, 150,000 iterations, a burn-in period of 75,000, an adaptation period of 75,000, and no thinning. We used Bayesian inference and MCMC because these methods offer advantages when fitting hierarchical models to model variation in ecological data44,45. We conducted these analyses in the RJAGS program (JAGS Version 4.0.0; https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/) in RStudio (R Version 3.5.0; https://cran.r-project.org/bin/windows/base/old/3.5.0/). We used uninformative priors for the y-intercept (i.e., rmax) and the slope (i.e., β) in order to allow solely the data to influence posterior estimates of these parameters. Informative priors were not necessary as long as parameter estimates from each chain converged. Convergence among chains was determined when the Gelman-Rubin diagnostic ((hat{R})) was less than 1.01, and through visual checks of trace and density plots46.
The estimate of rmax borrowed information among herds because this parameter should be similar among populations within a species22. Therefore, we modeled rmax for each herd (j) as a random effect following a normal distribution with (mu_{{r_{max} }} sim Normalleft( {0, 0.001} right)) and (sigma_{{r_{max} }} sim Uniformleft( {0, 100} right)). To model temporal variation in r for each herd, we included a zero-centered random effect which was also modeled following the normal distribution (gamma_{t,j} sim Normalleft( {0,sigma_{{gamma_{j} }} } right)), where (sigma_{{gamma_{j} }} sim Uniformleft( {0, 100} right)). The estimate of β did not borrow information among herds because this parameter can vary widely among herds18. The prior for β for each herd (j) followed the normal distribution (beta_{j} sim Normalleft( {0, 0.001} right)). To model temporal variation in β for each herd, we modified how we modeled β by using a normal distribution ({beta_{{delta }_{t,j}}} sim Normalleft( {mu_{{beta_{{delta }_{j} }}}} , {sigma_{{beta_{{delta }_{j}} }} } right)), where ({mu_{{beta_{{delta }_{j}} }}} sim Normalleft( {0, 0.001} right)) and ({sigma_{{beta_{{delta }_{j}} }}} sim Uniformleft( {0, 100} right)). Thus, there were four possible Ricker growth models for each herd; (1) no temporal variation in r and β,

$$ r_{t} = r_{max} + beta N_{t} + varepsilon , $$
(3)

(2) temporal variation in r,

$$ r_{t} = r_{max} + beta N_{t,} + gamma_{t} + varepsilon , $$
(4)

(3) temporal variation in β,

$$ r_{t} = r_{max} + {beta_{{delta }_{t}}} N_{t} + varepsilon , $$
(5)

and
(4) temporal variation in both rmax and β

$$ r_{t} = r_{max} + {beta_{{delta }_{t}}} N_{t} + gamma_{t} + varepsilon . $$
(6)

The residual variance was modeled as (varepsilon sim Uniformleft( {0,100} right)). We fit the model with no temporal variation (Eq. (3)) in either parameter to all seven time series simultaneously. All parameters except for rmax were estimated independently for each herd. For each time series of population survey data, we determined whether models with more parameters provided a better fit. We did so by fitting each possible growth model (Eqs. (4)–(6)) to each time series one at a time, while modeling all other time series with no temporal variation in rmax or β (Eq. (3)). The model with the lowest mean deviance from RJAGS by more than 2 was selected for that herd47.
Environmental and demographic stochasticity
We estimated fluctuation in abundance which can be attributed to demographic and environmental stochasticity for herds with different K for each herd. The stochasticity model was outlined by Ferguson and Ponciano9;

$$ Varleft( {N_{t – 1} } right) = Var_{dem} left( {N_{t – 1} } right) + Var_{r} left( {N_{t – 1} } right) + Var_{{upbeta }} left( {N_{t – 1} } right), $$
(7)

where (Varleft( {N_{t – 1} } right)) was total population stochasticity, (Var_{dem} left( {N_{t – 1,} } right)) was population abundance fluctuation due to demographic stochasticity, (Var_{r} left( {N_{t – 1} } right)) was population abundance fluctuation due to changes in r (i.e., density-independent environmental stochasticity), and (Var_{beta } left( {N_{t – 1} } right)) was population abundance fluctuation due to changes in β. The model assumes density-dependent survival following the Ricker model. Demographic stochasticity was calculated as follows;

$$ Var_{dem} left( {N_{t – 1} } right) = alpha N_{t – 1} e^{{ – beta_{Delta } left( {N_{t – 1} } right)}} left( {1 – e^{{ – beta_{Delta } left( {N_{t – 1} } right)}} } right) + sigma_{dem}^{2} N_{t – 1} e^{{ – 2beta_{Delta } left( {N_{t – 1} } right)}} $$
(8)

where (sigma_{dem}^{2}) was assumed to be equal to α9. Environmental stochasticity that is expressed as changes in r, otherwise known as density-independent or additive stochasticity, was calculated as follows;

$$ Var_{r} left( {N_{t – 1} } right) = sigma_{{beta_{Delta } }}^{2} alpha^{2} N_{t – 1}^{2} e^{{ – 2beta_{Delta } left( {N_{t – 1} } right)}} , $$
(9)

and environmental stochasticity that is expressed as changes in β, otherwise known as density-dependent or multiplicative stochasticity, was calculated as follows;

$$ Var_{{upbeta }} left( {N_{t – 1} } right) = sigma_{{beta_{Delta } }}^{2} alpha^{2} N_{t – 1}^{2} left( {N_{t – 1} } right)^{2} e^{{ – 2beta_{Delta } left( {N_{t – 1} } right)}} . $$
(10)

Population growth parameters from the selected Ricker growth model for each herd were used in these equations to estimate each of these sources of stochasticity for each herd across abundances ranging from five to above K. The relative total population stochasticity was expressed as the total population stochasticity at K for each herd divided by that herd’s K. More

1.
Mousseau, T. A. & Fox, C. W. The adaptive significance of maternal effects. Trends Ecol. Evol. 13, 403–407 (1998).
CAS  PubMed  Article  PubMed Central  Google Scholar 
2.
Hagan, H. R. A brief analysis of viviparity in insects. J. N. Y. Entomol. Soc. 56, 63–68 (1948).
Google Scholar 

3.
Resetarits, W. J. Jr. Oviposition site choice and life history evolution. Am. Zool. 36, 205–215 (1996).
Article  Google Scholar 

4.
Schwarzkopf, L. & Andrews, R. M. Are moms manipulative or just selfish? Evaluating the “maternal manipulation hypothesis” and implications for life-history studies of reptiles. Herpetologica 68, 147–159 (2012).
Article  Google Scholar 

5.
Bernardo, J. Maternal effects in animal ecology. Am. Zool. 36, 83–105 (1996).
Article  Google Scholar 

6.
Réale, D. & Roff, D. A. Quantitative genetics of oviposition behaviour and interactions among oviposition traits in the sand cricket. Anim. Behav. 64, 397–406 (2002).
Article  Google Scholar 

7.
McGaugh, S. E., Schwanz, L. E., Bowden, R. M., Gonzalez, J. E. & Janzen, F. J. Inheritance of nesting behaviour across natural environmental variation in a turtle with temperature-dependent sex determination. Proc. R. Soc. B Biol. Sci. 277, 1219–1226 (2010).
Article  Google Scholar 

8.
Seymour, R. S. & Ackerman, R. A. Adaptations to underground nesting in birds and reptiles. Am. Zool. 20, 437–447 (1980).
Article  Google Scholar 

9.
Booth, D. T. Influence of incubation temperature on hatchling phenotype in reptiles. Physiol. Biochem. Zool. 79, 274–281 (2006).
PubMed  Article  PubMed Central  Google Scholar 

10.
Deeming, D. C. in Temperature-Dependent Sex Determination in Vertebrates (eds Valenzuela, N. & Lance, V. A.) 33–41 (Smithsonian Books, 2004).

11.
Deeming, D. C. & Ferguson, M. in Egg Incubation: Its Effects on Embryonic Development in Birds and Reptiles (eds Deeming, D. C. & Ferguson, M. W. J.) 147–171 (Cambridge University Press, Cambridge, 1991).

12.
Schwarzkopf, L. & Brooks, R. J. Nest-site selection and offspring sex ratio in painted turtles, Chrysemys picta. Copeia 1987, 53–61 (1987).
Article  Google Scholar 

13.
Refsnider, J. M. & Janzen, F. J. Putting eggs in one basket: ecological and evolutionary hypotheses for variation in oviposition-site choice. Annu. Rev. Ecol. Evol. Syst. 41, 39–57 (2010).
Article  Google Scholar 

14.
Doody, J. S. et al. Nest site choice compensates for climate effects on sex ratios in a lizard with environmental sex determination. Evol. Ecol. 20, 307–330 (2006).
Article  Google Scholar 

15.
Ewert, M. A., Lang, J. W. & Nelson, C. E. Geographic variation in the pattern of temperature-dependent sex determination in the American snapping turtle (Chelydra serpentina). J. Zool. 265, 81–95 (2005).
Article  Google Scholar 

16.
Doody, J. S. Superficial lizards in cold climates: nest site choice along an elevational gradient. Austral. Ecol. 34, 773–779 (2009).
Article  Google Scholar 

17.
Doody, J. S. & Moore, J. A. Conceptual model for thermal limits on the distribution of reptiles. Herpetol. Conserv. Biol. 5, 283–289 (2010).
Google Scholar 

18.
Delmas, V., Bonnet, X., Girondot, M. & Prévot-Julliard, A.-C. Varying hydric conditions during incubation influence egg water exchange and hatchling phenotype in the red-eared slider turtle. Physiol. Biochem. Zool. 81, 345–355 (2008).
PubMed  Article  PubMed Central  Google Scholar 

19.
Fitch, H. S. Reproductive cycles in lizards and snakes. Univ. Kans. Mus. Nat. Hist. Misc. Publ. 52, 1–247 (1970).
Google Scholar 

20.
Gutzke, W. H., Packard, G. C., Packard, M. & Boardman, T. J. Influence of the hydric and thermal environments on eggs and hatchlings of painted turtles (Chrysemys picta). Herpetologica 43, 393–404 (1987).
Google Scholar 

21.
Muth, A. Physiological ecology of desert iguana (Dipsosaurus dorsalis) eggs: temperature and water relations. Ecology 61, 1335–1343 (1980).
Article  Google Scholar 

22.
Plumer, M. & Snell, H. Nest site selection and water relations of eggs in the snake, Opheodrys aestirus. Copeia 1988, 58–61 (1988).
Article  Google Scholar 

23.
Reedy, A. M., Zaragoza, D. & Warner, D. A. Maternally chosen nest sites positively affect multiple components of offspring fitness in a lizard. Behav. Ecol. 24, 39–46 (2013).
Article  Google Scholar 

24.
Socci, A. M., Schlaepfer, M. A. & Gavin, T. A. The importance of soil moisture and leaf cover in a female lizard’s (Norops polylepis) evaluation of potential oviposition sites. Herpetologica 61, 233–240 (2005).
Article  Google Scholar 

25.
Warner, D. A. & Andrews, R. M. Laboratory and field experiments identify sources of variation in phenotypes and survival of hatchling lizards. Biol. J. Lin. Soc. 76, 105–124 (2002).
Article  Google Scholar 

26.
Li, S. R. et al. Female lizards choose warm, moist nests that improve embryonic survivorship and offspring fitness. Funct. Ecol. 32, 416–423 (2018).
Article  Google Scholar 

27.
Warner, D. A., Jorgensen, C. F. & Janzen, F. J. Maternal and abiotic effects on egg mortality and hatchling size of turtles: temporal variation in selection over seven years. Funct. Ecol. 24, 857–866 (2010).
Article  Google Scholar 

28.
Black, C. P., Birchard, G. F., Schuett, G. W. & Black, V. D. in Respiration and Metabolism of Embryonic Vertebrates (ed Seymour, R. S.) 137–145 (Springer, Berlin, 1984).

29.
Hayes, W. K., Carter, R. L., Cyril, S. & Thornton, B. in Iguanas: Biology and Conservation (eds Alberts, A. C., Carter, R. L., Hayes, W. K., & Martins, E. P.) 232–257 (University of California Press, 2004).

30.
Iverson, J. B., Hines, K. N. & Valiulis, J. M. The nesting ecology of the Allen Cays rock iguana, Cyclura cychlura inornata in the Bahamas. Herpetol. Monogr. 18, 1–36 (2004).
Article  Google Scholar 

31.
Kam, Y.-C. Effects of simulated flooding on metabolism and water balance of turtle eggs and embryos. J. Herpetol. 28, 173–178 (1994).
Article  Google Scholar 

32.
Moll, E. O. & Legler, J. M. The life history of a neotropical slider turtle, Pseudemys scripta (Schoepff), in Panama. Bull. Los Angel.Cty. Mus.Nat. Hist. 11, 1–102 (1971).
Google Scholar 

33.
Tracy, C. R. Water relations of parchment-shelled lizard (Sceloporus undulatus) eggs. Copeia 3, 478–482 (1980).
Article  Google Scholar 

34.
Mortimer, J. A. The influence of beach sand characteristics on the nesting behavior and clutch survival of green turtles (Chelonia mydas). Copeia 1990, 802–817 (1990).
Article  Google Scholar 

35.
Platt, S. G. & Thorbjarnarson, J. B. Nesting ecology of the American crocodile in the coastal zone of Belize. Copeia 2000, 869–873 (2000).
Article  Google Scholar 

36.
Snell, H. L. & Tracy, C. R. Behavioral and morphological adaptations by Galapagos land iguanas (Conolophus subcristatus) to water and energy requirements of eggs and neonates. Am. Zool. 25, 1009–1018 (1985).
Article  Google Scholar 

37.
Thompson, M., Packard, G., Packard, M. & Rose, B. Analysis of the nest environment of tuatara Sphenodon punctatus. J. Zool. 238, 239–251 (1996).
Article  Google Scholar 

38.
Bodensteiner, B. L., Mitchell, T. S., Strickland, J. T. & Janzen, F. J. Hydric conditions during incubation influence phenotypes of neonatal reptiles in the field. Funct. Ecol. 29, 710–717 (2015).
Article  Google Scholar 

39.
Doody, J. S., James, H., Colyvas, K., Mchenry, C. R. & Clulow, S. Deep nesting in a lizard, déjà vu devil’s corkscrews: first helical reptile burrow and deepest vertebrate nest. Biol. J. Lin. Soc. 116, 13–26 (2015).
Article  Google Scholar 

40.
Doody, J. S. et al. Cryptic and complex nesting in the yellow-spotted monitor, Varanus panoptes. J. Herpetol. 48, 363–370 (2014).
Article  Google Scholar 

41.
Doody, J. S. et al. Deep, helical, communal nesting and emergence in the sand monitor: ecology informing paleoecology?. J. Zool. 305, 88–95 (2018).
Article  Google Scholar 

42.
Doody, J. S. et al. Deep communal nesting by yellow-spotted monitors in a desert ecosystem: indirect evidence for a response to extreme dry conditions. Herpetologica 74, 306–310 (2018).
Article  Google Scholar 

43.
Bureau of Meteorology. Average Annual, Seasonal and Monthly Rainfall, https://www.bom.gov.au/jsp/ncc/climate_averages/rainfall/index.jsp (2019).

44.
Cogger, H. Reptiles and Amphibians of Australia (CSIRO Publishing, 2014).

45.
Doody, J. S. et al. Chronic effects of an invasive species on an animal community. Ecology 98, 2093–2101 (2017).
PubMed  Article  PubMed Central  Google Scholar 

46.
Doody, J. S. et al. Invasive toads shift predator–prey densities in animal communities by removing top predators. Ecology 96, 2544–2554 (2015).
PubMed  Article  PubMed Central  Google Scholar 

47.
Shea, G. & Sadlier, R. An ovigerous argus monitor, Varanus panoptes panoptes. Herpetofauna 31, 132–133 (2001).
Google Scholar 

48.
Doody, J. S. et al. Impacts of the invasive cane toad on aquatic reptiles in a highly modified ecosystem: the importance of replicating impact studies. Biol. Invasions 16, 2303–2309 (2014).
Article  Google Scholar 

49.
Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304 (2004).
MathSciNet  Article  Google Scholar 

50.
Telemeco, R. S., Elphick, M. J. & Shine, R. Nesting lizards (Bassiana duperreyi) compensate partly, but not completely, for climate change. Ecology 90, 17–22 (2009).
PubMed  Article  PubMed Central  Google Scholar 

51.
Wilson, D. S. Nest-site selection: microhabitat variation and its effects on the survival of turtle embryos. Ecology 79, 1884–1892 (1998).
Article  Google Scholar 

52.
Refsnider, J., Bodensteiner, B., Reneker, J. & Janzen, F. Nest depth may not compensate for sex ratio skews caused by climate change in turtles. Anim. Conserv. 16, 481–490 (2013).
Article  Google Scholar 

53.
Morjan, C. L. Variation in nesting patterns affecting nest temperatures in two populations of painted turtles (Chrysemys picta) with temperature-dependent sex determination. Behav. Ecol. Sociobiol. 53, 254–261 (2003).
Article  Google Scholar 

54.
Refsnider, J. M. & Janzen, F. J. Behavioural plasticity may compensate for climate change in a long-lived reptile with temperature-dependent sex determination. Biol. Conserv. 152, 90–95 (2012).
Article  Google Scholar 

55.
Georges, A., Limpus, C. & Stoutjesdijk, R. Hatchling sex in the marine turtle Caretta caretta is determined by proportion of development at a temperature, not daily duration of exposure. J. Exp. Zool. 270, 432–444 (1994).
Article  Google Scholar 

56.
Barbault, R. Population dynamics and reproductive patterns of three African skinks. Copeia 1976, 483–490 (1976).
Article  Google Scholar 

57.
Brown, G. & Shine, R. Why do most tropical animals reproduce seasonally? Testing hypotheses on an Australian snake. Ecology 87, 133–143 (2006).
CAS  PubMed  Article  PubMed Central  Google Scholar 

58.
Van Dyke, J. U. in Reproductive Biology and Phylogeny of Lizards and Tuatara (ed Rheubert, J. L.) 121–155 (CRC Press, New York, 2014).

59.
James, C. & Shine, R. The seasonal timing of reproduction. Oecologia 67, 464–474 (1985).
ADS  PubMed  Article  PubMed Central  Google Scholar 

60.
Packard, G. C., Miller, K. & Packard, M. J. A protocol for measuring water potential in subterranean nests of reptiles. Herpetologica 48, 202–209 (1992).
Google Scholar 

61.
Taylor, J. A. & Tulloch, D. Rainfall in the wet-dry tropics: extreme events at Darwin and similarities between years during the period 1870–1983 inclusive. Aust. J. Ecol. 10, 281–295 (1985).
Article  Google Scholar 

62.
de Almeida Prado, C. P., Uetanabaro, M. & Lopes, F. S. Reproductive strategies of Leptodactylus chaquensis and L. podicipinus in the Pantanal Brazil. J. Herpetol. 34, 135–139 (2000).
Article  Google Scholar 

63.
Newton, I. Population limitation in birds: the last 100 years. Brit. Birds 100, 518–539 (2007).
Google Scholar 

64.
James, C. D. & Whitford, W. G. An experimental study of phenotypic plasticity in the clutch size of a lizard. Oikos 70, 49–56 (1994).
Article  Google Scholar 

65.
Jolly, C. J., Shine, R. & Greenlees, M. J. The impacts of a toxic invasive prey species (the cane toad, Rhinella marina) on a vulnerable predator (the lace monitor, Varanus varius). Biol. Invasions 18, 1499–1509 (2016).
Article  Google Scholar 

66.
Christian, K. in Varanoid Lizards of the World (eds Pianka, E. R. & King, D. R.) 423–429 (Indiana University Press, 2004).

67.
Christian, K. A., Corbett, L., Green, B. & Weavers, B. W. Seasonal activity and energetics of two species of varanid lizards in tropical Australia. Oecologia 103, 349–357 (1995).
ADS  PubMed  Article  PubMed Central  Google Scholar 

68.
Warner, D. A., Du, W.-G. & Georges, A. Introduction to the special issue—Developmental plasticity in reptiles: physiological mechanisms and ecological consequences. J. Exp. Zool. A Ecol. Int. Physiol. 329, 153–161 (2018).
Google Scholar 

69.
While, G. M. et al. Patterns of developmental plasticity in response to incubation temperature in reptiles. J. Exp. Zool. Part A Ecol. Integr. Physiol. 329, 162–176 (2018).
Google Scholar 

70.
Siepielski, A. M. et al. Precipitation drives global variation in natural selection. Science 355, 959–962 (2017).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar  More

General
This study was carried out in the Lijiang Forest Ecosystem Research Station, Yunnan Province, China during the period 13 July to 3 August 2019. This field station, which is operated by the Kunming Institute of Botany, Chinese Academy of Sciences, is located to the north of Lijiang on Yulong Snow Mountain at an elevation of 3200 m. Canopy vegetation is dominated by Pinus yunnanensis and Quercus variabilis, and Rhododendron decorum is a highly conspicuous understory shrub species, when in flower. Flowers occur in inflorescences with each plant typically having many inflorescences.
We chose, numbered and bagged one inflorescence on each of 25 plants of Rhododendron decorum, followed the state of individual flowers, and sampled nectar according to a protocol explained below. We selected plants, as encountered, that were flowering and within about 10 m of our walking path, which was along a road and foot track near the field station. We selected inflorescences, one per plant, with at least five unopened buds, and marked five of these buds with small lengths of differently coloured plastic drinking straws22. The different colours enabled us to distinguish flowers during nectar sampling and subsequent measurements of flower colour. All marked flowers were then checked daily to record flower state as bud, beginning to open, open-non-abscised, and open-abscised. Flowers were considered buds if there was no sign of petals unfolding, beginning to open if petals had begun to unfold, and open if petals had unfolded completely. Abscised flowers were clearly indicated by separation between the base of the petals and the rest of a flower, which was along a distinct abscission line (Fig. 1c). Inflorescences were bagged, using green mesh organza bags, to prevent any flower visitation and nectar removal. This species is self-incompatible19, so no pollination occurred.
We carried out two experiments, involving a total of 25 plants. One (Experiment A) involved 10 plants (numbered A-1 to A-10) and was carried out between 13 and 23 July 2019. Experiment B involved 15 plants (numbered B-1 to B-15) and was carried out between 24 July and 3 August 2019. Experiment B was carried out to increase sample sizes for flowers of all ages, and to provide information for relatively young flowers that was not provided by Experiment A (explained further below). Plants were numbered as encountered.
Collection of inflorescences
Inflorescences from Experiment A and Experiment B were collected for sampling of nectar according to the following protocol.
For all inflorescences in Experiment A (i.e., 10 inflorescences) and all in Experiment B, except numbers B-3, 6, 9, 12 & 15 (i.e., 10 inflorescences), each inflorescence was removed from its plant on the first day that abscission of any marked flower was observed. If any marked flower was observed to have abscised, its inflorescence was removed from its plant by breaking its subtending stem and taken to a nearby sheltered ‘nectar sampling station’ where nectar measurements were made. This occurred for flowers between 3 and 9 days of age, counting the first day that a flower was either open or beginning to open as age 1. In a small number of cases, flower abscission occurred when marked flowers were gently touched just prior to nectar sampling. Such flowers were also considered to have abscised.
In addition, five inflorescences from Experiment B (i.e., numbers B-3, 6, 9, 12 & 15) were similarly collected when they were four days old, regardless of whether any flowers had abscised. This provided nectar measurements for relatively young flowers (i.e., ages 1 to 4 days).
Nectar sampling
For collected inflorescences, almost all the marked flowers were open, and we sampled accumulated nectar in each marked and open flower as follows. Nectar was removed using micro-capillary tubes (Hirschmann microcapillary pipettes; 5 µl in Experiment A; 10 µl in Experiment B; both 32 mm long), with volume measured on the basis of nectar length along tube and subsequently converted to µl. When about 0.5 µl of nectar was obtained, this was expelled to a hand-held refractometer (i.e., Bellingham & Stanley, 0 to 50% brix, adjusted for small volumes) for measurement of sugar concentration as % wt/wt sucrose equivalents. These measurements were adjusted for ambient temperature (see Supplementary Information) using a formula developed from information supplied by the manufacturers of the refractometers we use23 and converted to wt/vol using the following formula24: Y = 0.00226 + 0.00937X + 0.0000585X2 where Y is sugar mass per unit volume (mg/µl) and X is % concentration wt/wt. The amount of sugar for a flower (in mg) was then calculated by multiplying nectar volume (µl) by sugar mass per unit volume (mg/ µl).
Nectar was sampled, for both abscised and non-abscised flowers, from where it accumulates after secretion (Fig. 1c). Nectar was sampled for non-abscised flowers from the base of the corolla between the ring of about 10–15 nectaries, around the base of the ovary, and adjacent flower petals. For flowers that had abscised, nectar was separately sampled from both the ring of nectaries and the inside lowest 5 mm of the flower petals, where some nectar becomes attached.
Some flowers were judged to have been affected by rain and their nectar concentration measurements were excluded from analyses. There were periods of rain during our study and occasionally nectar concentration readings of lower than 1.5% wt/wt were obtained (n = 6), and the nectar assumed to have been diluted by rainwater. These records were excluded from analyses. Fortunately, most flowers pointed downwards and were thus not affected by rain.
Flower colour and pigment
Flower colours were measured by means of a modified Panasonic GH-1 camera. The low-pass filter of the camera had been removed in order to increase the sensitivity for ultraviolet light. The camera body was combined to an Ultra-Achromatic-Takumar 1:4.5/85 lens made of fused quartz that transmits UV light. Since the modified camera is sensitive to ultraviolet and infrared light, a UV-/IR-Cut filter transmitting light between 400 nm and 700 only nm was used to capture a normal reference picture. In addition, a UV-picture was captured from the identical position using a Baader UV-filter that transmits near ultraviolet light only. A white Teflon disc reflecting equal amounts of light in a range of wavelength from 300 to 700 nm was used for manual white balance before taking pictures. Using Image J both pictures were split into the RGB color channels, and then a false color photo was merged using the green channel of the color picture as red, the blue channel of the color picture as green, and the blue channel of the UV picture as blue (see Supplementary Information). For more details see article by Verhoeven et al.25. Using IrfanView image’s histogram a uniform non-decomposed area (number of pixels  > 10,000) of the adaxial corolla on the false color picture was selected. The average intensity for the red, green and blue channel of the false color photos with values between 0 and 255 was used for color evaluation. Abscised and non-abscised flowers were photographed together enabling direct comparison of the colours of the flowers.
Pigment content was deduced from the sum of the values of the red, green and blue channel of the false color photos. Since abscised and non-abscised flowers both appear white to the human eye, the possible change in the content of a UV-absorbing pigment was checked by comparing the value for the blue channel in relation to the sum for the values of the green and red channels.
Recordings of the spectral reflectance were done with an abscised and an open, non-abscised flower from each of five inflorescences. Reflectance measurements were performed with an USB2000 + spectrophotometer (Ocean Optics) and illumination was provided by a DH-2000-BAL light-source (Ocean Optics), both connected via a coaxial fibre cable. All measurements were taken in an angle of 90° to the measuring spot with a pellet of barium sulphate used as white standard and a black piece of plastic used as black standard.
Analyses
We used the General Linear Model approach to determine relationships for all flowers between nectar attributes (i.e., volume—µl, concentration—wt/vol, sugar weight—µg) as dependent variables and flower age, whether abscised, experiment (i.e., A vs. B), and Plant ID as independent variables. We also treated Plant ID as an independent categorical variable, but nested within experiment.
We used ANOVA to evaluate relationships between reflectance intensity and whether flower abscised, across different false colours, with Kolmogorov–Smirnov test for normality and Tukey post-hoc comparisons between means. We took log intensity as the dependent variable in order to meet the normality assumption.
We compared spectral reflectance for abscised and open, non-abscised flowers on the basis of the average reflectance across all wavelengths. Here we assumed that the 1140 reflectance values for each flower could be combined into a single average measure and that this average measure adequately represented each flower. We compared the two groups of flowers with a Kolmogorov–Smirnov Two Sample Test.
All analyses were carried out using the software SYSTAT26. More