Ecology

Heterogeneous and affluence-influenced footprints of APAC
The amount of natural resource extractions, environmental emissions, and socio-economic influences associated with satisfying the final demand of average person vary significantly among the eight APAC countries/regions (Fig. 1). The variations correlate with differences in affluence levels (Supplementary Table 2), which is consistent with the findings of previous footprint research focusing on GHG emissions4, and blue water consumption42 across countries worldwide. While the footprints of lower-income countries (e.g., India, Indonesia, RoAP, and China) have also increased with poverty alleviation, yet, most of them are still below the global averages by 2015. In addition, based on the varying environmental footprints (i.e., blue water, energy, and GHG) of the region’s six countries assessed annually from 1995–2015, we find some evidence for the Environmental Kuznets Curve (EKC) hypothesis, i.e., environmental pollution first rises and then falls as economic development proceeds (Supplementary Fig. 1 and Supplementary Fig. 2). Previously, few studies have addressed the socio-economic implications driven by final demand, such as the labor inputs required (employment) and revenues generated (value added)1. Here we find that a country/region’s employment and value-added footprints are both positively correlated with its affluence level.
Fig. 1: The environmental-social-economic footprints of APAC countries/regions in 1995 and 2015.

The six footprint indicators fall into three categories and are presented in a natural resource, b local and global environmental threats, c socio-economic effects. The eight countries/regions are aligned on the y-axis following a descending order by average income per capita in 2015. Footprints are further broken down to the impacts occurred within the country/region, abroad in other APAC countries/regions, and abroad in non-APAC countries/regions. In each subplot, the dashed line indicates the world average level. GHG (CO2, N2O, and CH4) are measured in CO2 equivalent based on the 100-year Global Warming Potential (GWP100). Note, in the plots for value added, we use the insets to show the footprint compositions of the last four countries/region. China refers to Chinese mainland, and we call Taiwan, China as Taiwan for short in the rest.

Full size image

Affluence levels also affected the geo-compositions of the footprint indicators in the APAC region. Richer APAC economies showed high and increasing reliance on outsourcing blue water consumption, PM2.5 emissions, and labor abroad, especially within the APAC region. Specifically, 58, 56, and 52% of the blue water footprints of Japan, South Korea and Taiwan were traced to water consumed in other APAC economies in 2015 (Fig. 1), primarily through the intra-APAC trade of agricultural products, such as sugar cane/sugar beet and paddy rice (Supplementary Fig. 3). In contrast, the final demand of middle- to low-income APAC economies has been largely satisfied by domestic natural and labor resources over the two decades. Moreover, during the same period, the PM2.5 footprint of high-income APAC economies (Australia, Japan, South Korea, and Taiwan) that occurred in other APAC economies increased from ~17 to 40%. Yet, for the less wealthy countries (e.g., China, India, and Indonesia), 75−99% of PM2.5 footprints were indigenous, a considerable fraction of which are attributed to direct household emissions (Supplementary Fig. 3). The consumption of traditional fossil fuels and biomass (e.g., fuelwood and agriculture wastes) for home cooking and heating, and the open-air combustion of biomass, especially in rural areas of India and China has been considered a significant source of PM2.5 emissions43,44,45.
Among all the indicators, we find the value-added footprints of the APAC economies, i.e., the contributions of their final demand to global economic growth (domestic + intra-APAC + non-APAC) experienced the most significant increases. Such increases are especially significant for the middle- and low-income APAC economies. For example, China’s per capita value-added footprint increased by eight times during 1995−2015, followed by India (217%), Indonesia (113%), and RoAP (82%). The rates are much higher than those experienced by the high-income economies (28% on average).
The APAC economies have grown more interdependently linked during the past two decades. Depending on the footprint indicators, the foreign footprint shares traced to APAC countries (Abroad, APAC in Fig. 1) increased by 4–27% while the domestic shares decreased by 3–33%. Given that the compositional changes within the aggregated RoAP cannot be estimated, the variation of RoAP is not included in the ranges. The increased intra-APAC interdependencies are primarily attributable to the strengthened linkages among the six major APAC countries. In contrast, India’s environmental-social-economic footprints remained predominantly domestic (90–99% in 1995 and 82–97% in 2015). Previously, researchers highlighted that resource constraints have already become a bottleneck for India’s social and economic development, such as fresh water scarcity46 and energy deficiency47. Air quality deterioration caused by PM2.5 emissions has also made India under severe health burden48. Our results here confirm that India remains highly dependent on local resource use and labor-intensive production activities to sustain its socio-economic growth. Engaging in international trades has the potential of reducing the depletion of local scare resources (e.g., blue water and energy) through import while adding employment and added-value through exporting goods produced with abundant labor resources and low environmental impacts locally1,30. The trade-environment relationship is primarily rooted in the economic principle of competitive advantages among countries for international trade. Therefore, adopting the strategy of reducing trade barriers rationally, increasing the openness to the outside world to actively promote international economic cooperation may be one solution to alleviating the pressure of domestic resource depletion (i.e., water and energy) and environmental damages (i.e., carbon and air pollutions) during India’s economic growth.
Footprint outsourcing and disparity through intratrade
The linkages and imbalances among the APAC countries/regions, through the virtual flows of environmental-social-economic resources embedded in intra-APAC trade, are highlighted in Fig. 2 (2015) and Supplementary Fig. 4 (1995). The intra-APAC trade patterns observed over the past two decades confirm that developed economies (Japan, Australia, South Korea, and Taiwan) import natural resources and labor from less-developed regions where resources and labors are cheaper (China, India, and RoAP). In 2015, for the 37 billion m3 of net bilateral virtual water trade in the region, nearly 80% was associated with the exports from RoAP and India, while ~50% was driven by the final demand of Japan and South Korea, mainly embedded in a range of water-intensive agricultural crops and products. Yet, as mentioned above, India has already been threatened by serious water crises with low availabilities of safe drinking water49. For the 15 EJ net bilateral energy flow embedded in the 2015 intra-APAC trade, China was the main net supplier to the rest of the region, contributing 43%, followed by South Korea (27%), which possesses strong petrochemical and steel industries. China’s net energy outflows were predominately embedded in energy-intensive products (e.g., coal power, steel, and gasoline), 73% of which ultimately served RoAP’s final demand. Of the 488 Tg net bilateral flows of GHG emissions, 76% was supplied by China, while the final demand of RoAP contributed more than half of it (286 Tg), and the rest is attributed to the final demand of four higher-income economies (Australia, Japan, South Korea, and Taiwan, 123 Tg in total). 81% of the 677 thousand tons of net bilateral PM2.5 flows occurred in China and India. And 38% of the net flows were driven by the final demand of four higher-income economies. The social and economic effects of the intra-APAC trade are more nuanced than those related to natural resources and emissions. The Intra-APAC trade resulted in a net bilateral outpouring of 99 million people, i.e., labor resources, in 2015. RoAP was the largest net labor supplier, contributing 77%, while Australia, Japan, South Korea, and Taiwan had net inflows of labor to satisfy their final demands, equivalent to employing 11.4, 27.6, 9.3, and 5.6 million people from the rest of APAC, respectively. In terms of net flows of value added, the intraregion pattern appears significantly different and almost opposite from those of other indicators. By net outpouring value added to other APAC economies, China turned from the second trade surplus country in 1995 to the largest one within APAC in 2015, followed by the developed economies. Yet, on the per capita basis, South Korea, Australia, Japan and Taiwan achieved the most prominent economic gains through the intra-APAC trade.
Fig. 2: Net environmental-social-economic virtual flows of the intra-APAC trade in 2015.

The top five net flows are shown for each footprint indicator, and fall into three panels a natural resource, b local and global environmental threats, c socio-economic effects. The width of the arrows in each panel represents the magnitude of the net flow within the APAC region. The background colors indicate the specific net footprint (import-export) per capita of each region/country. The negative net footprint indicates net displacements (of resource use, emissions, labors, and economic value added) to other APAC countries/regions. Note, the background color and flows of Taiwan, China could not be shown on the map.

Full size image

Our results further demonstrate that the economic and environmental inequity owing to the intraregional trade worsened along time. From 1995 to 2015, the majority of the environmental externalities (i.e., more than 90% of the virtual water flows) were shifted from higher-income to lower-income countries, with a considerable worsening trend for PM2.5 (Supplementary Table 3). As for the economic gains associated with the intraregional trade, higher-income countries’ share increased from 38 to 59% from 1995 to 2015. At country level, China experienced the most significant transformation in intra-APAC trade, especially in the virtual trade of water and labor. For virtual blue water flows, China turned from the second largest net exporter of the region in 1995 to the third largest net importer in APAC (5.3 billion m3), largely due to the reverse of the import-export relationship with RoAP. Over the same period, the net virtual water export of India surged by 161%, despite the aggravating water stress concerns. China turned from the largest net supplier of employment (Supplementary Fig. 4) to the second net demander in APAC after two decades, while RoAP showed an opposite trend. For example, the labor- and resource-intensive textile and apparel production industries have been shifted from China to other less-developed countries, such as Vietnam and Bangladesh50. Such a role switching stemmed from the rising labor costs in China, which is expected to drive more low-end manufacturers to low-cost foreign economies in the coming years51,52. Overall, based on the consistently calculated virtual flows of multiple indicators, we highlight a growing environmental-socio-economic disparity within the APAC region, owing to the intraregional trade. The undesirable environmental externalities are primarily and increasingly shifted from higher to lower-income economies, while higher-income economies achieved more economic gains.
The role of APAC and intra-APAC trade in globalization
APAC is becoming an important player in the globalization process: as natural resource suppliers and manufacturers for the rest of the world, for managing global environmental emissions, and in the labor and monetary markets (Fig. 3). By 2015, APAC-related share of these categories had surpassed 50% (the red, yellow, and blue parts in Fig. 3). In contrast, the international trade without APAC countries (gray in Fig. 3) shrank for all the indicators. Moreover, for all the footprint indicators, the intra-APAC trade and non-APAC’s outsourcing to APAC (red and yellow parts in Fig. 3, respectively) account for a considerable and increasing fraction of the global trade. The former grew from 17 to 20% on average and the latter from 23 to 27% on average. APAC’s outsourcing to non-APAC countries (blue in Fig. 3) only grew slowly, with an average of 14 and 16% in 1995 and 2015, respectively. Earlier studies have proven that the world’s dominating embedded labor flows originate from developing countries, predominately to satisfy the final demand of wealthier economies23. Here, we further elucidate that non-APAC economies became more dependent on offshoring the resources, emissions and labor-intensive industries to APAC over the investigation period.
Fig. 3: Natural resources, environmental emissions, and socio-economic factors embedded in exports.

a natural resource, b local and global environmental threats, c socio-economic effects are shown in three panels. Note, the red and blue sections correspond to the same components in Fig. 1.

Full size image

Among the six indicators, we find APAC’s share in global value added was the smallest. In particular, APAC’s economic gains from exports (red + yellow) were significantly smaller (27–32%) than APAC’s resources, emissions, and labor embedded in its exports (40–48%). Intuitively, this may be due to the fact that the APAC region is dominated by population from developing countries, thus resource/labor-intensive and low value-added products dominate the region’s export. In comparison, non-APAC countries are relatively more skilled in producing products and with higher value added and lower intensities of resource and emissions4. Also, APAC’s relatively small share in the value-added dimensions of the global supply chains can be a result of the overall low resource efficiency of the region12. As resource extractions and environment emissions become increasingly outsourced to the region, where environmental regulations and efficiency measures are only emerging, more environmental impacts will likely be resulted on the global level, offsetting or even reversing the resource efficiency gains and climate change mitigation efforts achieved in developed countries.
On the positive note, from 1995 to 2015, APAC’s resource and environmental intensities declined substantially, both from the perspective of footprint, i.e., footprint per final expenditure, and from the perspective of trade, i.e., direct impacts/gross trade (see Materials and Methods). Such improvement is most noticeable for labor (Table 1). More specifically, over the past two decades, APAC improved faster than the world averages in blue water, PM2.5, and labor requirements per expenditure of final consumer, which decreased by 30, 39, and 35%, respectively, while lagging behind the global averages in reduction pace of energy and GHG intensity from footprint perspectives. The intra-APAC trade achieved a reduction in the intensities of energy, GHG and labor by 22, 42, and 34%, respectively. Relatively, this was higher than the world averages over the period. However, the PM2.5 embedded in intra-APAC trade per traded values even grew by 7%, whereas the same indicator decreased by 34% at the global scale. This can be explained primarily by the PM2.5 trade magnitude in APAC nearly tripling in 2015, much higher than that of the global average (0.8 time). This indicates that the PM2.5 issue has become graver for the APAC trade after two decades. The ratio of traded value added to total intratrade values measures the economic gains of trade. Value added is a general indicator of economic performance, yet, it could not evaluate all the advantages of trade, which may cause misleading results in some cases. Therefore, more comprehensive indicators need to be incorporated to assess its advantages.
Table 1 Intensity comparison between APAC and the global average from the footprint and trade perspectives.
Full size table

Despite APAC’s improvement, APAC’s resource and environmental intensities remain higher in general than the world averages in 2015. The comparison of intensity changes at the national and regional levels was also shown in Supplementary Fig. 5. This can be attributed to the characteristics of the APAC region and may not be altered soon. Specifically, the socio-economic development is supported by resource- and emission-intensive productions in primary industries, especially those linked with capital development and those to satisfy the export demand. The APAC region has a long way to go to achieve a greener and more sustainable trade pathway.
Our result further reveals that, China played a crucial role in APAC’s transformation to greener trade and greener consumption. Generally, without China’s improvement, all the footprint intensities of APAC would exhibit a further increase from 1995–2015, at 8% on average. By contrast, the value is −20% with China in the picture. For the trade intensities, the decline ranges would be much smaller (−1% on average) without China as opposed to with China (−24% on average). Another noteworthy finding is that when eliminating China from APAC, the intensities of footprint and trade in 1995 would be 20–30% lower, implying that China had turned from the lagger to the leader of efficiency improvement in the APAC region. More

Study area
This study was conducted at the Sierra Morena mountain range in Southwestern Spain (38° 22′ 50.64″ N, 4° 45′ 27.69″ W; Pozoblanco, Córdoba). This area has a continental-Mediterranean climate characterised by cold winters and hot, dry summers. Annual precipitation is 439 mm year−1 and mean annual temperature is 15.2 °C. The studied ecosystem is a “Mediterranean dehesa”, a savannah-type forest managed as a traditional sylvo-pastoral system, with an herbaceous stratum of native pasture and a tree layer of scattered oak trees (Quercus ilex L.). Tree density is 14.5 ± 1.3 trees ha−1, and the herbaceous layer is dominated by a diverse community of annual native species such as: Hordeum murinum (L.), Senecio vulgaris (L.), Bromus madritensis (L.), and Sinapis alba (L.), with a mean species richness of 9.6 ± 0.3 species m−2.
Experimental design
In September 2016, we selected three sites separated by 4.3 ± 1.9 km with similar tree density, topography, slope, orientation, and soil type. At each site, we sampled two areas: a partially shaded site under Q. ilex canopy (under tree canopy) and a nearby open grassland site (open grassland). These two habitat types are characteristic of the dehesa ecosystem in the region. Then, within each area (i.e. habitat type) at each site we established six permanent 4 × 6 m (1.2 m high) fenced plots (N = 36 plots). Each plot was divided into two subplots (N = 72 subplots), one subplot was subjected to a rainfall reduction treatment whereas the other was exposed to ambient rainfall. Rainfall was reduced by placing a 2.5 × 2.5 m rain-exclusion shelter over one of the subplots (Fig. 2), resulting in a 30% reduction in annual rainfall22 as predicted by the IPCC for Mediterranean regions23. In addition, within each subplot we sampled two ca. 1-m2 quadrats (144 quadrats in total, i.e. 2 quadrats × 72 subplots): one was subjected to a warming treatment and the other was not manipulated (i.e. ambient temperature). Warming was achieved by placing a 40 × 50 × 32 cm (ca. 0.65 m2) hexagonal open top chamber of methacrylate material without UV filter (Faberplast, Madrid) at the centre of the quadrat (Fig. 2), increasing air temperature by ca. 2 °C compared to ambient temperature (17.5 ± 0.1 vs. 15.7 ± 0.1 °C, respectively). This matches the temperature increase for the study region by climatic forecasting models (SRES A-2 model by the Intergovernmental Panel on Climate Change23). Ambient air temperature was slightly lower (approx. 0.5 °C) in the shaded habitat than in the open habitat, but the experimental temperature increase was of similar magnitude in both cases. Temperature in each quadrat was measured with data loggers. Based on the above set-up, the experiment followed a randomized split–split plot design replicated across three sites, with habitat type as the whole factor (plots as replicates), rainfall manipulation as the split factor (subplots), and temperature manipulation as the split–split factor (quadrats).
Figure 2

Experimental design. Picture showing open top chambers and rain-exclusion shelters used for our climate manipulations. Photo credit: Ignacio M. Pérez-Ramos.

Full size image

Sampling, measurements and chemical analyses
In April 2017 and April 2018, i.e. 1 and 2 years after establishment of climatic treatments, we identified all plant species and estimated their frequency within each quadrat. Then, in April 2018, after the second survey of species frequencies, we collected 3–4 fully expanded leaves from each plant for chemical analyses which were pooled to obtain a single sample per species and quadrat. Plant sampling was restricted to a central portion of 0.65 m2 in each quadrat (within chamber in the case of warming quadrats) as a function of the surface covered by and sampled within the hexagonal chamber placed in quadrats subjected to warming (see above). Due to lack of plant material, we only sampled 127 out of the 144 quadrats. In total, we collected 436 samples from 28 annual plant species (Table S1). Immediately after collection, we oven-dried leaves for 48 h at 40 °C, ground them with liquid nitrogen, and stored them at room temperature before conducting chemical analyses.
We determined total phenolic content colorimetrically by the Folin-Ciocalteu method24,25. Briefly, we extracted phenolics from 20 mg of plant tissue with 70% methanol in an ultrasonic bath for 15 min, followed by centrifugation. We determined total phenolic content colorimetrically by the Folin–Ciocalteu method in a Biorad 650 microplate reader (Bio-Rad Laboratories, PA, USA) at 740 nm, using tannic acid as standard. We performed three technical replicates of each sample in order to estimate variations due to the experimental procedure, and expressed concentrations based on dry weights (d.w.). We decided to use total phenolics as these could be measured across all species in order to address the goal of testing for community-level effects of climatic stressors on plant secondary chemistry.
Statistical analyses
For each quadrat, we calculated the CWMphenolics using the weighted.mean function in R software26. For this, we used phenolic values at the species level weighed by the abundance (i.e. frequency) of each species. We then performed linear mixed models using PROC MIXED in SAS 9.4 (SAS Institute, Cary, NC)27 to test for the effects of habitat type, rainfall manipulation, temperature manipulation, and their two-way and three-way interactions (all fixed factors) on CWMphenolics. We included site, the site × habitat type interaction, and the site × habitat type × rainfall manipulation interaction as random factors. The two- and three-way interaction interactions tested the main effect of habitat type and rainfall manipulation (respectively) with the appropriate error terms according to a split-split plot design. We log-transformed original values to achieve normality of residuals.
We also assessed whether climatic manipulations influenced the community-level expression of phenolic compounds through changes in plant composition with a PERMANOVA (‘vegan’ package in R) that included the effects of climatic treatment, year, and their interaction. More

Origin of hydrate-bound hydrocarbons
A relationship between C1/(C2 + C3) and C1 δ13C has been applied to identify the sources of hydrocarbons in submarine seeps24. Recently, this diagram was revised based on a large dataset25. As shown in Fig. 4a, hydrate-bound hydrocarbons at the Kedr MV have thermogenic and/or secondary microbial origins, whereas those of other gas hydrate sites (Malenky, Bolshoy, Malyutka, Peschanka P-2, Kukuy K-0, Kukuy K-2 and Goloustnoe; Fig. 1) in Lake Baikal demonstrate microbial or early mature thermogenic origins. The hydrate-bound C1 from all locations except those at the Kedr MV were interpreted to be of microbial origin via methyl-type fermentation23 according to Whiticar’s old diagram26; however, the revised diagram25 suggests early mature thermogenic gases (Fig. 4b). Those of the Kedr MV plot at the boundary of the thermogenic and secondary microbial origin zones. Low C1 and C2 δ13C at the Peschanka P-2 MV indicated that C1 and C2 are of microbial origin27,28, whereas Kedr MV shows high C1 and C2 δ13C indicating their thermogenic origin (Fig. 4c). At other sites, C1 and C2 δ13C suggested that gases are mainly of microbial origin (in terms of C1) with some thermogenic component (13C rich and higher concentration in C2).
Stable isotopes in hydrate-bound C1 at the Kedr-1 and Kedr-2 areas suggested its thermogenic origin. However, it is close to the field of secondary microbial C1 in Fig. 4b, and the data are plotted in the overlap between the fields of thermogenic and secondary microbial in Fig. 4a. Milkov29 mentioned that secondary microbial C1 is characterised by C1-rich dry gas, large C1 δ13C (between − 55‰ and − 35‰) and large CO2 δ13C (more than + 2 ‰). Although hydrate-bound and sediment gases in the Kedr MV were not C1 rich and contained 3%–15% of C2, C1 δ13C was around − 45‰, which agrees with the secondary microbial C1. Because some data of secondary microbial gas are plotted outside the field on the original graph25, we could include the gas data in the category of secondary microbial C1 in Fig. 4b.
Figure 6 shows the relationship between C1 δ13C and CO2 δ13C in the sediment gas obtained using headspace gas method. According to the genetic diagram25, gas hydrate cores are plotted at the zones of the thermogenic and secondary microbial origins, whereas the cores at the peripheral area are primary microbial. The headspace gas data of the hydrate-bearing cores in Fig. 6 seem to be plotted in the field of thermogenic gas (low CO2 δ13C), but the effect of light CO2 produced by methane oxidation in the subsurface layer also decreased CO2 δ13C as shown in Fig. 5. These results suggested that secondary microbial C1 mixes into thermogenic gas. Coal-bearing sediments exist around the Kedr area21,22, and secondary microbial C1 can also form from coal beds30. Hydrate-bound C1 of secondary microbial origin has been only reported at the Alaska North Slope31. This study is another case for it.
Figure 6

A diagram of headspace gases. CO2 δ13C plotted against C1 δ13C, based on the classification of Milkov and Etiope25.

Full size image

Formation process of the sII gas hydrates
As stated before, the crystallographic structure of gas hydrates at the Kedr MV is mainly due to the composition of thermogenic C2 in the volatile hydrocarbons. The concentration of C3, which is one of the sII-forming components, was two to three orders of magnitude smaller than that of C2, because biodegradation occurs and this preferentially reduces C3−5 of n-alkanes19,32, 33. The concentration of n-C4 was smaller than that of i-C4, whereas that of n-C5 was not detected (Table 1). C3 δ13C was around − 10‰, suggesting that light C3 is consumed by microbial activity. Assuming that sediment gas C3+ can be ignored, sediment gas ratio C1/C2 at the study area was 30 ± 17 (mean and standard deviation), and the concentration of C2 was ~ 3%. Such a composition of thermogenic gas is, therefore, considered to be supplied from a deep sediment layer, forming sI gas hydrates composed of mainly C1 and C211,12 in the lake floor sediment.
In the cases where sI gas hydrates plug and block migration pathways, upward fluid flow becomes more focused in other areas16. Once gas supply stops locally, gas hydrates begin to decompose, with the gas dissolving into gas–poor sediment pore water. In the system of C1 and C2, C2 is prone to be encaged in gas hydrate and decreases the equilibrium pressure of mixed-gas hydrate. Therefore, C2-rich gas hydrate forms in parallel with the decomposition of sI gas hydrate. The Colorado School of Mines Hydrate (CSMHYD) program34 showed that C2-rich sII gas hydrate (C2 concentration 17%) forms from mixed gas composed of C1 and C2 (C2 concentration 3%). The C2 concentration of hydrate-bound gas at the Kedr MV was ~ 14%, agreeing fairly well with the results of the CSMHYD program. Such secondary generation of gas hydrates can produce compositions and crystallographic structures that are different from the original crystals. A calorimetric study of synthetic C1 and C2 mixed-gas hydrate revealed that double peaks of heat flow correspond to the dissociation process of C1 and C2 mixed-gas hydrate, suggesting that C2-rich gas hydrate forms simultaneously from dissociated gas and showed that the second heat flow peak correspond to the dissociation of C2-rich gas hydrate18. The PXRD and solid-state 13C nuclear magnetic resonance techniques demonstrated that C2-rich sI gas hydrate forms in the dissociation process of C1 + C2 sII gas hydrate35.
Among twenty hydrate-bound cores in the Kedr area, four cores contained sI only, seven cores had sII only, and seven cores showed sII at the upper layer and sI at the lower layer, as observed at the Kukuy K-2 MV13,16,17. Furthermore, in the cores 2015St1GC15 and 2016St18GC2, gas hydrate structure had sI at the upper and lower layer, and sII at the middle layer. These results suggested that complex gas hydrate layers are composed of sI and sII in subsurface sediments as shown in the schematic illustration in Poort et al.16.
Depth profiles of C2 δ2H of gas hydrate cores from the Kedr MV are shown in Fig. 7. C2 δ2H of hydrate-bound gases varied between − 227‰ and − 206‰, with a grouping around − 210‰. C2 δ2H of sediment gases was also around − 210‰, indicating that C2 δ2H of the original thermogenic gas is − 210‰. As stated above, C2 δ2H of some cores showed low values at their base. Based on the isotopic fractionation of hydrogen in C2 during the formation of sI C2 hydrate36, δ2H of hydrate-bound C2 was 1.1‰ lower than that of residual C2. However, this is too small to explain the wide distribution in C2 δ2H shown in Fig. 7. On the other hand, Matsuda et al.37 reported that isotopic fractionation of hydrogen in C2 is dependent on the crystallographic structure: 1‰–2‰ for sI and ~ 10‰ for sII. Gas hydrates plotting around − 220‰ in C2 δ2H can be explained as a secondary generation of sII from dissociated gas hydrates, of which C2 δ2H was around − 210‰. However, some sII samples showed high C2 δ2H (around − 210‰), whereas some sI samples showed low C2 δ2H (around − 220‰). These results indicated that formation and dissociation processes of gas hydrates produce complicated isotopic profiles in C2 δ2H under non-equilibrium conditions.
Figure 7

Depth profiles of C2 δ2H of hydrate-bound and sediment gases. cmblf, centimetres below lake floor.

Full size image

Characteristics of hydrate-bound gases in sII
C3, i-C4, n-C4 and neo-C5 can be encaged in the larger hexadecahedral cages of sII1. n-C4 and neo-C5 can be encaged using a help gas (e.g. C1) to fill in the smaller dodecahedral cages of sII, because they cannot form pure n-C4 and neo-C5 hydrates, respectively. Figure 8 shows the concentration of C3, i-C4, n-C4, neo-C5 and i-C5 plotted against C2 concentration. The figure illustrates a clear division between sI (3–4%) and sII (14%) C2 concentrations. Data points between C2 concentrations of 5% and 13% were considered to have a mixture of sI and sII. Concentrations of C3, i-C4, n-C4 and neo-C5 had a positive correlation with the concentration of C2, and these concentrations in sII were 1 or 2 orders of magnitude larger than those in sI, suggesting that C3, i-C4, n-C4 and neo-C5 are encaged with C2 in the sII formation process.
Figure 8

Concentration of C3–5 against C2 concentration in the hydrate-bound gases.

Full size image

C3 values of 0.001%–0.01%, ~ 0.0001% of n-C4, and 0.0001%–0.01% of neo-C5 were also detected in sI hydrate-bound gas (Fig. 8), despite these hydrocarbons being unable to be encaged in sI. This can be explained by gases being adsorbed with sediments and gas hydrate crystals, which are then trapped in the grain boundary of polycrystalline gas hydrate crystals, and the gases are encaged if a small amount of sII crystals are present. For example, Uchida et al.38 examined natural gas hydrate retrieved at the Mackenzie Delta (onshore Canada) and detected C3 encaged in sII using Raman spectroscopy, although PXRD results suggested that the sample was sI and the major component of hydrate-bound gas was C1 (more than 99%).
neo-C5 is considered to form from the decomposition of gem-dimethylcycloalkanes derived from the terpenes of terrestrial organic matter39. It is easily enriched by preferential diffusion due to the nearly spherical molecules and its diffusion coefficient, which is higher than that of less branched isomers40. The sII hydrates retrieved at the Kukuy K-2 MV (central Baikal basin) contained 0.026–0.064% of neo-C5 in the volatile hydrocarbons13,14, and those at the Kedr MV had a maximum value of 0.054% of neo-C5 (Supplementary Information Table S1). On the contrary, in the case of natural gas hydrates retrieved at the Joetsu Basin (Japan Sea), neo-C5 was excluded and remained in sediment during the formation of sI gas hydrates from C1-rich gas41. The molecular size of i-C5 is considerably large to be encaged in the large cages of sII. Maximum concentration of i-C5 in the hydrate-bound gases was in several parts per million in both the fields of sI and sII (Fig. 8), indicating that i-C5 is not a hydrate-bound hydrocarbon and adsorbed with gas hydrate crystals and/or trapped in their grain boundary. More

General scheme
We reviewed the literature published until April 2018 on the economic impacts of invasive species. For reasons of feasibility (linguistic skills of the review team, restriction to a reasonable scale of the review), we conducted all searches in the English language assuming that a large body of knowledge (mostly from international peer-reviewed papers and reports) is written in English. The dates of each search process were systematically recorded. We used the following strategy for all repositories (Fig. 1), while also taking into consideration the specificity of their algorithms.
First, a literature search was performed using three online bibliographic sources successively to minimize the risk of omitting relevant materials (Fig. 1, step 1a): ISI Web of Science platform (https://webofknowledge.com/), Google Scholar database (https://scholar.google.com/) and the Google search engine (https://www.google.com/). We carefully composed appropriate search strings that were consensually retained as the most efficient among a set of potential candidates. A decision was taken following preliminary tests based on a handful of relevant articles provided by consulted subject experts on some taxonomic groups (amphibians, reptiles, fishes and ants). Final selection of search strings comprised those considered to have the largest potential to identify key references. Each search string was set to include a combination of two search terms, related to ‘invasive’ and ‘economics’. For both terms, we used a range of synonyms or related words. For example, for ‘invasive’ we used invasi*, invader or exotic; for ‘economics’, we used econom*, cost or monetary. In addition, the search string included exclusion terms to omit mismatches, for example, with studies from the field of medicine that are focused on pathologies or procedures that can be ‘invasive’ for patients. We complemented this search with documents gathered opportunistically (Fig. 1, step 1b). The potentially relevant materials derived from all these sources were combined in a single file and screened for duplicates. Second, retrieved documents were individually assessed at progressive levels (titles, then abstracts, keywords, and finally full text when abstracts were missing; Fig. 1, step 2) based on three criteria. Hence, materials were deemed relevant if (i) they matched with the linguistic competencies of the review team (i.e. written in English, or French where English language was restricted only to the title and/or abstract) for allowing reliable assessment, (ii) they contained at least one cost estimate (studies exclusively providing benefit estimates from direct use or exploitation of invasive species were excluded), and (iii) that this cost estimate is exclusively associated with invasive species (estimates merging non-invasive and invasive species, without the possibility of distinguishing the respective contribution of each group to the overall cost, were excluded). To ensure transparency and validity, each document was checked by two reviewers and in case of a disagreement between assessors, a third reviewer was involved. However, it was often difficult to judge from the topic whether the content of an article was relevant and so consequently many more articles were conservatively kept when final agreement was lacking among assessors.
Finally, relevant materials were scrutinized for data on economic costs (Table 1; Fig. 1, step 3). During this step, additional relevant materials were found as cited by the analysed materials. Obtained cost data were collated in a database and the costs were converted to a common and up to date currency (2017 US$), and then depicted by different descriptors. Categories extracted from relevant materials allow search of the database and data pre-selection to facilitating analysis of costs based on taxonomic groups, geographical areas, impacted sectors, types of costs, or other categories. The reliability of cost estimates and all associated information recorded in the final InvaCost database was systematically checked at least twice, and every ambiguous element was discussed to reach a consensus. We also checked all entries in the database to ensure that there were no obvious duplicate reports (i.e. multiple documents reporting the same cost estimate) or mistakes.
Hereafter, we specifically describe each of the steps made to generate InvaCost.
Literature search
Web of Science
We used the Web of Science (hereafter called WoS) to conduct a search for potentially relevant materials on 7 December 2017 (Fig. 1, step 1a). We applied the following search string: (econom* OR cost OR monetary OR dollar OR euro OR “sterling pound”) AND (invasi* OR alien OR non-indigenous OR nonindigenous OR nonnative OR non-native OR exotic OR introduced OR naturali* OR invader) NOT (cancer* OR cardio* OR surg* OR carcin* OR engineer* OR rotation OR ovar* OR polynom* OR purif* OR respirat* OR “invasive technique” OR carbon OR fuel OR therap* OR vehicle OR cell* OR drug OR fitness OR “operational research” OR banking OR liberalization). The terms were searched in the field code “Topic” which includes title, abstract and keywords, and which also comprises ‘Keywords Plus’ that are generated by WoS through an automatic computer algorithm, based on words and phrases that appear frequently in the titles of article’s bibliographic references and not necessarily in the main text of the article itself. To limit the search to relevant fields of research, we used the function ‘refine’ to exclude subject areas not related to economics and/or invasion biology.
We exported all records (n = 16,875) into an Excel worksheet30 (Table 1) to identify the relevant materials by a two-step procedure. First, we excluded the references identified only based on ‘Keywords Plus’, which were shown to be poor specific descriptors of the content of articles31. We also excluded references identified based on the presence of only a single search term in the topic, as we assumed that words related to both search terms (‘invasive’ and ‘economics’) should be mentioned at least once in the title, abstract and/or keywords of a relevant material. To identify these irrelevant materials within the references collected, we developed a script (see Code Availability) in the R programming language (R v.3.4.3)32. Subsequently, 10,592 references were kept for the next screening step based on the described criteria.
In the second step, the topic of every reference selected was checked manually to ensure potential relevance of its contents. This allowed the elimination of documents incorrectly identified as relevant, such as studies without a true monetary assessment, or those focusing on economic estimates not directly attributable to invasive species only. Finally, 1,333 documents were judged as relevant materials (Table 1) and moved to the final data collation step.
Google scholar
The Google Scholar database is a large source of grey as well as peer-reviewed literature. Nevertheless, we had to modify our approach in order to address inherent limitations of this database as a search tool (see Haddaway et al.33 for a comprehensive analysis). Typically, Google Scholar allows limited Boolean operators (no nesting using parentheses permitted) and search strings are limited to 256 characters. Additionally, only the first 1,000 search results can be viewed and the order in which results are returned is not disclosed. We also wanted to maximize novel information by avoiding too much overlap between the references collected with WoS and those gathered here.
In light of the above, we adapted our search string to generate the most efficient outcome, i.e. sufficiently pertinent to bring the most relevant items to the top of the result list while not unnecessarily large so as to limit the host of non-viewable results. Thus, the following search string was applied on 26 April 2018, using the advanced search facility to search for selected words anywhere in the article (see https://scholar.google.se/intl/en/scholar/help.html#searching for further details): dollars OR euros OR “USD” OR “EUR” OR “NZD” OR “AUD” OR “CAD” OR “GBP” OR “economic cost” OR “economic impact” OR “estimated cost” invasive species. We specified currencies for prioritising materials with monetary data in the top of the resulting list. These currencies were chosen as they were the most often used to express economic costs in the literature collected from the WoS. Nevertheless, any reference evoking economic costs in other currencies was expected to be also captured by some specific combinations of ‘economic’ terms in our search string that we would expect to be mentioned at least once in the full-text of relevant papers. In addition, we included the concomitant presence of ‘invasive’ and ‘species’ terms to restrict the outcomes to papers within the scope of our synthesis. Subsequently, we collected all viewable results (100 pages, n = 992 references of the 668,000 generated), thus going beyond the traditional and arbitrary sample size of first 50–100 results, which is frequently selected in many systematic reviews. We used a web-scraping programme (https://www.webscraper.io/) to extract all the titles’ references returned by the search in an Excel spreadsheet. Because we could not efficiently export the abstract for every reference, we screened them online to assess their potential relevance.
As a result of a search and relevance assessment within Google Scholar, the references, abstracts and specific bibliographic details of 432 documents were added to the sample for further analysis. After excluding duplicates with WoS retrieved references, 310 additional documents were included in the sample as potentially relevant materials (Table 1).
Google
We used the Google search engine to complete the standardised literature search. As when searching with Google Scholar, we took into account specific constraints related to the use of this search engine. Moreover, browsing through Google search results can be overwhelming due to the vast amount of information of highly variable quality. We attempted to implement a search strategy that could allow overcoming these limitations as much as possible. We used the following search string: economic species invasive OR nonnative OR alien OR exotic OR nonindigenous -disease -surgery -fungus -respiratory. We added four exclusion terms (disease, surgery, fungus, respiratory) identified during preliminary tests to restrict the number of irrelevant studies, associated with medical research. We did not use a range of economics-related terms, such as impact or cost, as they returned overly large numbers of mismatches.
The web search was conducted on 8 May 2018 by searching for specified terms within page titles of each document, in order to maximize the likelihood of identifying grey literature. We especially targeted grey literature because searches by the other two platforms mainly led to peer-reviewed publications. We assumed that documents published online by various governmental and non-governmental organisations (NGO), research centres and academic institutes are more likely to contain relevant data than other types of documents such as blogs and catalogues29. Therefore, we restricted our search to the documents located on governmental, academic and NGO webpages to ensure that explicit, traceable and expertise-based information was retrieved. We conducted independent searches for each type of webpage by specifying the type of web extension in the advanced search facility (.gov for governmental,.edu for academic, and.org for organisational webpages).
361 search hits were collected (document name, publishing year and URL of the main website homepage, if available) and stored in the database with the same host of dedicated information (Table 1). If the item analysed was a website homepage, we conducted on-line searches of potentially relevant materials within the website database(s), by filters if available, or by using the search bar with combinations of keywords. Websites that did not contain a database or search bar were searched manually. We then eliminated all duplicates resulting from references being listed on multiple websites, or due to typographical mistakes and/or incomplete records when reporting a reference within different repositories. A total of 119 potentially relevant materials was finally obtained (Table 1).
Targeted collection
Finally, we sourced other potentially relevant materials that did not originate from the above-described processes (Fig. 1, step 1b). On one side, we dedicated specific efforts on gathering cost estimates for particular taxa or areas for which data previously obtained seemed scarce. First, we made sure that some key species were adequately covered; for example, costs associated with invasive mosquito species responsible for much of the burden of mosquito-borne viral diseases worldwide (Aedes aegypti that mainly invaded the intertropical zone from the 15th-17th centuries, and Aedes albopictus for which the global dissemination was more recent34) were searched in a specific way using WoS and PubMed (https://www.ncbi.nlm.nih.gov/pubmed/) repositories (see supplementary file 1 for details on search strings and matching with PRISMA statements). Second, materials were also retrieved following requests to specialists (e.g. Aliens mailing list, https://list.auckland.ac.nz/sympa/info/aliens-l) to bridge gaps identified for Russia and China, two of the five largest countries for which available on-line data were particularly scarce. A typical message first summarized the objectives of our research project and second, requested recipients to provide relevant material and/or suggest further contacts in this regard. On the other side, we also compiled additional materials when establishing the methodology for the project (e.g. when testing different search string combinations at initial stages of the work), from the bibliographic alerts set up by the review team. All 1417 documents obtained from this process were entered in the database, with information on the person providing the document (Table 1;30). Subsequently, 150 documents identified as not previously retrieved were considered relevant for further, full-text screening (Table 1).
Extraction of cost estimates
The Online-only Table 1 comprises all the information of InvaCost that we mention further in this article, using simple quotation marks for ‘Columns’ of the database and italic letters for the different categories within each column. The full-text of each relevant material was scrutinized for any cost estimate that could be incorporated into InvaCost30. The final stage of inclusion/exclusion took place during this data extraction. When the screened documents reported cost estimates by citing sources that were not retrieved by our literature search, whenever possible we assessed the original sources of data in order to better characterize the reported cost. These novel information sources not initially captured by our literature search were then added to the collection list (Table 1). In such cases, we provided information on all documents that were consulted to trace back the original source (‘Previous materials’). In contrast, if no original cost data were found in the cited source, the document was discarded. For all reported costs where the original source was not available or accessible, we emphasized this in a dedicated column (‘Availability’).
Then, we first extracted raw cost data, i.e. how they appear in the material in local currency (‘Raw cost estimate local currency’). When multiple cost estimates were provided for a single instance, we calculated median values (e.g. different cost estimates according to several management scenarios dedicated to the same invasive population) and collated the minimum and maximum estimates provided (columns ‘Min/Max raw cost estimate local currency’). When costs were estimated at different time and/or spatial scales in the same material, we opted to choose – when possible – those estimate(s) that summarise(s) as effectively as possible the figure(s) shown in the study. If such an estimate was not obvious to identify throughout the full-text, we extracted every relevant cost estimate. In these latter cases where several cost estimates were provided in a single study, we also collated the minimum and maximum estimates provided.
Temporal information on the costs were also retrieved: the ‘Period of estimation’ as stated in the material and hence, when possible, the ‘Probable starting/ending year’ of the period of estimation and the ‘Time range’ (year if the estimate is given yearly or for a period up to one year, period if the estimate is given for a period exceeding a year). The ‘Occurrence’ column gives the status of the cost estimate as potentially ongoing (if the cost can be expected to continue beyond the period of estimation) or one-time (if the cost was deemed as unlikely to continue). For cost estimates provided without a clear indication on the timeframe considered, or covering periods shorter than a year, we considered them with a year ‘Time Range’ and a one-time ‘Occurrence’ to avoid the risk of overestimating the duration of collated costs. The ‘Raw cost estimate’– with complementary information on the ‘Time range’, ‘Period of estimation’ and ‘Occurrence’ – can be used to estimate total costs over a given period of time. We then transformed the raw cost estimates to cost estimates per year (‘Cost estimate per year’) by dividing the raw costs with a period ‘Time Range’ by the duration of the ‘Period of estimation’ (obtained from the difference between the ‘Probable ending year’ and ‘Probable starting year’). The raw costs with a year ‘Time Range’ were reported as they are, because they are already considered at the scale of a year.
Description of cost estimates in InvaCost
Each of the cost estimates recorded was characterized by a number of information, including (a) the reference from which the cost was extracted, (b) the taxonomy of the associated species, (c) the spatial and temporal coverage of the study, (d) the typology of each cost estimate and (e) the evaluation of the reliability of the estimation method(s). For most of the variables considered in InvaCost, a non-negligible part of the cost estimates was not attributable to a single existing category due to the lack of precise information provided by the authors or because they simultaneously belong to multiple categories. In such cases, we respectively reported them as either Diverse/Unspecified or as slash-separated lists of categories (e.g. Artiodactyla/Carnivora for the ‘Order’).
Details about the nature of the information retrieved as well as the choices made to characterize each cost are synthesized in Online-only Table 1:
(a) We provided bibliographic information on each reference (e.g. ‘Reference title’, ‘Authors’, ‘Publication year’). Others specific details (e.g. abstract, journal, download link) are given in a dedicated file30 with which the columns ‘Repository’ and ‘Reference ID’ of InvaCost allow correspondence of information.
(b) We normalised and harmonised all taxonomic information on the invasive species (‘Kingdom’ to ‘Species’ level) using the GBIF.org Backbone Taxonomy35. At this stage, spelling and other taxonomic errors were corrected. While each cost extracted was generally associated with a single invasive alien species, in some cases the data was related to multiple species without the possibility of disentangling species-specific costs. In this case, we mentioned either all species concerned if explicitly indicated by the author(s), or Diverse/Unspecified if not.
(c) We dedicated seven columns to describing the impacted area according to its environment (terrestrial and/or aquatic habitats), the temporal extent as mentioned earlier (e.g. ‘Period of estimation’, ‘Time range’) and the spatial coverage from the ‘Geographic region’ (e.g., Central America, South America, Oceania-Pacific Islands) – rather than the official continent for better accuracy – down to the exact site (‘Location’) when possible. Each area was related to its country of attachment, leading to some mismatches between the ‘Geographic region’ and ‘Official country’ columns due to the existence of countries with non-contiguous overseas territories. For instance, costs found from invaders in La Réunion (a French oversea department) were attributed to Africa as ‘Geographic region’ and France as ‘Country’, while France obviously belongs to European continent.
(d) We characterised the typology of each cost mainly based on the following descriptors. The ‘Implementation’ at the moment of the cost evaluation states whether the reported cost was observed (i.e. cost actually incurred by an invasive species within its invasive distribution area) or potential (i.e. not incurred but expected cost for an invasive species beyond its actual distribution area and/or predicted over time within or beyond its actual distribution area). The ‘Acquisition method’ provides information on how the cost data was obtained, i.e. report/estimation directly obtained or derived (using inference methods) from field-based information, or extrapolation relying on computational modelling. The ‘Impacted sectors’ indicates which activity, societal or market sectors were related to the cost estimate (see Table 2 for details). The ‘Type of cost’ ranges from the economic damages and losses incurred by an invasion (e.g. value of crop losses, damage repair) to different levels of means dedicated to the management of biological invaders (e.g. control, eradication, prevention).
Table 2 The different market and/or activity sectors mentioned in InvaCost.
Full size table

(e) Lastly, we evaluated the level of ‘Reliability’ of the methodology reported by the authors to provide cost estimates (Fig. 2). Prejudging the relevance of each cost estimate is not straightforward and could suffer from a high level of subjectivity. Here, we rather aimed to evaluate in the most objective manner whether the approach used for cost estimation was documented and traceable. Hence, materials that could not be accessed for full-text investigation were conservatively considered as of low reliability. Alternatively, each cost estimate recorded from any accessible material was qualitatively assessed as of high or low reliability following a procedure depending on the ‘Type of material’ analysed (peer-reviewed article or grey material; Fig. 2). Peer-reviewed articles and official documents (e.g. institutional or governmental reports) are likely validated by experts before publication. We assumed therefore that all cost estimates collected from these materials may likely be of high reliability. Conversely, for grey materials other than official reports, the attribution to one or other of these categories (high vs low reliability) was based on specific analysis of each cost estimate. We checked whether the method estimation was fully described, independently of its comprehensiveness, i.e. if the original sources or potential assumptions were properly documented or justified, and/or the calculation methodology was explicitly demonstrated. Here, we opted for a conservative strategy that might be not optimal, as depending mostly on the nature of the publication.
Fig. 2

Decision tree approach for assessing the reliability of the method used for estimating each cost. The colour of the boxes indicates which decision was taken: green when material was deemed as of high reliability, red when material was deemed as of low reliability, blue when taking any decision needs further investigation. The intended purpose of this process was not to evaluate the quality, relevance or realism of the studies performed for providing cost estimates, but rather to assess if the methodology (i) has been reviewed and validated by peers or experts prior any publication, or (ii) if not, whether this methodology was clearly stated and demonstrated.

Full size image

Beyond the factual elements included in the descriptors from (a) to (c), those presented in (d) and (e) (to which we can add the descriptor ‘Spatial scale’) are the result of a conceptual and analytical framework created based on our own experience. This experience was gained when collecting and getting acquainted with the diversity and complexity of situations one can find behind the “economic costs” linked to biological invasions, as well as the strategies used for estimating them. We think that the different subcategories identified therein (e.g. observed vs potential costs within the descriptor ‘Implementation’) should not be aggregated to limit potential confusions in future analysis. Also, we acknowledge that the possible sub-categories of these descriptors might be improved and adapted according to the scope of future analyses made using InvaCost. We are convinced that the descriptors thus defined and categorised may strongly help in this perspective.
Standardisation of cost data
Using definitions, data and indicators provided by the World Bank Open Data and the Organisation for Economic Cooperation and Development (OECD), we expressed all retrieved costs (raw costs and costs per year) in US dollars (US$) for the year 201730 using a multi-step procedure. We provided here two ways for standardising cost estimates according to the conversion factor: one based on the market exchange rate (local currency unit per US$, calculated as an annual average), and another based on the Purchasing Power Parity (PPP, local currency unit per US$, calculated as an annual average) that is the rate of currency conversion that standardises the purchasing power of different currencies by eliminating the differences in price levels between countries. Opting for one strategy or the other for further investigation or discussion is beyond the scope of this paper and will befall on the author(s) of future analyses made using InvaCost.
We first converted the cost estimates from local currencies to US$, by dividing the cost estimate with the official market exchange rate (https://data.worldbank.org/indicator/PA.NUS.FCRF?end=2017&start=1960) corresponding to the year of the cost estimation (‘Applicable year’, that is the year of the ‘Currency’ value, but not necessarily the year of the cost occurrence). The cost obtained in US$ of that year was then converted in 2017 US$ using an inflation factor that takes into account the evolution of the value of the US$ since the year of cost estimation. The inflation factor was computed by dividing the Consumer Price Index (CPI, which is a measure of the average change over time in the prices paid by consumers for a market basket of consumer goods and services; https://data.worldbank.org/indicator/FP.CPI.TOTL?end=2017&start=1960) of 2017 by the CPI of the year of the cost estimation.
As an alternative, we also converted costs to 2017 US$ value based on PPP instead of the classical market exchange rates in the initial conversion step. PPP values were primarily collected from data provided by the World Bank (https://data.worldbank.org/indicator/PA.NUS.PPP?end=2017&start=1990), or by the OECD (https://data.oecd.org/conversion/purchasing-power-parities-ppp.htm) when information was not retrievable through the World Bank database. For this purpose, we had to deal with published costs that were expressed in currency that was different from the country where the costs were estimated (e.g. published cost in African countries expressed in US or Canadian $). Thus, prior to using PPP as a conversion index, we had to perform a retro-conversion by multiplying the original cost estimate by the official market exchange rate (local currency unit per currency unit used). For PPP-based standardisation, it was not possible to perform the process for all cost estimates as PPP data do not exist for all countries and/or specific periods (we mentioned NA in the database when such information was missing).
In summary, we used the following formula to convert and standardise each cost estimate:

$${C}_{e}=left({{boldsymbol{M}}}_{{boldsymbol{V}}}/{{boldsymbol{C}}}_{{boldsymbol{F}}}right),times ,{{boldsymbol{I}}}_{{boldsymbol{F}}}$$

with Ce = Converted cost estimate (to 2017 US dollars based on exchange rate or Purchase Power Parity), MV = Cost estimate (either the ‘Raw cost estimate local currency’ extracted from analysed paper or the ‘Cost per year local currency’ transformed by us), CF = Conversion factor (either the official market exchange rate or the purchasing power parity, in US dollars), IF = Inflation factor since the year of cost estimation, calculated as CPI2017/CPIy with CPI corresponding to the Consumer Price Index and y corresponding to the year of the cost estimation (‘Applicable year’).
We thus provided four columns with the raw cost estimates or the cost estimates per year, expressed in 2017 USD based on the exchange rate or PPP.
Data summary
InvaCost currently contains 2419 cost estimates (1215 from peer-reviewed articles, 1204 from grey materials), collected from 849 references, of which 1769 estimates were deemed as of high reliability. In total, twenty currencies are reported in our database, the majority being US dollars, n = 1348 cost estimates. Not all cost estimates were successfully converted to 2017 US$ as (i) conversion data from official sources are available only since 1960 (cost estimates range from 1945 to 2017 in InvaCost) or simply not found for some years and countries, and/or (ii) cost data are sometimes simultaneously associated with several countries, constraining the PPP-based standardisations. Hence, respectively 2416 and 2126 estimates were successfully converted using market exchange rates and PPPs. Cost estimates are either direct reports/estimations (n = 2127) or values gathered from extrapolative computations (n = 292). At a taxonomic level, these estimates are associated with 343 species belonging to six kingdoms (Animalia, Bacteria, Chromista, Fungi, Plantae, Riboviria). InvaCost has global coverage (90 countries) and includes continental, insular and overseas territories. Data are associated with terrestrial as well as aquatic (freshwater, brackish and marine) environments. Costs were estimated at different spatial scales (continental (n = 35), country (n = 1111), global (n = 17), intercontinental (n = 9), regional (n = 67), site (n = 836), unit (n = 329)). The Table 3 summarises quantitative data and information reported in InvaCost for each geographic region considered (see also Supplementary file 2).
Table 3 Quantitative summary of information recorded in InvaCost according to the ‘Geographic region’ of the cost estimates.
Full size table

Possible applications
InvaCost is expected to help bridge the gap between a growing scientific understanding of invasion impacts and still inadequate management actions. This work is thus in line with the aims of a panel of decisions recently adopted by the Convention on Biological Diversity (Decision XIII/13, https://www.cbd.int/doc/decisions/cop-13/cop-13-dec-13-en.pdf) advocating the incorporation of invasion science knowledge into management planning. In addition to offer unique opportunities for future research, InvaCost will provide a strong quantitative and evidence-based support for impacts of invasive species reported in other databases such as the Global Register of Introduced and Invasive Species (GRIIS)20, helping refine information in this database. Also, invasive populations recorded in InvaCost but data deficient in the GRIIS should be ultimately classified in that database.
Additionally, InvaCost could be considered as another data-based component, adding novel and significant information on invader impacts categorised by the Socio-Economic Impact Classification of Alien Taxa (SEICAT)36. The latter is a classification system, applicable across a broad range of taxa and spatial scales, providing a consistent procedure for translating the broad range of measures and types of impacts into ranked levels of socio-economic impacts, assigning alien taxa on the basis of the best available evidence of their documented deleterious impacts. Quantitative support provided by InvaCost will strongly contribute to impact classification. Ultimately, integrating data from these diverse sources could allow a complete description of the overall impacts of biological invasions at regional and global scales.
Caveats and directions for further database improvement
Rather than claiming exhaustiveness of data collated, we highlight that InvaCost should be considered as the most current, standardised, accurate and globally representative repository of various economic losses and expenditures documented for the largest possible set of invaders. We are aware that our database can be improved in at least three ways.
First, InvaCost mostly does not include publications and reports not yet available in electronic format and/or using non-English language, leaving open the possibility of increasing data comprehensiveness and limiting potential biases. Indeed, local reports as well as research results from some countries (e.g., China, Russia) are likely to be published in non-English language37. Again, accessing grey literature is challenging as it is not systematically digitalised and/or included in well-curated bibliographic databases29. We strongly encourage future users of InvaCost to help gathering this currently unreachable information when possible. Furthermore, some mistakes might have occurred despite our best efforts when constructing InvaCost. In this regard, we advocate for regular public updates of InvaCost in order to improve it both quantitatively (by adding currently inaccessible or missed information) and qualitatively (if errors are identified).
Second, as the distribution and impacts of invaders are inherently dynamic for a number of reasons38, InvaCost should further consider the status of the species recorded for their economic impacts in order to improve both the relevance and the usefulness of the database. As an illustration, InvaCost likely includes invasive populations currently extirpated from particular areas after successful eradication campaign(s) as well as those still established but for which impacts are locally reduced as a result of management efforts. Attempting to obtain and integrate such information into InvaCost was beyond the scope of this work. Nonetheless, it should be reciprocally beneficial to establish connections between InvaCost and other databases such as the GRIIS that provides a harmonised, open source, multi-taxon database including verified information on the continued presence of introduced and invasive species for most countries20. In light of such additional information, the value of InvaCost will be its application for policy purposes, such as identification of exotic invaders that are currently associated with economic losses in particular areas. Also, crossing information between databases may allow the refinement of the descriptor ‘Spatial scale’ we propose here.
Third, we would recommend, for a future updated version of InvaCost that would require screening back all the materials, to improve the ‘Acquisition method’, ‘Implementation’ and ‘Reliability’ descriptors, to pay attention to the specificity of “avoided costs” and to create a new descriptor for ‘non-market values’. We detail these possibilities below.
Improving descriptors
An improved version of the ‘Acquisition method’ could lead to a subdivision of the extrapolation category into spatial, temporal and spatio-temporal extrapolation. This would allow simultaneous refinement from the currently binary ‘Implementation’ descriptor (observed vs potential) into several levels of certainty regarding the incurred cost (e.g. taking into consideration the temporality (past/current or predicted) of the onset of the cost and of the status of the invasive species in the study area). The next step for deeming the ‘Reliability’ of the cost estimates recorded in InvaCost would consist of assessing the repeatability of the methodology used, by adapting the approach previously developed by Bradshaw et al.14. The latter evidenced that assuming the reproducibility of published methods should not rely only on the nature of the materials and recognized the qualitative nature of the procedure, although applying this approach to InvaCost was constrained by the large sample size and high diversity in our database (Bradshaw et al.’s study focused on a single taxonomic class). Also, because InvaCost involves several collaborators and potential future contributors, consistent and objective criteria should be further defined to cope with the large array of materials, methods and situations encountered.
Avoided costs
Introducing certain actions against biological invasions leads to avoided costs. Such avoided costs are sometimes evaluated, for instance to examine the relevance of different potential actions or to assess the effectiveness of an action that was taken. However, avoided costs cover a great variety of situations and require a careful consideration for future analysis, even if they do not have to be analysed separately from the other economic costs gathered in InvaCost. For instance, in the case of hypothetical actions, avoided costs can be considered as minimum estimates of the “real” costs (if they are unknown). However, in the case of completed or planned actions, the reported data should be the original costs (if known) minus the avoided costs, because the latter do no longer exist. Some avoided costs are probably already included in InvaCost but they are likely underestimated because keywords such as “savings” or “benefits” were not included in the search strings. Also, even if they are sometimes mentioned as “benefits” in the literature, care should be taken not to confuse these avoided-costs with the benefits incurred by direct use or exploitation of invasive species. The latter have been ignored in InvaCost since they were relatively few (and beyond of the scope of this database), but might constitute a twin project.
A new ‘Non-market values’ descriptor
The means dedicated to preventing or managing an invasion (e.g. manual removal of invasive plants) and certain economic losses and damage due to an invasion (e.g. the value of crop losses or the repair costs of damaged infrastructures due to an invasive insect) are observable on markets. However, some costs are not observable on markets but can be translated in monetary terms using several valuation methods – for instance, the willingness to pay for the conservation of a native species that is impacted by an invasive species is considered as the value given by a group of people to preserving the native species (i.e. the value that would be lost if this native species was impacted). We recognize the importance of informing the public about “non-market values”, as giving an economic value to ecosystems or biodiversity can be a way of recognising and taking them into account in public decision-making processes39, but attention should be paid to the issues linked to their assessment40,41. Among others, the different methods for assessing non-market values do not necessarily capture the same aspects of the values, so the resulting estimates might be different. Moreover, the very principle of giving a value to “benefit from nature” through economic valuation is not necessarily acknowledged by the entirety of scientific and civil communities39,42. For future analysis, the ‘non-market values’ should not be systematically aggregated with the other economic costs gathered in InvaCost. It is to note that while some non-market values are probably already included in InvaCost within the losses and damage ‘Type of cost’, the loss of non-market values is probably largely underestimated in the database because they were not the primary focus of InvaCost and therefore the related keywords were not included in the search strings.
These possible ways of improvement call for completion and/or refinement of existing entries as well as integration of newly published or acquired data by future contributors in InvaCost, with the aim to consolidate its long-term relevance (cf. Usage Note paragraph). More

Tree-ring data
We used tree-ring records from over 210,000 trees of 110 species, distributed globally in habitats ranging from the tropics to the Arctic region over more than 70,000 sites (Supplementary Fig. 1, Supplementary Table 1). The largest publicly available data source from which we used data is the International Tree-Ring Data Bank (ITRDB, https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring). These were complemented with other datasets to maximize the number of records for each species and to fill in spatial gaps. A particularly large tree-ring dataset used in our analyses is the National Forestry Inventory data from the Ministère des Forêts de la Faune et des Parcs from Quebec, Canada52,53 (hereafter NFI-Quebec). This data consists of a complete set of ring-width data from 156,711 trees from 79.381 sites across the province of Quebec. Field tree-ring data were collected according to specific standard protocols52,53, which consisted of selection of up to nine trees in each plot, with 3–5 trees ( >91 mm diameter at breast height, DBH) randomly selected, 1–2 selected from the largest trees, and 1–2 from trees closest to the mean tree diameter of the plot52,53. From selected trees one core per tree was collected. Tropical tree-ring data were compiled from the ITRDB and from unpublished and published records26,27,54,55. For species with larger sample sizes, we distinguished between tree-ring data from trees that died naturally before the moment of sampling and trees that were alive at the moment of sampling, allowing us to test the assumption that living tree ages can be used to estimate trees natural lifespans (see section “Trade-off estimates and assessment of possible artefacts”). Part of these dead tree data were obtained from the ITRDB by selecting tree-ring records of which the last measured ring width was dated to before AD1900. We assumed that most of these trees must have been dead at the time of sampling as no records were collected before 1900. In addition, we compiled published dead tree data from refs. 21,56, and used subfossil tree-ring data from refs. 57,58. Supplementary Table 1 provides an overview of the datasets, and full details of each dataset are available online as supplemental info.
Various data controls and selection procedures were used to assure high confidence in our dataset. Where possible we tried to identify duplicate records, i.e., multiple cores taken from the same individual tree. This is only a problem for ITRDB, but not for the NFI-Quebec dataset where only one core per tree existed, or for datasets from co-authors. We thus merged ITRDB records that had identical ID′s except for the last character of their ID (e.g. 01a, 01b, or ID1-1, ID1-2, etc). From the ITRDB, we only used species for which we could obtain data from a minimum of 3 different sites with at least 20 records each, and only selected species that had a minimum total of 100 separate ring width series. We excluded those sites from the ITRDB that showed relative even age structures, and are thus unlikely to represent old-growth populations that provide robust estimates of trees’ maximum lifespans. To this end, we calculated for each ITRDB site the coefficient of variation in tree ages (CVAge = StandDevAge/ MeanAge × 100) and excluded sites with a CVAge lower than 10%. A large subset of ITRDB-data from 46 species has previously been inspected for data quality by co-author S. Voelker59. In this subset of data, each ring width series was manually re-aligned by cambial age (i.e., ring number from pith), providing more reliable estimates of tree ages. For the datasets that were not acquired from the ITRDB, we used slightly different criteria. From NFI-Quebec, we used all available sites, excluding those that were classified with evidence of recent management (commercial thinning or clear-cutting) and where fire or insect disturbances destroyed more than 25% of the forest cover.
For estimates of species-level early growth rates and lifespans (cf. Fig. 1a), we included only species with a minimum of 30 records, as lower sample size is unlikely to provide good approximations of tree life spans. This resulted in the inclusion of 110 species, with a median sample size of 305 trees and 12 sites per species. To assess within-species relationships between early growth and tree lifespan, we included 82 species. As a general rule, we included only species with more than a total of 150 trees and from at least 3 sites. About half of our species had more than 300 tree records (see Supplementary Information).
To assess what minimum sample size is needed to get a representative estimate of the true maximum age of a species or a site, and to evaluate how sample size affected estimates of trade-offs between early growth and longevity, we randomly resampled 500 times varying sample sizes—from 25 to 600 trees—from a subset of 11,752 Picea mariana trees from NFI-Quebec sites located north of 50.7°N. Comparison of the maximum ages of these random subsets of trees with the true observed maximum ages shows that a sample size of 100 trees results in 99.4% of the cases in maximum age estimates larger than the 95th percentile of the original dataset, and in 67.2% of the cases in ages larger than the 99th percentile original age (Supplementary Fig. 2a, b). As more than 70% of the species had at least 100 trees we thus assumed that for most species, the estimates of their lifespan were close to true lifespans. We used this same approach to assess how sample size affected the estimate of the trade-offs (i.e., estimation of the negative exponential decay constant; see next section for details). This analysis showed that sample size of 300 trees (corresponding to median sample sizes for trade-off analysis), leads to mean errors in the estimated slope of 12% (Supplementary Fig. 2c). Thus, for most species we achieve relatively accurate estimates of the trade-off strength. Low sample sizes for some species will nevertheless result in small errors of the mean slope, but we expect that positive and negative errors will cancel out against each other. Indeed, we do not observe a specific bias towards over- or under-estimation for low sample sizes, as the mean exponential decay constant for a simulated sample size of 150 trees is very similar to that observed (i.e., −0.409 versus −0.399).
Trade-off estimates and assessment of possible artefacts
The strength of the trade-offs between growth and tree lifespan was assessed for each species using a 95th quantile regression between mean early growth rate and the natural logarithm of age using the QUANTREG package in R60, as

$$begin{array}{*{20}{c}} {{{log}}(A_{(95{mathrm{th}},{mathrm{quantile}})}) = a + b cdot overline {{mathrm{RW}}} } \ {{mathrm{or}}} \ {{mathrm{Lifespan}} approx A_{left( {95{mathrm{th}},{mathrm{quantile}}} right)} = {mathrm{exp}},left( {a + b cdot overline {{mathrm{RW}}} } right)} end{array}$$
(1)

where A is age of the tree, (overline {{mathrm{RW}}}) is the mean ring width over the first 10 years. The constant b describes the negative exponential decay constant (i.e., exponential rate of decrease of tree lifespan with increasing early growth rate). This quantile regression fit results in similar estimates of the maximum ages of trees as the 95th percentile ages in binned early growth rate categories (see Supplementary Fig. 3a–c). Note that in contrast the maximum diameter does not vary strongly between slow and fast-growing trees (Supplementary Fig. 3d). We chose a relative short period, the first ten years, for estimating early growth as our study included some relative short-lived species. Previous studies have found similar results when using longer periods (50 years)56, and we expect no substantial difference using different early growth periods as tree growth is usually strongly auto-correlated in time61.
To assess trade-off strengths within species, we calculated the mean decay constant (b, Eq. 2) for each species using relative age, A/max(A), and relative mean early ring width, (overline {{mathrm{RW}}})/max((overline {{mathrm{RW}}})). Maximum of A and (overline {{mathrm{RW}}}) are species level maxima. The mean slope calculation across all species was weighted by the cube root of the sample size to account for the large differences between species in sample size, and confidence of the trade-off estimates.
While these relationships suggest true trade-offs, they may also be affected (or even driven by) the approaches or analytical methods used here. In particular, we here evaluate the effect of the following four possible artefacts on our results; (1) the use of living trees to estimate tree lifespans, (2) effect of recent growth increases on early growth-age relationship, (3) effects of pith offsets and wood decay on early growth-age relationship, (4) sampling artefacts, such as disproportionate selection of large trees.
(1)
Use of living trees: our analysis includes mostly trees that were sampled when still alive, and may thus not be representative for the true lifespan trees may achieve. To assess to which degree use of living trees may affect our results, we analyse and compare the trade-off strengths of trees that died before 1900 to living trees for 12 species with sufficient data availability (minimum of 150 dead and 150 living trees). As the slopes of dead and living trees do not differ significantly (Supplementary Fig. 5f, g, paired t-test exponential decay coefficient, t = −0.1095, p = 0.915, n = 12), we conclude that 95th quantile regressions on living trees can be used to approximate tree lifespan.

(2)
Effect of recent growth increases: recent growth stimulation of trees due e.g. to CO2 fertilisation, warming in higher latitudes, and/or nitrogen deposition, may result in observation of a trade-off. This is because recent increases in growth will lead to higher early growth rates for young trees compared to old trees, resulting in a negative relationship between early growth and tree age. The comparison of trade-off strength of dead versus living trees provides strong evidence that this effect does not drive the trade-off. In addition to this, we use a data driven forest simulation (see section “Examining the effects of growth stimulation on forest dynamics”) to assess how growth increases affect estimations of the trade-off strength. In this simulation, we used the actual tree-ring data to simulate realistic growth increases of Picea mariana tree-ring trajectories in response to high latitude warming. By sampling from these trajectories at the end of the growth increase period (i.e., year 350), and in a period without any recent growth increases (i.e., year 600, see Fig. 3e), we establish that growth increases result in only a small over-estimation of the trade-off, decreasing the exponential decay coefficient from −0.37 to −0.44 (see Supplementary Fig. 9c). Thus, it is unlikely that recent growth stimulation is the cause for the negative relation between early growth and tree lifespan.

(3)
Pith offset: tree-ring data, especially those acquired from ITRDB, may miss the innermost sections due to incomplete cores, decayed centres, or imperfect increment borer alignment. Missing rings will result in underestimation of tree ages and inaccurate early growth rates estimation and could thus affect the estimated relationship between early growth and lifespan. However, ring widths in most species decrease with tree age and size17, and even trees showing constant wood production with age, will show decreasing ring width because of geometry. Thus early growth in these samples will underestimate true growth rates and would most likely weaken the observed trade-off, rather than strengthening it. A comparison of species present in both the NFI-Quebec and the ITRDB datasets confirms this. The NFI-Quebec dataset was less affected by pith offset problems, as the trees were carefully screened and trees with substantial differences between cumulative ring widths and field diameters were excluded. Yet, we find that slopes were more negative for NFI-Quebec compared to ITRDB (mean b of −0.25 for Quebec vs. −0.10 for ITRDB, two-sided paired t.test, t = 2.49, p = 0.047, n = 7 species) and pith offsets thus do not explain the relationship. This comparison also shows that estimates of the strength of the trade-offs between early growth and longevity inferred from ITRDB data are probably conservative, as the Quebec data can be considered to be of higher quality, and were collected according to standard protocols. In contrast, data from the ITRDB may contain incomplete series and were collected for unknown purposes, and these issues probably weaken trade-offs in the ITRDB.

(4)
Sampling biases: one potential bias in our dataset may arise due to the tendency of tree-ring studies to sample predominantly large trees in the field (i.e., big tree selection bias62,63,64). This may result in a negative relationship between early growth and tree age, as young slow-growing trees tend to be underrepresented in the tree-ring sample (i.e., have not reached the field minimum size threshold yet), compared to fast-growing young trees that are much larger, and therefore more likely to be sampled. This effect would reduce the number of trees with slow early growth and young ages in the tree-ring sample (i.e., trees in the lower left-hand corner of the early growth-lifespan graphs, cf. Supplementary Fig. 6a), and results in overestimation of the 95th percentile age estimates for slow-growing trees. Our approach to estimate to which degree this bias affected our estimates of growth-lifespan trade-offs was as follows. We first used the tree-ring NFI-Quebec data of Picea mariana, combined with plot data from Quebec to reconstruct a new artificial tree-ring dataset with a size frequency distribution identical to the population size frequency distribution for this species in Quebec (Supplementary Fig. 6b). For each tree of Picea mariana sampled for their tree-rings we know the early growth rate and age, and also the complete diameter- and age-trajectory up to the year of sampling. From these data, we resampled for each size class (in bin widths of 2 cm) the same number of trees as that observed in the field. By doing this we filled in the lacking data of trees smaller than 91 mm, and created a new artificial tree-ring dataset that had an identical size structure to that observed in the field. We know the mean growth rate over the first ten years and the age at which each individual tree reached the diameter of their respective size class, and could thus reconstruct the early growth rate versus tree age graphs for the full population, including the smaller size classes which were missing from our original tree-ring sample. We then compared the early growth -lifespan relationship for the complete population to that of the trees larger than 91 mm, mimicking the NFI-Quebec field collection protocol. This shows that the exponential decrease is marginally larger (b = −0.505 compared to −0.470 for trees >91 mm) and that the intercept is lower (159 years compared to 220 for trees >91 mm) when sampling all trees compared to only trees with diameters >91 mm (Supplementary Fig. 6). Hence, the use of a minimum size threshold (91 mm) in the NFI from Quebec results in a slight underestimation of the trade-off (by ~7%) for the Picea mariana dataset. We also resampled from this artificial dataset the 10% largest trees, to mimic a hypothetical standard tree-ring sampling scenario that only samples the largest trees. Such a sampling scenario resulted in a decay constant of −0.432, thus again causing a small underestimation of the true trade-off. This simulation proves that the trade-off is not a result of a sampling bias.

Possible environmental drivers of the trade-offs
We evaluated whether the observed trade-off between early growth and tree lifespans could be caused by covariance of growth and lifespan with climate, soil or competition. Temperature variation for example reduces tree growth and lifespan in various species (cf. Fig. 2a, b). To this end, we calculated site-level mean early growth rates and the maximum tree age for a set of species covering different geographic regions (North America, Europe and Quebec). For Quebec, we combined multiple nearby locations to obtain a minimum of 30 trees per site, as sample sizes were low for each location. Site-level mean annual temperature and precipitation was obtained from WorldClim65. We then assessed for nine different species the effect of temperature and precipitation on site-level mean early life growth and maximum tree age using major axis regression from the package smart-366. These analyses confirm that early life growth is positively related to temperature for all nine species studied, and that lifespan decreases significantly with temperature for seven out of nine species (see Fig. 2a, b). Using linear mixed effect models with species as random factor (nlme-package-R67), we find that across all nine species, mean early life growth increases on average by 0.11 mm for each degree temperature increase, while lifespan decreases by 13 years for each degree temperature increase. Precipitation has no significant effect on early life growth or tree lifespan.
To disentangle whether lifespan decreases are a direct effect of temperature increases, or due to increases in early life growth, mixed effect models were run for all nine species that simultaneously included temperature and mean early life growth rate as explanatory variables for variation in tree lifespan. To account for species differences in growth and age, we used relative mean early ring width ((overline {RW})/max((overline {RW}))), and relative maximum age (A/max(A)), and used species as a random factor with random intercepts for both early life growth and temperature. This analysis shows that mean early life growth is a stronger predictor of tree lifespans than temperature (t-value early growth = −7.2, p  , D_0} end{array}} right..$$
(4)

Here b = 0.025, k = 3 × 10−11 and D0 = 91 mm, which is the minimum sampling diameter used for the NFI-Quebec. The rationale to set the mortality below D0 to zero in this model is that the trees sampled for tree rings did all survive to this diameter (as only trees with diameters >91 mm were cored). Thus, only by setting mortality to zero for trees , A_0} end{array}} right.,$$
(5)

with a = 0.021 and b = 0.0000015. Analogous to the diameter-dependent mortality model, we set mortality to zero for trees younger than 74 years (A0), which is the age at which the average Picea mariana tree reaches 91 mm in diameter (D0).
We next created a 600 year sequence of annually seeded, 1250 member, tree cohorts. Each member of a cohort is a randomly selected tree diameter growth trajectory derived from the tree-ring cores, extended in time to an age of 500 years. Short growth trajectories were extended by using the mean growth of the 10 oldest trees that had a similar ring width in the first ten years of growth. To this end all early mean ring width was grouped into six equal early ring width classes.
All trees of this so-constructed tree cohort set live, by construction, exactly 600 years. Realistic age structures were realised by sequentially (year-on-year and tree-by-tree) assigning death to trees where a random number generator identified those individuals  smaller than μ(D) or μ(A). To test the realism of this procedure, we compared the predicted and observed tree age versus early ring width relationship—or i.e. the growth-longevity trade-off. The slope of the relationship for the diameter-dependent mortality model μ(D) is very close to observed (Supplementary Fig. 9a), and is thus a realistic representation of the observed mortality process and justifies its use to examine the effect of a growth stimulation on standing stocks. In contrast, we find that the age-dependent mortality model μ(A) does not result in a significant trade-off between early growth and tree lifespan (Supplementary Fig. 9b), providing an ideal comparison for models that fail to incorporate the observed trade-offs.
To mimic growth stimulation, we boosted growth of trajectories from year t0 = 300 year onwards of the 600 year sequence of cohorts, while exposing the trajectories over the entire 600 year period to the mortality algorithm just described. We stimulated growth rate, (overline {{mathrm{RW}}}) (mm year−1) from year t0 = 300 year onwards according to

$$overline {{mathrm{RW}}_{stim}} left( t right) = overline {{mathrm{RW}}} left( t right) cdot exp left( {lambda left( A right) cdot delta Tleft( t right)} right),$$
(6)

and

$$delta Tleft( t right) = left{ {begin{array}{*{20}{c}} {frac{{overline {dT} }}{{dt}} cdot left( {t – t_0} right),} & {t – t_0 ,150 yrs) as younger trees are more sensitive to temperature increases than older trees (Supplementary Fig. 9d), and use an exponential model of the form

$$overline {{mathrm{RW}}} left( {A,,T} right) = a cdot {mathrm{exp}}left( {lambda left( A right) cdot delta T} right),$$
(8)

where a is a constant, λ(A) is the exponential increase rate for age class [A, A+ δA] and δT = (T − 0) (°C). We then used the exponential growth increase with temperature for each age band, λ(A), to estimate the relationship between tree age and λ(A) (Supplementary Fig. 9e). In all cases the age modulation of the stimulus is

$$lambda left( A right) = left{ {begin{array}{*{20}{c}} {0.0000132 cdot A^2-0.00291 cdot A + 0.229,} & {A , < ,135} \ {0.07,} & {A , ge ,135} end{array}} right..$$ (9) Finally, we compared the effect of growth stimulation on forests dynamics against simulation without growth increases (baseline) for the two different mortality models. We also performed one simulation where we multiplied the full growth series by 2, as a representation of the effect of growth stimulation on a faster-growing species. Finally, we evaluated for each simulation the mean ring width growth, the stem mortality rate, the age of the largest (75th percentile) trees that died and the change in the total basal area stock over time for the full population. For the stem mortality and the basal area stocks, we calculate and present the change in dynamics of the growth stimulation scenario relative to the baseline scenario without growth stimulation. All analyses and simulations were performed using R-studio, version 0.99.90370. Maps in SI figures were produced using the ggmap function from the R-package ‘ggplot2’71 More

Affiliations

State Key Laboratory for Biology of Plant Diseases and Insect Pests, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China
Kris A. G. Wyckhuys & Yanhui Lu

Fujian Agriculture and Forestry University, Fuzhou, China
Kris A. G. Wyckhuys

University of Queensland, Brisbane, Queensland, Australia
Kris A. G. Wyckhuys & Michael J. Furlong

Chrysalis Consulting, Hanoi, Vietnam
Kris A. G. Wyckhuys

Zhejiang University, Hangzhou, China
Wenwu Zhou

CABI, Wallingford, UK
Matthew J. W. Cock & Frances E. Williams

USDA ARS, Maricopa, AZ, USA
Steven E. Naranjo

Secretariat of the Pacific Community SPC, Suva, Fiji
Atumurirava Fereti

Authors
Kris A. G. Wyckhuys

Yanhui Lu

Wenwu Zhou

Matthew J. W. Cock

Steven E. Naranjo

Atumurirava Fereti

Frances E. Williams

Michael J. Furlong

Corresponding authors
Correspondence to Kris A. G. Wyckhuys or Yanhui Lu or Wenwu Zhou. More

1.
Perry, A. L., Low, P. J., Ellis, J. R. & Reynolds, J. D. Climate change and distribution shifts in marine fishes. Science 308, 1912–1915 (2005).
CAS  Google Scholar 
2.
Dulvy, N. K. et al. Climate change and deepening of the North Sea fish assemblage: a biotic indicator of warming seas. J. Appl. Ecol. 45, 1029–1039 (2008).
Google Scholar 

3.
Sunday, J. M., Bates, A. E. & Dulvy, N. K. Thermal tolerance and the global redistribution of animals. Nat. Clim. Change 2, 686–690 (2012).
Google Scholar 

4.
Pinsky, M. L., Worm, B., Fogarty, M. J., Sarmiento, J. L. & Levin, S. A. Marine taxa track local climate velocities. Science 341, 1239–1242 (2013).
CAS  Google Scholar 

5.
Hutchins, L. W. The bases for temperature zonation in geographical distribution. Ecol. Monogr. 17, 325–335 (1947).
Google Scholar 

6.
Pineda, J., Reyns, N. B. & Starczak, V. R. Complexity and simplification in understanding recruitment in benthic populations. Pop. Ecol. 51, 17–32 (2009).
Google Scholar 

7.
Morgan, S. G., Shanks, A. L., MacMahan, J. H., Reniers, A. J. H. M. & Feddersen, F. Planktonic subsidies to surf-zone and intertidal communities. Annu. Rev. Mar. Sci. 10, 345–369 (2018).
Google Scholar 

8.
Gaylord, B. & Gaines, S. D. Temperature or transport? Range limits in marine species mediated solely by flow. Am. Nat. 155, 769–789 (2000).
Google Scholar 

9.
García Molinos, J., Burrows, M. T. & Poloczanska, E. S. Ocean currents modify the coupling between climate change and biogeographical shifts. Sci. Rep. 7, 1332 (2017).
Google Scholar 

10.
Kumagai, N. H. et al. Ocean currents and herbivory drive macroalgae-to-coral community shift under climate warming. Proc. Natl Acad. Sci. USA 115, 8990–8995 (2017).
Google Scholar 

11.
Harley, C. D. G. et al. The impacts of climate change in coastal marine systems. Ecol. Lett. 9, 228–241 (2006).
Google Scholar 

12.
Strathmann, M. F. Reproduction and Development of Marine Invertebrates of the Northern Pacific Coast (Univ. of Washington Press, 1987).

13.
Thorson, G. Reproductive and larval ecology of marine bottom invertebrates. Biol. Rev. 25, 1–45 (1950).
CAS  Google Scholar 

14.
Olive, P. J. W. Annual breeding cycles in marine invertebrates and environmental temperature: probing the proximate and ultimate causes of reproductive synchrony. J. Therm. Biol. 20, 79–90 (1995).
Google Scholar 

15.
Philippart, C. J. M. et al. Climate-related changes in recruitment of the bivalve Macoma balthica. Limnol. Oceanogr. 48, 2171–2185 (2003).
Google Scholar 

16.
Asch, R. G. Climate change and decadal shifts in the phenology of larval fishes in the California Current Ecosystem. Proc. Natl Acad. Sci. USA 112, E4065–E4074 (2015).
CAS  Google Scholar 

17.
Shearman, R. K. & Lentz, S. J. Long-term sea surface temperature variability along the U.S. East Coast. J. Phys. Oceanogr. 40, 1004–1017 (2010).
Google Scholar 

18.
Saba, V. S. et al. Enhanced warming of the northwest Atlantic Ocean under climate change. J. Geophys. Res. Oceans 121, 118–132 (2016).
Google Scholar 

19.
Castelao, R., Glenn, S. & Schofield, O. Temperature, salinity, and density variability in the central Middle Atlantic Bight. J. Geophys. Res. 115, C10005 (2010).
Google Scholar 

20.
Richaud, B., Kwon, Y.-O., Joyce, T. M., Fratantoni, P. S. & Lentz, S. J. Surface and bottom temperature and salinity climatology along the continental shelf off the Canadian and U.S. East Coasts. Cont. Shelf Res. 124, 165–181 (2016).
Google Scholar 

21.
Roughgarden, J., Gaines, S. & Possingham, H. Recruitment dynamics in complex life cycles. Science 241, 1460–1466 (1988).
CAS  Google Scholar 

22.
Connolly, S. R., Menge, B. A. & Roughgarden, J. A latitudinal gradient in recruitment of intertidal invertebrates in the northeast Pacific Ocean. Ecology 82, 1799–1813 (2001).
Google Scholar 

23.
Ma, H., Grassle, J. P. & Chant, R. J. Vertical distribution of bivalve larvae along a cross-shelf transect during summer upwelling and downwelling. Mar. Biol. 149, 1123–1138 (2006).
Google Scholar 

24.
Shanks, A. L. & Brink, L. Upwelling, downwelling, and cross-shelf transport of bivalve larvae: test of a hypothesis. Mar. Ecol. Prog. Ser. 302, 1–12 (2005).
Google Scholar 

25.
Drake, P. T., Edwards, C. A., Morgan, S. G. & Dever, E. P. Influence of larval behavior on transport and population connectivity in a realistic simulation of the California Current System. J. Mar. Res. 71, 317–350 (2013).
Google Scholar 

26.
Shanks, A. L. & Morgan, S. G. Testing the intermittent upwelling hypothesis: upwelling, downwelling, and subsidies to the intertidal zone. Ecol. Monogr. 88, 22–35 (2018).
Google Scholar 

27.
Menge, B. A. & Menge, D. N. L. Testing the intermittent upwelling hypothesis: comment. Ecology 100, e02476 (2019).
Google Scholar 

28.
Lentz, S. J. Seasonal variations in the circulation over the Middle Atlantic Bight continental shelf. J. Phys. Oceanogr. 38, 1486–1500 (2008).
Google Scholar 

29.
Gong, D., Kohut, J. T. & Glenn, S. M. Seasonal climatology of wind-driven circulation on the New Jersey Shelf. J. Geophys. Res. 115, C04006 (2010).
Google Scholar 

30.
Whitney, M. M. & Garvine, R. W. Wind influence on a coastal buoyant outflow. J. Geophys. Res. 110, C03014 (2005).
Google Scholar 

31.
Largier, J. L. Considerations in estimating larval dispersal distances from oceanographic data. Ecol. Appl. 13, S71–S89 (2003).
Google Scholar 

32.
Byers, J. E. & Pringle, J. M. Going against the flow: retention, range limits and invasions in advective environments. Mar. Ecol. Prog. Ser. 313, 27–41 (2006).
Google Scholar 

33.
Fuchs, H. L., Gerbi, G. P., Hunter, E. J. & Christman, A. J. Waves cue distinct behaviors and differentiate transport of congeneric snail larvae from sheltered versus wavy habitats. Proc. Natl Acad. Sci. USA 115, E7532–E7540 (2018).
CAS  Google Scholar 

34.
Sunday, J. M., Bates, A. E. & Dulvy, N. K. Global analysis of thermal tolerance and latitude in ectotherms. Proc. R. Soc. B 278, 1823–1830 (2011).
Google Scholar 

35.
Wilson, R. J. et al. Changes to the elevational limits and extent of species ranges associated with climate change. Ecol. Lett. 8, 1138–1146 (2005).
Google Scholar 

36.
Freeman, B. G., Scholer, M. N., Ruiz-Guttierrez, V. & Fitzpatrick, J. W. Climate change causes upslope shifts and mountaintop extirpations in a tropical bird community. Proc. Natl Acad. Sci. USA 115, 11982–11987 (2018).
CAS  Google Scholar 

37.
Free, C. M. et al. Impacts of historical warming on fisheries production. Science 363, 979–983 (2019).
CAS  Google Scholar 

38.
Young, I. R. & Ribal, A. Multiplatform evaluation of global trends in wind speed and wave height. Science 364, 548–552 (2019).
CAS  Google Scholar 

39.
Ocean Biogeographic Information System (Intergovernmental Oceanographic Commission of UNESCO, 2018); www.iobis.org

40.
Tingley, M. W. & Beissinger, S. R. Detecting range shifts from historical species occurrences: new perspectives on old data. Trends Ecol. Evol. 24, 625–633 (2009).
Google Scholar 

41.
Wigley, R. L. & Theroux, R. B. Atlantic Continental Shelf and Slope of the United States; Macrobenthic Invertebrate Fauna of the Middle Atlantic Bight Region; Faunal Composition and Quantitative Distribution Professional Paper No. 529-N (USGS, 1981).

42.
Amante, C. & Eakins, B. W. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis Technical Memorandum NESDIS NGDC-24 (National Geophysical Data Center, NOAA, 2009); https://doi.org/10.7289/V5C8276M

43.
Loarie, S. R. et al. The velocity of climate change. Nature 462, 1052–1057 (2009).
CAS  Google Scholar 

44.
Kang, D. & Curchitser, E. N. Gulf stream eddy characteristics in a high-resolution ocean model. J. Geophys. Res. Oceans 118, 4474–4487 (2013).
Google Scholar 

45.
Narváez, D. A. et al. Long-term dynamics in Atlantic surfclam (Spisula solidissima) populations: the role of bottom water temperature. J. Mar. Sys. 141, 136–148 (2015).
Google Scholar 

46.
Chen, Z., Curchitser, E., Chant, R. & Kang, D. Seasonal variability of the cold pool over the Mid-Atlantic Bight continental shelf. J. Geophys. Res. Oceans 123, 8203–8226 (2018).
Google Scholar 

47.
D’Errico, J. inpaint_nans (MATLAB Central File Exchange, 2019); https://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint_nans

48.
Gypaets trigradient2 (GitHub, 2020); https://www.github.com/Gypaets/trigradient2

49.
Yeager, S. & NCAR Staff The Climate Data Guide: COREv2 Air-Sea Surface Fluxes (UCAR, 2016); https://climatedataguide.ucar.edu/climate-data/corev2-air-sea-surface-fluxes

50.
National Water Information System Data (USGS, 2016); http://waterdata.usgs.gov/nwis/

51.
Lentz, S. J. Observations and a model of the mean circulation over the Middle Atlantic Bight continental shelf. J. Phys. Oceanogr. 38, 1203–1221 (2008).
Google Scholar 

52.
Shanks, A. L. Pelagic larval duration and dispersal distance revisited. Biol. Bull. 216, 373–385 (2009).
Google Scholar  More

1.
USDA. Coffee Annual Coffee. https://gain.fas.usda.gov/RecentGAINPublications/LOCK-UPREPORT_Pretoria_SouthAfrica-Republicof_10-29-2009.pdf (2019).
2.
Carvalho Guarçoni, R. et al. Influence of solar radiation and wet processing on the final quality of Arabica coffee. J. Food Qual. https://doi.org/10.1155/2018/6408571 (2018).
Article  Google Scholar 

3.
Iamanaka, B. T. et al. Reprint of ‘The mycobiota of coffee beans and its influence on the coffee beverage’. Food Res. Int. 61, 33–38. https://doi.org/10.1016/j.foodres.2014.05.023 (2014).
Article  Google Scholar 

4.
Barnes, E. C., Jumpathong, J., Lumyong, S., Voigt, K. & Hertweck, C. Daldionin, an unprecedented binaphthyl derivative, and diverse polyketide congeners from a fungal orchid endophyte. Chem. A Eur. J. 22, 4551–4555. https://doi.org/10.1002/chem.201504005 (2016).
Article  CAS  Google Scholar 

5.
Descroix, F. & Snoeck, J. Environmental factors suitable for coffee cultivation. In Coffee: Growing, Processing, Sustainable Production 164–177, https://doi.org/10.1002/9783527619627.ch6 (2008).

6.
De Bruyn, F. et al. Exploring the impacts of postharvest processing on the microbiota and metabolite profiles during green coffee bean production. Am. Soc. Microbiol. https://doi.org/10.1128/AEM.02398-16 (2016).
Article  Google Scholar 

7.
Hamdouche, Y. et al. Discrimination of post-harvest coffee processing methods by microbial ecology analyses. Food Control 65, 112–120. https://doi.org/10.1016/j.foodcont.2016.01.022 (2016).
Article  CAS  Google Scholar 

8.
Zhao, Q. et al. Long-term coffee monoculture alters soil chemical properties and microbial communities. Sci. Rep. 8, 1–11. https://doi.org/10.1038/s41598-018-24537-2 (2018).
ADS  Article  CAS  Google Scholar 

9.
Júnior, P. P. et al. Agroecological coffee management increases arbuscular mycorrhizal fungi diversity. PLoS ONE 14, 1–19. https://doi.org/10.1371/journal.pone.0209093 (2019).
Article  CAS  Google Scholar 

10.
Melloni, R. et al. Sistemas Agroflorestais cafeeiro-araucária e seu efeito na microbiota do solo e seus processos. Ciência Florest. 28, 784–795. https://doi.org/10.5902/1980509832392 (2018).
Article  Google Scholar 

11.
Oliveira, M. N. V. et al. Endophytic microbial diversity in coffee cherries of Coffea arabica from southeastern Brazil. Can. J. Microbiol. 59, 221–230. https://doi.org/10.1139/cjm-2012-0674 (2013).
Article  PubMed  CAS  Google Scholar 

12.
Nasanit, R. & Satayawut, K. Microbiological study during coffee fermentation of Coffea arabica var chiangmai 80 in Thailand. Kasetsart J. Nat. Sci. 49, 32–41 (2015).
CAS  Google Scholar 

13.
Evangelista, S. R. et al. Improvement of coffee beverage quality by using selected yeasts strains during the fermentation in dry process. Food Res. Int. 61, 183–195. https://doi.org/10.1016/j.foodres.2013.11.033 (2014).
Article  CAS  Google Scholar 

14.
Pereira, G. V. D. M. et al. Potential of lactic acid bacteria to improve the fermentation and quality of coffee during on-farm processing. Int. J. Food Sci. Technol. 51, 1689–1695. https://doi.org/10.1111/ijfs.13142 (2016).
Article  CAS  Google Scholar 

15.
Sahu, N., Duraisamy, V., Sahu, A., Lal, N. & K. Singh, S. Strength of microbes in nutrient cycling: A key to soil health. In Agriculturally Important Microbes for Sustainable Agriculture 69–86, https://doi.org/10.1007/978-981-10-5589-8_4 (2017).

16.
Zhang, S. J. et al. Following coffee production from cherries to cup: Microbiological and metabolomic analysis of wet processing of Coffea arabica. Appl. Environ. Microbiol. 85, 1–22. https://doi.org/10.1128/AEM.02635-18 (2019).
Article  CAS  Google Scholar 

17.
Ramos, C. L. et al. Determination of dynamic characteristics of microbiota in a fermented beverage produced by Brazilian Amerindians using culture-dependent and culture-independent methods. Int. J. Food Microbiol. 140, 225–231. https://doi.org/10.1016/j.ijfoodmicro.2010.03.029 (2010).
Article  PubMed  CAS  Google Scholar 

18.
Faoro, H. et al. Influence of soil characteristics on the diversity of bacteria in the Southern Brazilian Atlantic Forest. Appl. Environ. Microbiol. 76, 4744–4749. https://doi.org/10.1128/AEM.03025-09a (2010).
Article  PubMed  PubMed Central  CAS  Google Scholar 

19.
Defelipo, B. V. & Ribeiro, A. C. Análise química do solo (metodologia). Bol. Extensão 28, 1–26 (1997).
Google Scholar 

20.
Walters, W. et al. Improved bacterial 16S rRNA gene (V4 and V4–5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. Am. Soc. Microbiol. https://doi.org/10.1128/msystems.00009-15 (2015).
Article  Google Scholar 

21.
Pylro, V. S. et al. Data analysis for 16S microbial profiling from different benchtop sequencing platforms. J. Microbiol. Methods 107, 30–37. https://doi.org/10.1016/j.mimet.2014.08.018 (2014).
Article  PubMed  CAS  Google Scholar 

22.
Edgar, R. C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998. https://doi.org/10.1038/nmeth.2604 (2013).
Article  CAS  Google Scholar 

23.
Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. https://doi.org/10.1038/nmeth.f.303 (2010).
Article  PubMed  PubMed Central  CAS  Google Scholar 

24.
Edgar, R. C. UCHIME2: Improved chimera prediction for amplicon sequencing. BioRxiv https://doi.org/10.1101/074252 (2016).
Article  Google Scholar 

25.
Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. https://doi.org/10.1093/nar/gks1219 (2013).
Article  PubMed  PubMed Central  CAS  Google Scholar 

26.
Lozupone, C., Lladser, M. E., Knights, D., Stombaugh, J. & Knight, R. UniFrac: An effective distance metric for microbial community comparison. ISME J. 5, 169–172. https://doi.org/10.1038/ismej.2010.133 (2011).
Article  PubMed  Google Scholar 

27.
Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461. https://doi.org/10.1093/bioinformatics/btq461 (2010).
Article  PubMed  PubMed Central  CAS  Google Scholar 

28.
Bengtsson-Palme, J. et al. Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data. Methods Ecol. Evol. 4, 914–919. https://doi.org/10.1111/2041-210X.12073 (2013).
Article  Google Scholar 

29.
Nilsson, R. H. et al. The UNITE database for molecular identification of fungi: Handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 47, D259–D264. https://doi.org/10.1093/nar/gky1022 (2019).
Article  PubMed  CAS  Google Scholar 

30.
Oksanen, J. et al. Community Ecology Package. 1–296, https://cran.r-project.org/web/packages/vegan/vegan.pdf (2019).

31.
R Core Team. R: A Language and Environment for Statistical Computing. https://www.r-project.org/ (2018).

32.
Borcard, D. et al. Canonical ordination. In Numerical Ecology with R 153–225, https://doi.org/10.1007/978-1-4419-7976-6_6 (2011).

33.
Gomes, D. G. E. et al. Bats perceptually weight prey cues across sensory systems when hunting in noise. Science 353, 1277–1280. https://doi.org/10.1126/science.aaf7934 (2016).
ADS  Article  PubMed  CAS  Google Scholar 

34.
Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey data. PLOS Comput. Biol. 8, 1–11. https://doi.org/10.1371/journal.pcbi.1002687 (2012).
Article  CAS  Google Scholar 

35.
Watts, S. C., Ritchie, S. C., Inouye, M. & Holt, K. E. FastSpar: Rapid and scalable correlation estimation for compositional data. Bioinformatics 35, 1064–1066. https://doi.org/10.1093/bioinformatics/bty734 (2019).
Article  PubMed  CAS  Google Scholar 

36.
Csardi, G. & Nepusz, T. The igraph software package for complex network research. InterJournal Complex Syst 1695, 1–9 (2006).
Google Scholar 

37.
Avelino, J. et al. Effects of slope exposure, altitude and yield on coffee quality in two altitude terroirs of Costa Rica, Orosi and Santa María de Dota. J. Sci. Food Agric. 85, 1869–1876. https://doi.org/10.1002/jsfa.2188 (2005).
Article  CAS  Google Scholar 

38.
Wei, L., Wai, M., Curran, P., Yu, B. & Quan, S. Coffee fermentation and flavor—An intricate and delicate relationship. Food Chem. 185, 182–191. https://doi.org/10.1016/j.foodchem.2015.03.124 (2015).
Article  CAS  Google Scholar 

39.
Fulthorpe, R., Martin, A. R. & Isaac, M. E. Root endophytes of coffee ( Coffea arabica): Variation across climatic gradients and relationships with functional traits. Phytobiomes J. 4, 27–39. https://doi.org/10.1094/PBIOMES-04-19-0021-R (2020).
Article  Google Scholar 

40.
Chu, H. et al. Effects of slope aspects on soil bacterial and arbuscular fungal communities in a boreal forest in China. Pedosphere 26, 226–234. https://doi.org/10.1016/S1002-0160(15)60037-6 (2016).
Article  Google Scholar 

41.
Karungi, J. et al. Elevation and cropping system as drivers of microclimate and abundance of soil macrofauna in coffee farmlands in mountainous ecologies. Appl. Soil Ecol. 132, 126–134. https://doi.org/10.1016/J.APSOIL.2018.08.003 (2018).
Article  Google Scholar 

42.
Ferreira, W. P. M., Queiroz, D. M., Silvac, S. A., Tomaz, R. S. & Corrêa, P. C. Effects of the orientation of the mountainside, altitude and varieties on the quality of the coffee beverage from the “Matas de Minas” region, Brazilian Southeast. Am. J. Plant Sci. 7, 1291–1303. https://doi.org/10.4236/ajps.2016.78124 (2016).
Article  Google Scholar 

43.
Velmourougane, K. Impact of organic and conventional systems of coffee farming on soil properties and culturable microbial diversity. Scientifica 1–9, 2016. https://doi.org/10.1155/2016/3604026 (2016).
Article  CAS  Google Scholar 

44.
Siles, J. A. & Margesin, R. Abundance and diversity of bacterial, archaeal, and fungal communities along an altitudinal gradient in alpine forest soils: What are the driving factors?. Soil Microbiol. 72, 207–220. https://doi.org/10.1007/s00248-016-0748-2 (2016).
Article  Google Scholar 

45.
Frank, A., Saldierna Guzmán, J. & Shay, J. Transmission of bacterial endophytes. Microorganisms 5, 70. https://doi.org/10.3390/microorganisms5040070 (2017).
Article  PubMed Central  CAS  Google Scholar 

46.
Haile, M. & Kang, W. H. The role of microbes in coffee fermentation and their impact on coffee quality. J. Food Qual. 2019, 6. https://doi.org/10.1155/2019/4836709 (2019).
Article  CAS  Google Scholar 

47.
Decazy, F. et al. Quality of different Honduran coffees in relation to several environments. J. Food Sci. 68, 2356–2361. https://doi.org/10.1111/j.1365-2621.2003.tb05772.x (2003).
Article  CAS  Google Scholar 

48.
de Melo Pereira, G. V. et al. Conducting starter culture-controlled fermentations of coffee beans during on-farm wet processing: Growth, metabolic analyses and sensorial effects. Food Res. Int. 75, 348–356. https://doi.org/10.1016/j.foodres.2015.06.027 (2015).
Article  PubMed  CAS  Google Scholar 

49.
Zhang, W. et al. Microbial diversity in two traditional bacterial douchi from Gansu province in northwest China using Illumina sequencing. PLoS ONE 13, 1–16. https://doi.org/10.1371/journal.pone.0194876 (2018).
Article  CAS  Google Scholar 

50.
Tolessa, K., D’heer, J., Duchateau, L. & Boeckx, P. Influence of growing altitude, shade and harvest period on quality and biochemical composition of Ethiopian specialty coffee. J. Sci. Food Agric. 97, 2849–2857. https://doi.org/10.1002/jsfa.8114 (2017).
Article  PubMed  CAS  Google Scholar 

51.
Batista, D. et al. Legitimacy and implications of reducing Colletotrichum kahawae to subspecies in plant pathology. Front. Plant Sci. 7, 1–9. https://doi.org/10.3389/fpls.2016.02051 (2017).
Article  CAS  Google Scholar 

52.
Wei, Z. et al. Trophic network architecture of root-associated bacterial communities determines pathogen invasion and plant health. Nat. Commun. 6, 1–9. https://doi.org/10.1038/ncomms9413 (2015).
ADS  Article  CAS  Google Scholar  More