Philippa Kaur

ParticipantsThe Ministry of Health, Labour, and Welfare of Japan has been conducting The Longitudinal Survey of Newborns in the 21st Century since 2001 to establish strategies to counter the declining birthrate in Japan. The survey targeted all babies born in Japan between January 10 and 17 or between July 10 and 17 of 2001. Baseline questionnaires were sent to a total of 53,575 families when eligible babies reached the age of 6 months and 47,015 families initially completed the baseline questionnaire (88% response rate). These respondents were mailed follow-up questionnaires to investigate medical conditions and behaviors when children reached the ages of 1.5, 2.5, 3.5, 4.5, 5.5, 7, 8, 9, 10, 11, 12, 13, 14, and 15 years20,21,22,23. Birth record data from Vital Statistics of Japan are also linked for each child participating in the study. The current study included data for children/families who responded both to the baseline questionnaire and the fifteenth questionnaire at age 15 years.The baseline survey at age 6 months included questions regarding children’s perinatal status as well as household and socioeconomic factors such as parental academic attainment, parental smoking status, and daycare attendance. The subsequent annual surveys starting at age 1.5 years included questions regarding each child’s height, weight and health status. We excluded 2382 children born before 37 weeks of pregnancy and one child with responses only for the baseline survey and the survey at age 15 years. A total of 26,778 children (315,581 data points) were included in the final analysis. A total of 11,141 children (41.61%) had responses to all 15 questionnaires between the ages of 6 months and 15 years, and responses to more than 12 questionnaires were available for the majority (91.94%) of children (Fig. 1, Table S1).Figure 1Flowchart of study participants.Full size imageMeasuresWe calculated BMI based on each participant’s reported annual height and weight. Each participant’s annual BMI was converted to a BMI Z-score using smoothed L, M, and S values for BMI standards from a representative population of Japanese children24. Briefly, the LMS (lambda–mu–sigma) method is a method proposed by Cole et al. to monitor changes in the skewness of the distribution during childhood as a way of constructing normalized growth standards25. Participants were then classified into four BMI categories based on the World Health Organization (WHO) criteria26: underweight (BMI standard deviation [SD] score of − 5 or more but less than − 2), normal weight (BMI SD score of − 2 or more but less than 1), overweight (BMI SD score of 1 or more but less than 2), and obese (BMI SD score of 2 or more but less than 5). The definitions of overweight and obesity were different for children under 5 years of age: a BMI Z-score of 2 SD or more was categorized as overweight and a BMI Z-score of 3 SD or more was categorized as obese. BMI category at age 15 years was the main outcome of interest in the current study.We also calculated annual height growth for each participant by subtracting the height reported at the previous survey from that reported in the current survey. For annual height growth between 5.5 and 7 years of age, this value was multiplied by 2/3 because of the 1.5-year interval between surveys.Statistical analysesWe first compared baseline characteristics among the four BMI categories (underweight, normal weight, overweight and obese) at age 15 years. To evaluate potential selection bias resulting from losses to follow-up, we also compared the baseline characteristics of children included in the analysis and those of children lost to follow-up through to the fifteenth survey (at age 15 years).We retrospectively examined annual aggregate categorical changes in individuals of the four BMI categories (groups) at age 15 years. For each group, the proportion of each BMI category at each survey between the ages of 1.5 and 14 years was calculated. In addition, we prospectively calculated the proportion of children in each BMI category at each survey between the ages of 1.5 and 14 years who eventually became underweight, normal weight, overweight, or obese at age 15 years. Note that these analyses were based on aggregate data and do not describe individual BMI changes and were performed using only the data obtained without imputation of missing values.Under the assumption that missing data were missing at random, mixed effect models with natural cubic regression splines were applied to calculate the trajectories of BMI Z-scores and annual BMI Z-score changes through age 15 years for participants of each BMI category at age 15 years. Knots at seven locations were placed in percentiles of age to yield a sufficient number of measurements between each consecutive knot (age 1.5, 3.5, 5.5, 8.5, 11, 13 and 15 years), as recommended by Harrell27. The mixed effect model is useful for describing population average growth trajectories and individual growth trajectories even when data are not available for all children at all ages28,29,30,31. Briefly, the population average growth trajectory was modeled with fixed effects, while the individual variability is represented as random effects.After fitting individual BMI trajectories using a mixed-effects model with natural cubic spline function, we estimated individual adiposity rebound timing as the age where the first derivative of the trajectory reached its minimum and the second derivative was positive32. Children were then classified into five categories (1.5–2.5 years, 3.5–4.5 years, 5.5–7 years, 8–10 years, and 11 years or older) for analysis of adiposity rebound timing33,34. The distribution of adiposity rebound timing was calculated for individuals of each BMI status at age 15 years overall and by gender.Finally, we modelled annual height change and its associations with BMI status at age 15 years separately for each gender using mixed-effects models with natural cubic regression splines.All statistical analyses were performed using Stata version 16 (StataCorp LLC, College Station, TX, USA). This study was approved by the Institutional Review Board at Okayama University Graduate School of Medicine, Dentistry, and Pharmaceutical Sciences (No.1506-073) and was conducted in accordance with the 1964 Helsinki Declaration and Ethical Guidelines for Medical and Health Research Involving Human Subjects. Informed consent was obtained by the opt-out method on the university’s website. More

1.Rivadeneira, M. M. et al. Testing the abundant-centre hypothesis using intertidal porcelain crabs along the Chilean coast: Linking abundance and life-history variation. J. Biogeogr. 37, 486–498 (2010).
Google Scholar 
2.Hutchins, L. W. The bases for temperature zonation in geographical distribution. Ecol. Monogr. 17, 325–335 (1947).
Google Scholar 
3.Lewis, J. R. Latitudinal trends in reproduction, recruitment and population characteristics of some rocky littoral molluscs and cirripedes. Hydrobiologia 142, 1–13 (1986).
Google Scholar 
4.Bernardo, J. The particular maternal effect of propagule size, especially egg size: Patterns, models, quality of evidence and interpretations. Am. Zool. 36, 216–236 (1996).
Google Scholar 
5.Thorson, G. Reproductive and larval ecology of marine bottom invertebrates. Biol. Rev. 25, 1–45 (1950).CAS 
PubMed 

Google Scholar 
6.Marshall, D. J., Pettersen, A. K. & Cameron, H. A global synthesis of offspring size variation, its eco-evolutionary causes and consequences. Funct. Ecol. 32, 1436–1446 (2018).
Google Scholar 
7.Des Roches, S. et al. The ecological importance of intraspecific variation. Nat. Ecol. Evol. 2, 57–64 (2018).PubMed 

Google Scholar 
8.Violle, C. et al. Let the concept of trait be functional!. Oikos 116, 882–892 (2007).
Google Scholar 
9.Sides, C. B. et al. Revisiting Darwin’s hypothesis: Does greater intraspecific variability increase species’ ecological breadth?. Am. J. Bot. 101, 56–62 (2014).PubMed 

Google Scholar 
10.Moran, E. V., Hartig, F. & Bell, D. M. Intraspecific trait variation across scales: Implications for understanding global change responses. Glob. Change Biol. 22, 137–150 (2016).ADS 

Google Scholar 
11.Violle, C. et al. The return of the variance: Intraspecific variability in community ecology. Trends Ecol. Evol. 27, 244–252 (2012).PubMed 

Google Scholar 
12.Stark, J., Lehman, R., Crawford, L., Enquist, B. J. & Blonder, B. Does environmental heterogeneity drive functional trait variation? A test in montane and alpine meadows. Oikos 126, 1650–1659 (2017).
Google Scholar 
13.Stearns, S. C. The Evolution of Life Histories. xii, 249p. No. 575 S81 (Oxford, Oxford University, 1992).14.Vance, R. R. On reproductive strategies in marine benthic invertebrates. Am. Nat. 107, 339–352 (1973).
Google Scholar 
15.Levitan, D. R. Gamete traits influence the variance in reproductive success, the intensity of sexual selection, and the outcome of sexual conflict among congeneric sea urchins. Evolution 62, 1305–1316 (2008).PubMed 

Google Scholar 
16.Lavorel, S. & Garnier, E. Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail. Funct. Ecol. 16, 545–556 (2002).
Google Scholar 
17.Pineda, M. C. et al. Tough adults, frail babies: An analysis of stress sensitivity across early life-history stages of widely introduced marine invertebrates. PLoS ONE 7, e46672 (2012).ADS 
CAS 
PubMed 
PubMed Central 

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

Google Scholar 
19.Stein, A., Gerstner, K. & Kreft, H. Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecol. Lett. 17, 866–880 (2014).PubMed 

Google Scholar 
20.Foo, S. A. & Byrne, M. Marine gametes in a changing ocean: Impacts of climate change stressors on fecundity and the egg. Mar. Environ. Res. 128, 12–24 (2017).CAS 
PubMed 

Google Scholar 
21.Dahlhoff, E. P. Biochemical indicators of stress and metabolism: Applications for marine ecological studies. Annu. Rev. Physiol. 66, 183–207 (2004).CAS 
PubMed 

Google Scholar 
22.Soudant, P. et al. Comparison of the lipid class and fatty acid composition between a reproductive cycle in nature and a standard hatchery conditioning of the Pacific Oyster Crassostrea gigas. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 123, 209–222 (1999).
Google Scholar 
23.Lester, S. E., Gaines, S. D. & Kinlan, B. P. Reproduction on the edge: Large-scale patterns of individual performance in a marine invertebrate. Ecology 88, 2229–2239 (2007).PubMed 

Google Scholar 
24.Helmuth, B., Mieszkowska, N., Moore, P. & Hawkins, S. J. Living on the edge of two changing worlds: Forecasting the responses of rocky intertidal ecosystems to climate change. Annu. Rev. Ecol. Evol. Syst. 37, 373–404 (2006).
Google Scholar 
25.Sagarin, R. D., Barry, J. P., Gilman, S. E. & Baxter, C. H. Climate-related change in an intertidal community over short and long time scales. Ecol. Monogr. 69, 465–490 (1999).
Google Scholar 
26.Dubois, S., Retière, C. & Olivier, F. Biodiversity associated with Sabellaria alveolata (Polychaeta: Sabellariidae) reefs: Effects of human disturbances. J. Mar. Biol. Assoc. UK 82, 817–826 (2002).
Google Scholar 
27.Jones, A. G., Dubois, S. F., Desroy, N. & Fournier, J. Interplay between abiotic factors and species assemblages mediated by the ecosystem engineer Sabellaria alveolata (Annelida: Polychaeta). Estuar. Coast. Shelf Sci. 200, 1–18 (2018).ADS 

Google Scholar 
28.Bonifazi, A. et al. Macrofaunal biodiversity associated with different developmental phases of a threatened Mediterranean Sabellaria alveolata (Linnaeus, 1767) reef. Mar. Environ. Res. 145, 97–111 (2019).CAS 
PubMed 

Google Scholar 
29.Holt, T. J., Biogenic Reefs. An Overview of Dynamic and Sensitivity Characteristics for Conservation Management of Marine SACs. UK Marine SACs Project (1998).30.Crisp, D. The effects of the severe winter of 1962–1963 on marine life in Britain. J. Anim. Ecol. 33, 165–210 (1964).
Google Scholar 
31.Firth, L. B. et al. Historical comparisons reveal multiple drivers of decadal change of an ecosystem engineer at the range edge. Ecol. Evol. 5, 3210–3222 (2015).PubMed 
PubMed Central 

Google Scholar 
32.Firth, L. B. et al. Specific niche requirements underpin multidecadal range edge stability, but may introduce barriers for climate change adaptation. Divers. Distrib. 27, 668–683 (2021).
Google Scholar 
33.Wethey, D. S. et al. Response of intertidal populations to climate: Effects of extreme events versus long term change. J. Exp. Mar. Biol. Ecol. 400, 132–144 (2011).
Google Scholar 
34.Sagarin, R. D. & Gaines, S. D. Geographical abundance distributions of coastal invertebrates: Using one-dimensional ranges to test biogeographic hypotheses. J. Biogeogr. 29, 985–997 (2002).
Google Scholar 
35.Rahman, M. A., Henderson, S., Miller-Ezzy, P., Li, X. X. & Qin, J. G. Immune response to temperature stress in three bivalve species: Pacific oyster Crassostrea gigas, Mediterranean mussel Mytilus galloprovincialis and mud cockle Katelysia rhytiphora. Fish Shellfish Immunol. 86, 868–874 (2019).CAS 
PubMed 

Google Scholar 
36.Osada, M., Nishikawa, M. & Nomura, T. Involvement of prostaglandins in the spawning of the scallop, Patinopecten yessoensis. Comp. Biochem. Physiol. C 94, 595–601 (1989).
Google Scholar 
37.Stanley, D. W. & Howard, R. W. The biology of prostaglandins and related eicosanoids in invertebrates: Cellular, organismal and ecological actions. Am. Zool. 38, 369–381 (1998).CAS 

Google Scholar 
38.Pernet, F., Tremblay, R., Comeau, L. & Guderley, H. Temperature adaptation in two bivalve species from different thermal habitats: Energetics and remodelling of membrane lipids. J. Exp. Biol. 210, 2999–3014 (2007).PubMed 

Google Scholar 
39.Muir, A. P., Nunes, F. L. D., Dubois, S. F. & Pernet, F. Lipid remodelling in the reef-building honeycomb worm, Sabellaria alveolata, reflects acclimation and local adaptation to temperature. Sci. Rep. 6, 35669 (2016).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
40.Hulbert, A. & Else, P. L. Membranes as possible pacemakers of metabolism. J. Theor. Biol. 199, 257–274 (1999).ADS 
CAS 
PubMed 

Google Scholar 
41.Brokordt, K. B., Himmelman, J. H., Nusetti, O. A. & Guderley, H. E. Reproductive investment reduces recuperation from exhaustive escape activity in the tropical scallop Euvola zizac. Mar. Biol. 137, 857–865 (2000).CAS 

Google Scholar 
42.Levitan, D. R. & Roitberg, B. D. Optimal egg size in marine invertebrates: Theory and phylogenetic analysis of the critical relationship between egg size and development time in echinoids. Am. Nat. 156, 175–192 (2000).PubMed 

Google Scholar 
43.Moran, A. L. & McAlister, J. S. Egg size as a life history character of marine invertebrates: Is it all it’s cracked up to be?. Biol. Bull. 216, 226–242 (2009).PubMed 

Google Scholar 
44.Marshall, D. J. & Burgess, S. C. Deconstructing environmental predictability: Seasonality, environmental colour and the biogeography of marine life histories. Ecol. Lett. 18, 174–181 (2015).PubMed 

Google Scholar 
45.Racault, M.-F., Le Quéré, C., Buitenhuis, E., Sathyendranath, S. & Platt, T. Phytoplankton phenology in the global ocean. Ecol. Indic. 14, 152–163 (2012).
Google Scholar 
46.Henson, S., Cole, H., Beaulieu, C. & Yool, A. The impact of global warming on seasonality of ocean primary production. Biogeosciences 10, 4357–4369 (2013).ADS 

Google Scholar 
47.Morim, J. et al. Robustness and uncertainties in global multivariate wind-wave climate projections. Nat. Clim. Change 9, 711–718 (2019).ADS 

Google Scholar 
48.Stillman, J. H. Heat waves, the new normal: Summertime temperature extremes will impact animals, ecosystems, and human communities. Physiology 34, 86–100 (2019).CAS 
PubMed 

Google Scholar 
49.McCarthy, D., Young, C. & Emson, R. Influence of wave-induced disturbance on seasonal spawning patterns in the sabellariid polychaete Phragmatopoma lapidosa. Mar. Ecol. Prog. Ser. 256, 123–133 (2003).ADS 

Google Scholar 
50.Aviz, D., Pinto, A. J. A., Ferreira, M. A. P., Rocha, R. M. & Rosa Filho, J. S. Reproductive biology of Sabellaria wilsoni (Sabellariidae: Polychaeta), an important ecosystem engineer on the Amazon coast. J. Mar. Biol. Assoc. UK https://doi.org/10.1017/S0025315416001776 (2016).Article 

Google Scholar 
51.Bowman, R. S. & Lewis, J. Annual fluctuations in the recruitment of Patella vulgata L. J. Mar. Biol. Assoc. U. K. 57, 793–815 (1977).
Google Scholar 
52.Sagarin, R. D. & Somero, G. N. Complex patterns of expression of heat-shock protein 70 across the southern biogeographical ranges of the intertidal mussel Mytilus californianus and snail Nucella ostrina. J. Biogeogr. 33, 622–630 (2006).
Google Scholar 
53.Wernberg, T. et al. An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat. Clim. Change 3, 78–82 (2013).ADS 

Google Scholar 
54.Firth, L. B., Knights, A. M. & Bell, S. S. Air temperature and winter mortality: Implications for the persistence of the invasive mussel, Perna viridis in the intertidal zone of the south-eastern United States. J. Exp. Mar. Biol. Ecol. 400, 250–256 (2011).
Google Scholar 
55.Seabra, R., Wethey, D. S., Santos, A. M. & Lima, F. P. Side matters: Microhabitat influence on intertidal heat stress over a large geographical scale. J. Exp. Mar. Biol. Ecol. 400, 200–208 (2011).
Google Scholar 
56.Meneghesso, C. et al. Remotely-sensed L4 SST underestimates the thermal fingerprint of coastal upwelling. Remote Sens. Environ. 237, 111588 (2020).ADS 

Google Scholar 
57.Marshall, D. J. & Keough, M. J. The evolutionary ecology of offspring size in marine invertebrates. in Advances in Marine Biology, 1–60. https://doi.org/10.1016/S0065-2881(07)53001-4 (Elsevier, 2007).58.Albert, C. H. et al. A multi-trait approach reveals the structure and the relative importance of intra- vs. interspecific variability in plant traits: Intra- vs. interspecific variability in plant traits. Funct. Ecol. 24, 1192–1201 (2010).
Google Scholar 
59.Olofsson, H., Ripa, J. & Jonzén, N. Bet-hedging as an evolutionary game: The trade-off between egg size and number. Proc. R. Soc. B Biol. Sci. 276, 2963–2969 (2009).
Google Scholar 
60.Osovitz, C. J. & Hofmann, G. E. Marine macrophysiology: Studying physiological variation across large spatial scales in marine systems. Comp. Biochem. Physiol. A. Mol. Integr. Physiol. 147, 821–827 (2007).PubMed 

Google Scholar 
61.Clarke, A. Reproduction in the cold: Thorson revisited. Invertebr. Reprod. Dev. 22, 175–183 (1992).
Google Scholar 
62.Hawkins, S. J. et al. Distinguishing globally-driven changes from regional- and local-scale impacts: The case for long-term and broad-scale studies of recovery from pollution. Mar. Pollut. Bull. 124, 573–586 (2017).CAS 
PubMed 

Google Scholar 
63.Dahlhoff, E. P., Stillman, J. H. & Menge, B. A. Physiological community ecology: Variation in metabolic activity of ecologically important rocky intertidal invertebrates along environmental gradients. Integr. Comp. Biol. 42, 862–871 (2002).PubMed 

Google Scholar 
64.Nunes, F. L. D., Rigal, F., Dubois, S. F. & Viard, F. Looking for diversity in all the right places? Genetic diversity is highest in peripheral populations of the reef-building polychaete Sabellaria alveolata. Mar. Biol. 168, 63 (2021).
Google Scholar 
65.Bush, L. E. Stability and Variability of the Ecosystem Engineer Sabellaria alveolata on Differing Temporal and Spatial Scales (Bangor University, 2016).
Google Scholar 
66.Lourenço, C. R., Nicastro, K. R., McQuaid, C. D., Krug, L. A. & Zardi, G. I. Strong upwelling conditions drive differences in species abundance and community composition along the Atlantic coasts of Morocco and Western Sahara. Mar. Biodivers. 50, 15 (2020).
Google Scholar 
67.Ritchie, H. & Marshall, D. J. Fertilisation is not a new beginning: Sperm environment affects offspring developmental success. J. Exp. Biol. 216, 3104–3109 (2013).PubMed 

Google Scholar 
68.Dubois, S., Comtet, T., Retière, C. & Thiébaut, E. Distribution and retention of Sabellaria alveolata larvae (Polychaeta: Sabellariidae) in the Bay of Mont-Saint-Michel, France. Mar. Ecol. Prog. Ser. 346, 243–254 (2007).ADS 
CAS 

Google Scholar 
69.Costello, D. P., Henley, C., & Marine Biological Laboratory (Woods Hole, Mass.). Methods for obtaining and handling marine eggs and embryos [by] Donald P. Costello and Catherine Henley. ([s.n.], 1971). https://doi.org/10.5962/bhl.title.1020.70.Gruet, Y. Aspects morphologiques et dynamiques de constructions de l’Annélide polychete Sabellaria alveolata (Linne). Rev. Trav. Inst. Pêch. Marit. 36, 131–161 (1972).
Google Scholar 
71.Saulquin, B., Gohin, F. & Garrello, R. Regional objective analysis for merging high-resolution MERIS, MODIS/Aqua, and SeaWiFS Chlorophyll-a data from 1998 to 2008 on the European Atlantic Shelf. IEEE Trans. Geosci. Remote Sens. 49, 143–154 (2011).ADS 

Google Scholar 
72.Gohin, F. Annual cycles of chlorophyll-a, non-algal suspended particulate matter, and turbidity observed from space and in-situ in coastal waters. Ocean Sci. 7, 705–732 (2011).ADS 
CAS 

Google Scholar 
73.Seabra, R., Wethey, D. S., Santos, A. M. & Lima, F. P. Understanding complex biogeographic responses to climate change. Sci. Rep. 5, 12930 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
74.Schlegel, R. W., Darmaraki, S., Benthuysen, J. A., Filbee-Dexter, K. & Oliver, E. C. J. Marine cold-spells. Progress Oceanogr. 198, 102684 (2021).
Google Scholar 
75.Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. Prog. Oceanogr. 141, 227–238 (2016).ADS 

Google Scholar 
76.Schlegel, R. W., Oliver, E. C., Hobday, A. J. & Smit, A. J. Detecting marine heatwaves with sub-optimal data. Front. Mar. Sci. 6, 737 (2019).
Google Scholar 
77.Egbert, G. D., Erofeeva, S. Y. & Ray, R. D. Assimilation of altimetry data for nonlinear shallow-water tides: Quarter-diurnal tides of the Northwest European Shelf. Cont. Shelf Res. 30, 668–679 (2010).ADS 

Google Scholar 
78.Burrows, M., Harvey, R. & Robb, L. Wave exposure indices from digital coastlines and the prediction of rocky shore community structure. Mar. Ecol. Prog. Ser. 353, 1–12 (2008).ADS 

Google Scholar 
79.Wessel, P. & Smith, W. H. F. A global, self-consistent, hierarchical, high-resolution shoreline database. J. Geophys. Res. Solid Earth 101, 8741–8743 (1996).
Google Scholar 
80.Seers, B. fetchR: Calculate Wind Fetch. R Package Version 2-1 (2017).81.Guillaume, A. S., Monro, K. & Marshall, D. J. Transgenerational plasticity and environmental stress: Do paternal effects act as a conduit or a buffer?. Funct. Ecol. 30, 1175–1184 (2016).
Google Scholar 
82.Curd, A. et al. Connecting organic to mineral: How the physiological state of an ecosystem-engineer is linked to its habitat structure. Ecol. Indic. 98, 49–60 (2019).CAS 

Google Scholar 
83.Gruet, Y. & Lassus, P. Contribution a l’etude de la biologie reproductive d’une population naturelle de l’Annelide Polychete, Sabellaria alveolata (Linnaeus). Ann. Inst. Oceanogr. Monaco 59, 127–140 (1983).
Google Scholar 
84.Hazel, J. The role of alterations in membrane lipid composition in enabling physiological adaptation of organisms to their physical environment. Prog. Lipid Res. 29, 167–227 (1990).CAS 
PubMed 

Google Scholar 
85.Hochachka, P. W. & Somero, G. N. Biochemical Adaptation: Mechanism and Process in Physiological Evolution (Oxford University Press, 2002).
Google Scholar 
86.Abele, D. & Puntarulo, S. Formation of reactive species and induction of antioxidant defence systems in polar and temperate marine invertebrates and fish. Comp. Biochem. Physiol. A Mol. Integr. Physiol. 138, 405–415 (2004).PubMed 

Google Scholar 
87.Folch, J., Lees, M. & Stanley, G. H. S. A simple method for the isolation and purification of total lipides from animal tissues. J. Biol. Chem. 226, 497–509. http://www.jbc.org/content/226/1/497 (1957).CAS 
PubMed 

Google Scholar 
88.Sieracki, C., Sieracki, M. & Yentsch, C. An imaging-in-flow system for automated analysis of marine microplankton. Mar. Ecol. Prog. Ser. 168, 285–296 (1998).ADS 

Google Scholar 
89.Pasteels, J. J. Etude au microscope électronique de la réaction corticale. II. La réaction corticale de l’oeuf vierge de Sabellaria alveolata. J. Embryol. Exp. Morphol. 13, 327–339 (1965).CAS 
PubMed 

Google Scholar 
90.Doledec, S. & Chessel, D. Co-inertia analysis: An alternative method for studying species-environment relationships. Freshw. Biol. 31, 277–294 (1994).
Google Scholar 
91.Robert, P. & Escoufier, Y. A unifying tool for linear multivariate statistical methods: The RV-coefficient. J. R. Stat. Soc. Ser. C Appl. Stat. 25, 257–265 (1976).MathSciNet 

Google Scholar 
92.Legendre, P. & Legendre, L. Ecological resemblance. in Developments in Environmental Modelling Chapter 7, Vol. 24, 265–335 (Elsevier, 2012).93.Borcard, D., Legendre, P. & Drapeau, P. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055 (1992).
Google Scholar 
94.Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D. Variation partitioning of species data matrices: Estimation and comparison of fractions. Ecology 87, 2614–2625 (2006).PubMed 

Google Scholar 
95.Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).
Google Scholar 
96.Messier, J., McGill, B. J. & Lechowicz, M. J. How do traits vary across ecological scales? A case for trait-based ecology: How do traits vary across ecological scales?. Ecol. Lett. 13, 838–848 (2010).PubMed 

Google Scholar 
97.Rao, C. R. The use and interpretation of principal component analysis in applied research. Sankhyā Indian J. Stat. Ser. A (1961-2002) 26, 329–358 (1964).98.R Core Team: A language and environment for statistical computing. Available from: https://www.R-project.org/ (2018).99.Oksanen, J. et al. Package ‘vegan’. Commun. Ecol. Package Version 2, 1–295 (2013).
Google Scholar  More

1.Horning, N., Robinson, J. A., Sterling, E. J., Turner, W. & Spector, S. Remote sensing for ecology and conservation. Techniques in Ecology & Conservation Series (Oxford University Press, 2010).2.Li, J. et al. A review of remote sensing for environmental monitoring in China. Remote Sens. 12, 1130 (2020).ADS 

Google Scholar 
3.Carter, W. D. & Paulson, R. W. Introduction to monitoring dynamic environmental phenomena of the world using satellite data collection systems. (U.S. Geological Survey, 1979).4.Nurdin, S., Mustapha, M. A. & Lihan, T. The relationship between sea surface temperature and chlorophyll-a concentration in fisheries aggregation area in the archipelagic waters of spermonde using satellite images. AIP Conf. Proc. 1571, 466–472 (2013).ADS 

Google Scholar 
5.Ward, D., Phinn, S. R. & Murray, A. T. Monitoring growth in rapidly urbanizing areas using remotely sensed data. Prof. Geogr. 52, 371–386 (2000).
Google Scholar 
6.Singh, A. Review article: Digital change detection techniques using remotely-sensed data. Int. J. Remote Sens. 10, 989–1003 (1989).
Google Scholar 
7.Dewan, A. M. & Yamaguchi, Y. Land use and land cover change in Greater Dhaka, Bangladesh: Using remote sensing to promote sustainable urbanization. Appl. Geogr. 29, 390–401 (2009).
Google Scholar 
8.Green, K., Kempka, D. & Lackey, L. Using remote sensing to detect and monitor land-cover and land-use change. Photogramm. Eng. Remote Sens. 60, 331–337 (1994).
Google Scholar 
9.Nagendra, H. Using remote sensing to assess biodiversity. Int. J. Remote Sens. 22, 2377–2400 (2001).
Google Scholar 
10.Rosenqvist, Å., Milne, A., Lucas, R., Imhoff, M. & Dobson, C. A review of remote sensing technology in support of the Kyoto Protocol. Environ. Sci. Policy 6, 441–455 (2003).
Google Scholar 
11.Liu, J. A process-based boreal ecosystem productivity simulator using remote sensing inputs. Remote Sens. Environ. 62, 158–175 (1997).ADS 

Google Scholar 
12.Colwell, R. R. Global climate and infectious disease: The cholera paradigm. Science 274, 2025–2031 (1996).ADS 
PubMed 
CAS 

Google Scholar 
13.Escobar, L. E. et al. A global map of suitability for coastal Vibrio cholerae under current and future climate conditions. Acta Trop. 149, 202–211 (2015).PubMed 

Google Scholar 
14.Watts, N. et al. The 2019 report of The Lancet Countdown on health and climate change: Ensuring that the health of a child born today is not defined by a changing climate. Lancet 394, 1836–1878 (2019).PubMed 

Google Scholar 
15.Alesheikh, A. A., Ghorbanali, A. & Nouri, N. Coastline change detection using remote sensing. Int. J. Environ. Sci. Technol. 4, 61–66 (2007).
Google Scholar 
16.Specter, C. & Gayle, D. Managing technology transfer for coastal zone development: Caribbean experts identify major issues. Int. J. Remote Sens. 11, 1729–1740 (1990).
Google Scholar 
17.Green, E. P., Mumby, P. J., Edwards, A. J. & Clark, C. D. A review of remote sensing for the assessment and management of tropical coastal resources. Coast. Manag. 24, 1–40 (1996).
Google Scholar 
18.NASA. MODIS (Moderate Resolution Imaging Spectroradiometer). https://modis.gsfc.nasa.gov/about/ (2021).19.Kilpatrick, K. A. et al. A decade of sea surface temperature from MODIS. Remote Sens. Environ. 165, 27–41 (2015).ADS 

Google Scholar 
20.Esaias, W. E. et al. An overview of MODIS capabilities for ocean science observations. IEEE Trans. Geosci. Remote Sens. 36, 1250–1265 (1998).ADS 

Google Scholar 
21.Donlon, C. J. et al. Toward improved validation of satellite SST measurements for climate research. J. Clim. 15, 353–369 (2002).ADS 

Google Scholar 
22.Minnett, P. J. Satellite infrared scanning radiometers — AVHRR and ATSR/M. in Microwave Remote Sensing for Oceanographic and Marine Weather-Forecast Models 141–163 (Springer Netherlands, 1990).23.Hillger, D. et al. First-Light Imagery from Suomi NPP VIIRS. Bull. Am. Meteorol. Soc. 94, 1019–1029 (2013).ADS 

Google Scholar 
24.O’Brien, J. From MODIS to VIIRS – Making the Switch for Air Quality Professionals. NASA Earth Science/Applied Science https://appliedsciences.nasa.gov/our-impact/news/modis-viirs-making-switch-air-quality-professionals (2020).25.Minnett, P. J., Evans, R. H., Podestá, G. P. & Kilpatrick, K. A. Sea-surface temperature from Suomi-NPP VIIRS: Algorithm development and uncertainty estimation. in SPIE 9111, Ocean Sensing and Monitoring VI (eds. Hou, W. W. & Arnone, R. A.) 91110C (2014).26.Drusch, M. et al. Sentinel-2: ESA’s Optical High-Resolution Mission for GMES Operational Services. Remote Sens. Environ. 120, 25–36 (2012).ADS 

Google Scholar 
27.Donlon, C. et al. The global ocean data assimilation experiment high-resolution sea surface temperature pilot project. Bull. Am. Meteorol. Soc. 88, 1197–1214 (2007).ADS 

Google Scholar 
28.NOAA. Ocean Facts: Why do scientists measure sea surface temperature? https://oceanservice.noaa.gov/facts/sea-surface-temperature.html (2020).29.Wei, G. F., Tang, D. L. & Wang, S. Distribution of chlorophyll and harmful algal blooms (HABs): A review on space based studies in the coastal environments of Chinese marginal seas. Adv. Sp. Res. 41, 12–19 (2008).ADS 
CAS 

Google Scholar 
30.O’Reilly, J. E. et al. Ocean color chlorophyll algorithms for SeaWiFS. J. Geophys. Res. Ocean. 103, 24937–24953 (1998).ADS 

Google Scholar 
31.Hu, C., Lee, Z. & Franz, B. Chlorophyll a algorithms for oligotrophic oceans: A novel approach based on three-band reflectance difference. J. Geophys. Res. Ocean. 117, C01011 (2012).ADS 

Google Scholar 
32.Vezzulli, L. et al. Climate influence on Vibrio and associated human diseases during the past half-century in the coastal North Atlantic. Proc. Natl. Acad. Sci. 113, E5062–E5071 (2016).PubMed 
PubMed Central 
CAS 

Google Scholar 
33.Lipp, E. K., Huq, A. & Colwell, R. R. Effects of global climate on infectious disease: The Cholera model. Clin. Microbiol. Rev. 15, 757–770 (2002).PubMed 
PubMed Central 

Google Scholar 
34.Grimes, J. D. et al. Viewing marine bacteria, their activity and response to environmental drivers from orbit: Satellite remote sensing of bacteria. Microb. Ecol. 67, 489–500 (2014).PubMed 
PubMed Central 

Google Scholar 
35.Shen, L., Xu, H. & Guo, X. Satellite remote sensing of harmful algal blooms (HABs) and a potential synthesized framework. Sensors 12, 7778–803 (2012).36.Hayashi, M., Jin, F. & Stuecker, M. F. Dynamics for El Niño-La Niña asymmetry constrain equatorial-Pacific warming pattern. Nat. Commun. 11, 1–10 (2020).
Google Scholar 
37.Hughes, T. P. et al. Spatial and temporal patterns of mass bleaching of corals in the Anthropocene. Science 359, 80–83 (2018).38.Minnett, P. J. et al. Sea-surface temperature measurements from the moderate-resolution imaging spectroradiometer (MODIS) on Aqua and Terra. in IEEE International Geoscience and Remote Sensing Symposium Proceedings. 2004 7, 4576–4579 (2004).39.Minnett, P. J. The validation of sea surface temperature retrievals from spaceborne infrared radiometers. in Oceanography from Space (Springer Netherlands, 2010).40.Minnett, P. J. & Corlett, G. K. A pathway to generating climate data records of sea-surface temperature from satellite measurements. Deep Sea Res. Part II Top. Stud. Oceanogr. 77–80, 44–51 (2012).ADS 

Google Scholar 
41.Castaneda-Guzman, M., Mantilla-Saltos, G., Murray, K. A., Settlage, R. & Escobar, L. E. A database of global coastal conditions. Figshare https://doi.org/10.6084/m9.figshare.c.5660263.v1 (2021).42.R Core Team. R: A Language and Environment for Statistical Computing. (2020).43.NOAA. National Oceanic and Atmospheric Administration (NOAA) Coastal Watch. https://coastwatch.pfeg.noaa.gov/erddapinfo/ (2021).44.Castaneda-Guzman, M., Mantilla-Saltos, G., Murray, K. A., Settlage, R. & Escobar, L. E. Methods and code. Figshare https://doi.org/10.6084/m9.figshare.13708642.v4 (2021).45.Stanford. Best practices for file formats. https://library.stanford.edu/research/data-management-services/data-best-practices/best-practices-file-formats (2021).46.UCAR Community Programs. Network Common Data Form (NetCDF). https://www.unidata.ucar.edu/software/netcdf/ (2021).47.Michna, P. & Woods, M. RNetCDF: Interface to ‘NetCDF’ Datasets. (2019).48.Hijmans, R. J. raster: Geographic Data Analysis and Modeling. (2020).49.ArcGIS. What is a raster data? https://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/what-is-raster-data.htm (2021).50.United Nations. United Nations Convention on the Law of the Sea. 1833 U.N.T.S. 397 (1982).51.Tilstone, G. H. et al. Assessment of MODIS-Aqua chlorophyll-a algorithms in coastal and shelf waters of the eastern Arabian Sea. Cont. Shelf Res. 65, 14–26 (2013).ADS 

Google Scholar 
52.Hoge, F. E. et al. Validation of Terra-MODIS phytoplankton chlorophyll fluorescence line height. I. Initial airborne Lidar results. Appl. Opt. 42, 2767-2771 (2003).ADS 
PubMed 

Google Scholar 
53.Remer, L. A. Validation of MODIS aerosol retrieval over ocean. Geophys. Res. Lett. 29, 8008 (2002).ADS 

Google Scholar 
54.Gentemann, C. L. Three way validation of MODIS and AMSR-E sea surface temperatures. J. Geophys. Res. Ocean. 119, 2583–2598 (2014).ADS 

Google Scholar 
55.Fang, H., Wei, S. & Liang, S. Validation of MODIS and CYCLOPES LAI products using global field measurement data. Remote Sens. Environ. 119, 43–54 (2012).ADS 

Google Scholar 
56.Hosoda, K., Murakami, H., Sakaida, F. & Kawamura, H. Algorithm and validation of sea surface temperature observation using MODIS sensors aboard terra and aqua in the western North Pacific. J. Oceanogr. 63, 267–280 (2007).
Google Scholar 
57.Hao, Y. et al. Validation of MODIS sea surface temperature product in the coastal waters of the Yellow Sea. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10, 1667–1680 (2017).ADS 

Google Scholar 
58.Sims, D. A. et al. On the use of MODIS EVI to assess gross primary productivity of North American ecosystems. J. Geophys. Res. Biogeosciences 111 (2006).59.Miles, T. N. & He, R. Temporal and spatial variability of Chl-a and SST on the South Atlantic Bight: Revisiting with cloud-free reconstructions of MODIS satellite imagery. Cont. Shelf Res. 30, 1951–1962 (2010).ADS 

Google Scholar 
60.Ma, S., Zhang, X., Ding, C., Han, W. & Lu, Y. Comparison of the spatiotemporal variation of Chl-a in the East China Sea and Bohai Sea based on long time series satellite data. in 2021 9th International Conference on Agro-Geoinformatics (Agro-Geoinformatics) 1–6 (2021).61.Watts, N. et al. The 2020 report of The Lancet Countdown on health and climate change: Responding to converging crises. Lancet 6736 (2020).62.Moradi, M. & Kabiri, K. Spatio-temporal variability of SST and Chlorophyll-a from MODIS data in the Persian Gulf. Mar. Pollut. Bull. 98, 14–25 (2015).PubMed 
CAS 

Google Scholar 
63.Golder, M. R. et al. Chlorophyll-a, SST and particulate organic carbon in response to the cyclone Amphan in the Bay of Bengal. J. Earth Syst. Sci. 130, 157 (2021).ADS 
CAS 

Google Scholar 
64.Minnett, P. J., Evans, R. H., Kearns, E. J. & Brown, O. B. Sea-surface temperature measured by the Moderate Resolution Imaging Spectroradiometer (MODIS). in IEEE International Geoscience and Remote Sensing Symposium vol. 2, 1177–1179 (IEEE, 2002).65.Qin, H., Chen, G., Wang, W., Wang, D. & Zeng, L. Validation and application of MODIS-derived SST in the South China Sea. Int. J. Remote Sens. 35, 4315–4328 (2014).
Google Scholar 
66.Saulquin, B., Gohin, F. & Garrello, R. Regional Objective Analysis for Merging High-Resolution MERIS, MODIS/Aqua, and SeaWiFS Chlorophyll-a Data From 1998 to 2008 on the European Atlantic Shelf. IEEE Trans. Geosci. Remote Sens. 49, 143–154 (2011).ADS 

Google Scholar 
67.Chen, J. & Quan, W. An improved algorithm for retrieving chlorophyll-a from the Yellow River Estuary using MODIS imagery. Environ. Monit. Assess. 185, 2243–2255 (2013).PubMed 

Google Scholar 
68.Hanafin, J. A. & Minnett, P. J. Thermal profiling of the sea surface skin layer using FTIR measurements. in Gas Transfer at Water Surfaces 161–166 (Blackwell Publishing, 2002).69.Wong, E. W. & Minnett, P. J. The response of the ocean thermal skin layer to variations in incident infrared radiation. J. Geophys. Res. Ocean. 123, 2475–2493 (2018).ADS 

Google Scholar 
70.Ward, B. Near-surface ocean temperature. J. Geophys. Res. 111, C02004 (2006).ADS 

Google Scholar 
71.Kilpatrick, K. A., Podestá, G. P. & Evans, R. Overview of the NOAA/NASA advanced very high resolution radiometer Pathfinder algorithm for sea surface temperature and associated matchup database. J. Geophys. Res. Ocean. 106, 9179–9197 (2001).ADS 

Google Scholar 
72.Hollstein, A., Segl, K., Guanter, L., Brell, M. & Enesco, M. Ready-to-use methods for the detection of clouds, cirrus, snow, shadow, water and clear sky pixels in Sentinel-2 MSI images. Remote Sens. 8, 666 (2016).ADS 

Google Scholar 
73.Luo, B., Minnett, P. J., Gentemann, C. & Szczodrak, G. Improving satellite retrieved night-time infrared sea surface temperatures in aerosol contaminated regions. Remote Sens. Environ. 223, 8–20 (2019).ADS 

Google Scholar 
74.Moore, T. S., Campbell, J. W. & Dowell, M. D. A class-based approach to characterizing and mapping the uncertainty of the MODIS ocean chlorophyll product. Remote Sens. Environ. 113, 2424–2430 (2009).ADS 

Google Scholar 
75.Pieri, M. et al. Assessment of three algorithms for the operational estimation of [CHL] from MODIS data in the Western Mediterranean Sea. Eur. J. Remote Sens. 48, 383–401 (2015).
Google Scholar 
76.Tilstone, G. H. et al. Performance of Ocean Colour Chlorophyll-a algorithms for Sentinel-3 OLCI, MODIS-Aqua and Suomi-VIIRS in open-ocean waters of the Atlantic. Remote Sens. Environ. 260, 112444 (2021).ADS 

Google Scholar  More

1.Schmidt MWI, Torn MS, Abiven S, Dittmar T, Guggenberger G, Janssens IA, et al. Persistence of soil organic matter as an ecosystem property. Nature. 2011;478:49–56.CAS 
PubMed 

Google Scholar 
2.Bond-Lamberty B, Bailey VL, Chen M, Gough CM, Vargas R. Globally rising soil heterotrophic respiration over recent decades. Nature. 2018;560:80–3.CAS 
PubMed 

Google Scholar 
3.Bradford MA. Thermal adaptation of decomposer communities in warming soils. Front Microbiol. 2013;4:1–16.
Google Scholar 
4.Cavicchioli R, Ripple WJ, Timmis KN, Azam F, Bakken LR, Baylis M, et al. Scientists’ warning to humanity: microorganisms and climate change. Nat Rev Microbiol. 2019;17:569–86.CAS 
PubMed 
PubMed Central 

Google Scholar 
5.Jansson JK, Hofmockel KS. Soil microbiomes and climate change. Nat Rev Microbiol. 2020;18:35–46.CAS 
PubMed 

Google Scholar 
6.Liu L, Greaver TL. A global perspective on belowground carbon dynamics under nitrogen enrichment. Ecol Lett. 2010;13:819–28.PubMed 

Google Scholar 
7.Knicker H. Soil organic N – An under-rated player for C sequestration in soils? Soil Biol Biochem. 2011;43:1118–29.CAS 

Google Scholar 
8.Soong JL, Fuchslueger L, Marañon-Jimenez S, Torn MS, Janssens IA, Peñuelas J, et al. Microbial carbon limitation: The need for integrating microorganisms into our understanding of ecosystem carbon cycling. Glob Chang Biol. 2020;26:1953–61.
Google Scholar 
9.Mooshammer M, Wanek W, Hämmerle I, Fuchslueger L, Hofhansl F, Knoltsch A, et al. Adjustment of microbial nitrogen use efficiency to carbon:nitrogen imbalances regulates soil nitrogen cycling. Nat Commun. 2014;5:1–7.
Google Scholar 
10.Geisseler D, Horwath WR, Joergensen RG, Ludwig B. Pathways of nitrogen utilization by soil microorganisms – a review. Soil Biol Biochem. 2010;42:2058–67.CAS 

Google Scholar 
11.Wang X, Wang C, Cotrufo MF, Sun L, Jiang P, Liu Z, et al. Elevated temperature increases the accumulation of microbial necromass nitrogen in soil via increasing microbial turnover. Glob Chang Biol. 2020;26:5277–89.PubMed 

Google Scholar 
12.Simpson AJ, Simpson MJ, Smith E, Kelleher BP. Microbially derived inputs to soil organic matter: Are current estimates too low? Environ Sci Technol. 2007;41:8070–6.CAS 
PubMed 

Google Scholar 
13.Kuypers MMM, Marchant HK, Kartal B. The microbial nitrogen-cycling network. Nat Rev Microbiol. 2018;16:263–76.CAS 
PubMed 

Google Scholar 
14.Walker TWN, Kaiser C, Strasser F, Herbold CW, Leblans NIW, Woebken D, et al. Microbial temperature sensitivity and biomass change explain soil carbon loss with warming. Nat Climate Change. 2018;8:885–9.CAS 

Google Scholar 
15.Marañón-Jiménez S, Peñuelas J, Richter A, Sigurdsson BD, Fuchslueger L, Leblans NIW, et al. Coupled carbon and nitrogen losses in response to seven years of chronic warming in subarctic soils. Soil Biol Biochem. 2019;134:152–61.
Google Scholar 
16.Nguyen TTH, Myrold DD, Mueller RS. Distributions of extracellular peptidases across prokaryotic genomes reflect phylogeny and habitat. Front Microbiol. 2019;10:1–14.
Google Scholar 
17.Zimmerman AE, Martiny AC, Allison SD. Microdiversity of extracellular enzyme genes among sequenced prokaryotic genomes. ISME J. 2013;7:1187–99.CAS 
PubMed 
PubMed Central 

Google Scholar 
18.Beier S, Bertilsson S. Bacterial chitin degradation-mechanisms and ecophysiological strategies. Front Microbiol. 2013;4:1–12.
Google Scholar 
19.Kielak AM, Cretoiu MS, Semenov AV, Sørensen SJ, Van, Elsas JD. Bacterial chitinolytic communities respond to chitin and pH alteration in soil. Appl Environ Microbiol. 2013;79:263–72.CAS 
PubMed 
PubMed Central 

Google Scholar 
20.Weintraub MN, Schimel JP. Seasonal protein dynamics in Alaskan arctic tundra soils. Soil Biol Biochem. 2005;37:1469–75.CAS 

Google Scholar 
21.Boer VM, De Winde JH, Pronk JT, Piper MDW. The genome-wide transcriptional responses of Saccharomyces cerevisiae grown on glucose in aerobic chemostat cultures limited for carbon, nitrogen, phosphorus, or sulfur. J Biol Chem. 2003;278:3265–74.CAS 
PubMed 

Google Scholar 
22.Kolkman A, Daran-Lapujade P, Fullaondo A, Olsthoorn MMA, Pronk JT, Slijper M, et al. Proteome analysis of yeast response to various nutrient limitations. Mol Syst Biol. 2006;2:1–16.
Google Scholar 
23.Silberbach M, Hüser A, Kalinowski J, Pühler A, Walter B, Krämer R, et al. DNA microarray analysis of the nitrogen starvation response of Corynebacterium glutamicum. J Biotechnol. 2005;119:357–67.CAS 
PubMed 

Google Scholar 
24.Merrick MJ, Edwards RA. Nitrogen control in bacteria. Microbiol Rev. 1995;59:604–22.CAS 
PubMed 
PubMed Central 

Google Scholar 
25.Daebeler A, Abell GCJ, Bodelier PLE, Bodrossy L, Frampton DMF, Hefting MM, et al. Archaeal dominated ammonia-oxidizing communities in Icelandic grassland soils are moderately affected by long-term N fertilization and geothermal heating. Front Microbiol. 2012;3:1–14.
Google Scholar 
26.Yeager CM, Kornosky JL, Housman DC, Grote EE, Belnap J, Kuske CR. Diazotrophic community structure and function in two successional stages of biological soil crusts from the colorado plateau and Chihuahuan Desert. Appl Environ Microbiol. 2004;70:973–83.CAS 
PubMed 
PubMed Central 

Google Scholar 
27.Malik AA, Swenson T, Weihe C, Morrison EW, Martiny JBH, Brodie EL, et al. Drought and plant litter chemistry alter microbial gene expression and metabolite production. ISME J. 2020;14:2236–47.CAS 
PubMed 
PubMed Central 

Google Scholar 
28.Tveit A, Schwacke R, Svenning MM, Urich T. Organic carbon transformations in high-Arctic peat soils: Key functions and microorganisms. ISME J. 2013;7:299–311.CAS 
PubMed 

Google Scholar 
29.Geisen S, Tveit AT, Clark IM, Richter A, Svenning MM, Bonkowski M, et al. Metatranscriptomic census of active protists in soils. ISME J. 2015;9:2178–90.CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Urich T, Lanzén A, Qi J, Huson DH, Schleper C, Schuster SC. Simultaneous assessment of soil microbial community structure and function through analysis of the meta-transcriptome. PLoS ONE. 2008;3:1–13.
Google Scholar 
31.Kallenbach CM, Frey SD, Grandy AS. Direct evidence for microbial-derived soil organic matter formation and its ecophysiological controls. Nat Commun. 2016;7:1–10.
Google Scholar 
32.Walker TWN, Janssens IA, Weedon JT, Sigurdsson BD, Richter A, Peñuelas J, et al. A systemic overreaction to years versus decades of warming in a subarctic grassland ecosystem. Nat Ecol Evol. 2020;4:101–8.PubMed 

Google Scholar 
33.Sigurdsson BD, Wallander H, Gunnarsdóttir GE, Richter A, Sigurðsson P, Leblans NIW, et al. Geothermal ecosystems as natural climate change experiments: the ForHot research site in Iceland as a case study. Icelandic Agric Sci. 2016;29:53–71.
Google Scholar 
34.Söllinger A, Séneca J, Dahl MB, Prommer J, Verbruggen E, Sigurdsson BD, et al. Downregulation of the microbial protein biosynthesis machinery in response to weeks, years and decades of soil warming. 2021 Research Square preprint. https://doi.org/10.21203/rs.3.rs-132190/v235.Leblans N. Natural gradients in temperature and nitrogen: Iceland represents a unique environment to clarify long-term global change effects on carbon dynamics. Joint doctoral dissertation. Antwerp University and Agricultural University of Iceland, Reykjavik, Iceland; 2016:1–229.36.Angel R, Claus P, Conrad R. Methanogenic archaea are globally ubiquitous in aerated soils and become active under wet anoxic conditions. ISME J. 2012;6:847–62.CAS 
PubMed 

Google Scholar 
37.Hyatt D, Chen G-L, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:5–11.
Google Scholar 
38.Gillespie CS. Fitting heavy tailed distributions: the poweRlaw Package. J Stat Softw. 2015;64:1–16.
Google Scholar 
39.El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2018;47:427–32.
Google Scholar 
40.Eddy SR. Accelerated profile HMM searches. PLOS Comput Biol. 2011;7:e1002195.CAS 
PubMed 
PubMed Central 

Google Scholar 
41.Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat Methods. 2011;8:785–6.CAS 
PubMed 

Google Scholar 
42.Bendtsen JD, Kiemer L, Fausbøll A, Brunak S. Non-classical protein secretion in bacteria. BMC Microbiol. 2005;5:1–13.
Google Scholar 
43.Yu NY, Wagner JR, Laird MR, Melli G, Rey S, Lo R, et al. PSORTb 3.0: improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes. Bioinformatics. 2010;26:1608–15.CAS 
PubMed 
PubMed Central 

Google Scholar 
44.Orsi WD. MetaProt: an integrated database of predicted proteins for improved annotation of metaomic datasets. Open Data LMU. 2020. https://doi.org/10.5282/ubm/data.18345.Lombard V, Golaconda Ramulu H, Drula E, Coutinho PM, Henrissat B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 2013;42:490–5.
Google Scholar 
46.Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59–60.PubMed 

Google Scholar 
47.Oksanen AJ, Blanchet FG, Kindt R, Legen- P, Minchin PR, Hara RBO, et al. vegan: Community Ecology Package. 2019. https://cran.r-project.org/package=vegan48.Lê S, Josse J, Husson F. FactoMineR: an R package for multivariate analysis. J Stat Softw. 2008;25:1–18.
Google Scholar 
49.Kolde R. pheatmap: pretty heatmaps. 2019. https://cran.r-project.org/package=pheatmap50.Noll L, Zhang S, Zheng Q, Hu Y, Wanek W. Wide-spread limitation of soil organic nitrogen transformations by substrate availability and not by extracellular enzyme content. Soil Biol Biochem. 2019;133:37–49.CAS 
PubMed 
PubMed Central 

Google Scholar 
51.Schimel JP, Bennett J. Nitrogen mineralization: challenges of a changing paradigm. Ecology. 2004;85:591–602.
Google Scholar 
52.Wild B, Ambus P, Reinsch S, Richter A. Resistance of soil protein depolymerization rates to eight years of elevated CO2, warming, and summer drought in a temperate heathland. Biogeochemistry. 2018;140:255–67.CAS 

Google Scholar 
53.Wanek W, Mooshammer M, Blöchl A, Hanreich A, Richter A. Determination of gross rates of amino acid production and immobilization in decomposing leaf litter by a novel 15N isotope pool dilution technique. Soil Biol Biochem. 2010;42:1293–302.CAS 

Google Scholar 
54.Liang C, Schimel JP, Jastrow JD. The importance of anabolism in microbial control over soil carbon storage. Nat Microbiol. 2017;2:1–6.
Google Scholar 
55.Vranova V, Rejsek K, Formanek P. Proteolytic activity in soil: a review. Appl Soil Ecol. 2013;70:23–32.
Google Scholar 
56.Schimel JP, Weintraub MN. The implications of exoenzyme activity on microbial carbon and nitrogen limitation in soil: a theoretical model. Soil Biol Biochem. 2003;35:549–63.CAS 

Google Scholar 
57.Rawlings ND, Waller M, Barrett AJ, Bateman A. MEROPS: The database of proteolytic enzymes, their substrates and inhibitors. Nucleic Acids Res. 2014;42:503–9.
Google Scholar 
58.Vollmer W, Joris B, Charlier P, Foster S. Bacterial peptidoglycan (murein) hydrolases. FEMS Microbiol Rev. 2008;32:259–86.CAS 
PubMed 

Google Scholar 
59.Vermassen A, Leroy S, Talon R, Provot C, Popowska M, Desvaux M. Cell wall hydrolases in bacteria: Insight on the diversity of cell wall amidases, glycosidases and peptidases toward peptidoglycan. Front Microbiol. 2019;10:1–27.
Google Scholar 
60.Donhauser J, Qi W, Bergk-Pinto B, Frey B. High temperatures enhance the microbial genetic potential to recycle C and N from necromass in high-mountain soils. Glob Chang Biol. 2020;27:1365–86.61.Vollmer W, Blanot D, De Pedro MA. Peptidoglycan structure and architecture. FEMS Microbiology Reviews. 2008;32:149–67.CAS 
PubMed 

Google Scholar 
62.Semchenko M, Leff JW, Lozano YM, Saar S, Davison J, Wilkinson A, et al. Fungal diversity regulates plant-soil feedbacks in temperate grassland. Science Adv. 2018;4.63.Saary P, Mitchell AL, Finn RD. Estimating the quality of eukaryotic genomes recovered from metagenomic analysis. Genome Biol. 2020;21:244.PubMed 
PubMed Central 

Google Scholar  More

1.Russell JB, Cook GM. Energetics of bacterial growth: balance of anabolic and catabolic reactions. Microbiol Rev. 1995;59:48–62.CAS 
PubMed 
PubMed Central 

Google Scholar 
2.Klumpp S, Scott M, Pedersen S, Hwa T. Molecular crowding limits translation and cell growth. PNAS. 2013;110:16754–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
3.Pedersen S. Escherichia coli ribosomes translate in vivo with variable rate. EMBO J. 1984;3:2895–8.CAS 
PubMed 
PubMed Central 

Google Scholar 
4.Ran W, Higgs PG. Contributions of speed and accuracy to translational selection in bacteria. PLoS One. 2012;7:e51652.CAS 
PubMed 
PubMed Central 

Google Scholar 
5.Vieira-Silva S, Rocha E. The systemic imprint of growth and its uses in ecological (meta)genomics. PLoS Genet. 2009;6:1–15.
Google Scholar 
6.Roller BRK, Stoddard SF, Schmidt TM. Exploiting rRNA operon copy number to investigate bacterial reproductive strategies. Nat Microbiol. 2016;1:1–7.
Google Scholar 
7.Buskirk AR, Green R. Ribosome pausing, arrest and rescue in bacteria and eukaryotes. Philos Trans R Soc B. 2017;372:20160183–11.
Google Scholar 
8.Wohlgemuth I, Brenner S, Beringer M, Rodnina MV. Modulation of the rate of peptidyl transfer on the ribosome by the nature of substrates. J Biol Chem. 2008;283:32229–35.CAS 
PubMed 

Google Scholar 
9.Pavlov MY, Watts RE, Tan Z, Cornish VW, Ehrenberg M, Forster AC. Slow peptide bond formation by proline and other N-alkylamino acids in translation. PNAS. 2009;106:50–54.CAS 
PubMed 

Google Scholar 
10.Mandal A, Mandal S, Park MH. Genome-wide analyses and functional classification of proline repeat-rich proteins: potential role of eIF5A in eukaryotic evolution. PLoS One. 2014;9:e111800–13.PubMed 
PubMed Central 

Google Scholar 
11.Adzhubei AA, Sternberg MJE, Makarov AA. Polyproline-II helix in proteins: structure and function. J Mol Biol. 2013;425:2100–32.CAS 
PubMed 

Google Scholar 
12.Elam WA, Schrank TP, Campagnolo AJ, Hilser VJ. Evolutionary conservation of the polyproline II conformation surrounding intrinsically disordered phosphorylation sites. Protein Sci. 2013;22:405–17.PubMed 

Google Scholar 
13.Ball LJ, Kühne R, Schneider-Mergener J, Oschkinat H. Recognition of proline-rich motifs by protein-protein-interaction domains. Angew Chem Int Ed Engl. 2005;44:2852–69.CAS 
PubMed 

Google Scholar 
14.Starosta AL, Lassak J, Peil L, Atkinson GC, Virumäe K, Tenson T, et al. Translational stalling at polyproline stretches is modulated by the sequence context upstream of the stall site. Nucleic Acids Res. 2014;42:10711–9.CAS 
PubMed 
PubMed Central 

Google Scholar 
15.Woolstenhulme CJ, Guydosh NR, Green R, Buskirk AR. High-precision analysis of translational pausing by ribosome profiling in bacteria lacking EFP. Cell Rep. 2015;11:13–21.CAS 
PubMed 
PubMed Central 

Google Scholar 
16.Hersch SJ, Elgamal S, Katz A, Ibba M, Navarre WW. Translation initiation rate determines the impact of ribosome stalling on bacterial protein synthesis. J Biol Chem. 2014;289:28160–71.CAS 
PubMed 
PubMed Central 

Google Scholar 
17.Lassak J, Wilson DN, Jung K. Stall no more at polyproline stretches with the translation elongation factors EF‐P and IF‐5A. Mol Microbiol. 2016;99:219–35.CAS 
PubMed 

Google Scholar 
18.Yanagisawa T, Sumida T, Ishii R, Takemoto C, Yokoyama S. A paralog of lysyl-tRNA synthetase aminoacylates a conserved lysine residue in translation elongation factor P. Nature. 2010;17:1136–43.CAS 

Google Scholar 
19.Park J-H, Johansson HE, Aoki H, Huang BX, Kim H-Y, Ganoza MC, et al. Post-translational modification by beta-lysylation is required for activity of Escherichia coli elongation factor P (EF-P). J Biol Chem. 2012;287:2579–90.CAS 
PubMed 

Google Scholar 
20.Lassak J, Keilhauer E, Fürst M, Wuichet K, Gödeke J, Starosta AL, et al. Arginine-rhamnosylation as new strategy to activate translation elongation factor P. Nat Chem Biol. 2015;11:266–70.CAS 
PubMed 
PubMed Central 

Google Scholar 
21.Tollerson R, Witzky A, Ibba M. Elongation factor P is required to maintain proteome homeostasis at high growth rate. PNAS. 2018;115:1–6.
Google Scholar 
22.Peng WT, Banta LM, Charles TC, Nester EW. The chvH locus of Agrobacterium encodes a homologue of an elongation factor involved in protein synthesis. J Bacteriol. 2001;183:36–45.CAS 
PubMed 
PubMed Central 

Google Scholar 
23.Rajkovic A, Hummels KR, Witzky A, Erickson S, Gafken PR, Whitelegge JP, et al. Translation control of swarming proficiency in Bacillus subtilis by 5-amino-pentanolylated elongation factor P. J Biol Chem. 2016;291:10976–85.CAS 
PubMed 
PubMed Central 

Google Scholar 
24.Navarre WW, Zou SB, Roy H, Xie JL, Savchenko A, Singer A, et al. PoxA, YjeK, and elongation factor P coordinately modulate virulence and drug resistance in Salmonella enterica. Mol Cell. 2010;39:209–21.CAS 
PubMed 
PubMed Central 

Google Scholar 
25.Hummels KR, Kearns DB. Suppressor mutations in ribosomal proteins and FliY restore Bacillus subtilis swarming motility in the absence of EF-P. PLoS Genet. 2019;15:e1008179–27.CAS 
PubMed 
PubMed Central 

Google Scholar 
26.Rajkovic A, Erickson S, Witzky A, Branson OE, Seo J, Gafken PR, et al. Cyclic rhamnosylated elongation factor P establishes antibiotic resistance in Pseudomonas aeruginosa. MBio. 2015;6:1–9.
Google Scholar 
27.Yanagisawa T, Takahashi H, Suzuki T, Masuda A, Dohmae N, Yokoyama S. Neisseria meningitidis translation elongation factor P and its active-site arginine residue are essential for cell viability. PLoS One. 2016;11:e0147907–27.PubMed 
PubMed Central 

Google Scholar 
28.Krafczyk R, Qi F, Sieber A, Mehler J, Jung K, Frishman D, et al. Proline codon pair selection determines ribosome pausing strength and translation efficiency in bacteria. Commun Biol. 2021;4:1–11.
Google Scholar 
29.Qi F, Motz M, Jung K, Lassak J, Frishman D. Evolutionary analysis of polyproline motifs in Escherichia coli reveals their regulatory role in translation. PLoS Comput Biol. 2018;14:e1005987–19.PubMed 
PubMed Central 

Google Scholar 
30.Karlin S, Mrázek J, Campbell A, Kaiser D. Characterizations of highly expressed genes of four fast-growing bacteria. J Bacteriol. 2001;183:5025–40.CAS 
PubMed 
PubMed Central 

Google Scholar 
31.Dethlefsen L, Schmidt TM. Performance of the translational apparatus varies with the ecological strategies of bacteria. J Bacteriol. 2007;189:3237–45.CAS 
PubMed 
PubMed Central 

Google Scholar 
32.Weissman JL, Hou S, Fuhrman JA. Estimating maximal microbial growth rates from cultures, metagenomes, and single cells via codon usage patterns. PNAS. 2021;118:1–10.
Google Scholar 
33.Hersch SJ, Wang M, Zou SB, Moon K-M, Foster LJ, Ibba M, et al. Divergent protein motifs direct elongation factor P-mediated translational regulation in Salmonella enterica and Escherichia coli. MBio. 2013;4:1–10.
Google Scholar 
34.Pinheiro B, Scheidler CM, Kielkowski P, Schmid M, Forné I, Ye S, et al. Structure and function of an elongation factor P subfamily in Actinobacteria. Cell Rep. 2020;30:4332–42. e5.CAS 
PubMed 

Google Scholar 
35.Chen I-MA, Markowitz VM, Chu K, Palaniappan K, Szeto E, Pillay M, et al. IMG/M: integrated genome and metagenome comparative data analysis system. Nucleic Acids Res. 2017;45:D507–D516.CAS 
PubMed 

Google Scholar 
36.Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.CAS 
PubMed 
PubMed Central 

Google Scholar 
37.Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics. 2020;36:1925–7.CAS 

Google Scholar 
38.Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7:e1002195–16.CAS 
PubMed 
PubMed Central 

Google Scholar 
39.Madin JS, Nielsen DA, Brbic M, Corkrey R, Danko D, Edwards K, et al. A synthesis of bacterial and archaeal phenotypic trait data. Sci Data. 2020;7:1–8.
Google Scholar 
40.Novembre JA. Accounting for background nucleotide composition when measuring codon usage bias. Mol Biol Evol. 2002;19:1390–4.CAS 
PubMed 

Google Scholar 
41.Erdos G, Dosztányi Z. Analyzing protein disorder with IUPred2A. Curr Protoc Bioinforma. 2020;70:1–15.
Google Scholar 
42.Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.CAS 
PubMed 
PubMed Central 

Google Scholar 
43.Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol. 2009;26:1641–50.CAS 
PubMed 
PubMed Central 

Google Scholar 
44.Pagel M. Inferring the historical patterns of biological evolution. Nature. 1999;401:877–84.CAS 
PubMed 

Google Scholar 
45.Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2011;3:217–23.
Google Scholar 
46.Orme D, Freckleton R, Thomas G, Petzoldt T, Fritz S, Isaac N, et al. caper: comparative analyses of phylogenetics and evolution in R. 2018; https://CRAN.R-project.org/package=caper.47.Wickham H. ggplot2: elegant graphics for data analysis. 2016. Springer-Verlag New York.48.Symonds MRE, Blomberg SP. A primer on phylogenetic generalized least squares. In: Garamszegi L (eds). Modern phylogenetic comparative methods and their application in evolutionary biology. (Springer, Berlin, Heidelberg, 2014) pp. 105–30.49.Watanabe K, Suzuki Y. Protein thermostabilization by proline substitutions. J Mol Catal B Enzym. 1998;4:167–80.CAS 

Google Scholar 
50.Sabath N, Ferrada E, Barve A, Wagner A. Growth temperature and genome size in bacteria are negatively correlated, suggesting genomic streamlining during thermal adaptation. Genome Biol Evol. 2013;5:966–77.PubMed 
PubMed Central 

Google Scholar 
51.Goldman BS, Nierman WC, Kaiser D, Slater SC, Durkin AS, Eisen JA, et al. Evolution of sensory complexity recorded in a myxobacterial genome. PNAS. 2006;103:15200–5.CAS 
PubMed 
PubMed Central 

Google Scholar 
52.Long AM, Hou S, Ignacio-Espinoza JC, Fuhrman JA. Benchmarking microbial growth rate predictions from metagenomes. ISME J. 2020;15:1–13.
Google Scholar 
53.Rocha EPC. Codon usage bias from tRNA’s point of view: redundancy, specialization, and efficient decoding for translation optimization. Genome Res. 2004;14:2279–86.CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Klappenbach JA, Dunbar JM, Schmidt TM. rRNA operon copy number reflects ecological strategies of bacteria. Appl Environ Microbiol. 2000;66:1328–33.CAS 
PubMed 
PubMed Central 

Google Scholar 
55.Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40:D109–D114.CAS 
PubMed 

Google Scholar 
56.Perez J, Castaneda-García A, Jenke-Kodama H, Muller R, Munoz-Dorado J. Eukaryotic-like protein kinases in the prokaryotes and the myxobacterial kinome. PNAS. 2008;105:15950–5.CAS 
PubMed 
PubMed Central 

Google Scholar 
57.Shi L, Pigeonneau N, Ravikumar V, Dobrinic P, Macek B, Franjevic D, et al. Cross-phosphorylation of bacterial serine/threonine and tyrosine protein kinases on key regulatory residues. Front Microbiol. 2014;5:1–13.CAS 

Google Scholar 
58.Jakob U, Kriwacki R, Uversky VN. Conditionally and transiently disordered proteins: awakening cryptic disorder to regulate protein function. Chem Rev. 2014;114:6779–805.CAS 
PubMed 
PubMed Central 

Google Scholar 
59.Starosta AL, Lassak J, Peil L, Atkinson GC, Woolstenhulme CJ, Virumäe K, et al. A conserved proline triplet in Val-tRNA synthetase and the origin of elongation factor P. Cell Rep. 2014;9:476–83.CAS 
PubMed 
PubMed Central 

Google Scholar 
60.Nariya H, Inouye S. A protein Ser/Thr kinase cascade negatively regulates the DNA-binding activity of MrpC, a smaller form of which may be necessary for the Myxococcus xanthus development. Mol Microbiol. 2006;60:1205–17.CAS 
PubMed 

Google Scholar 
61.Stein EA, Cho K, Higgs PI, Zusman DR. Two Ser/Thr protein kinases essential for efficient aggregation and spore morphogenesis in Myxococcus xanthus. Mol Microbiol. 2006;60:1414–31.CAS 
PubMed 

Google Scholar 
62.Iakoucheva LM, Radivojac P, Brown CJ, OConnor TR, Sikes JG, Obradovic Z, et al. The importance of intrinsic disorder for protein phosphorylation. Nucleic Acids Res. 2004;32:1037–49.CAS 
PubMed 
PubMed Central 

Google Scholar 
63.Elsen S, Swem LR, Swem DL, Bauer CE. RegB/RegA, a highly conserved redox-responding global two-component regulatory system. Microbiol Mol Biol R. 2004;68:263–79.CAS 

Google Scholar 
64.Tawa P, Stewart RC. Kinetics of CheA autophosphorylation and dephosphorylation reactions. Biochemistry. 1994;33:7917–24.CAS 
PubMed 

Google Scholar 
65.Yoshida T, jian CaiS, Inouye M. Interaction of EnvZ, a sensory histidine kinase, with phosphorylated OmpR, the cognate response regulator. Mol Microbiol. 2002;46:1283–94.CAS 
PubMed 

Google Scholar 
66.Cho M-H, Wrabl JO, Taylor J, Hilser VJ. Hidden dynamic signatures drive substrate selectivity in the disordered phosphoproteome. PNAS. 2020;117:1–11.
Google Scholar  More

Study sites: site description and climatic limitIn this study, we focused on seven flux tower sites located in Alaska, United States, including six AmeriFlux sites and one site of the Korea Polar Research Institute (KOPRI) (Fig. S1A, Table S1). Over the seven study sites, the annual mean temperature was between −10.09 and −0.55 °C and the annual total precipitation ranged from 287 to 540 during 2001–2018 based on the North American Regional Reanalysis (NARR41, 0.3-degree resolution every 3 h). The annual mean temperature increased from 2001 to 2018 at all sites, with rates between 0.5 and 2.2 °C per decade (p  0.05 at all sites). Our study sites were mostly dominated by wet sedges, grasses, moss, lichens, and dwarf shrubs. For example, the dominant plants at the US-Atq site (at a higher latitude) are herbaceous sedges (Carex aquatilis, Eriophorum russeolum, and Eriophorum angustifolium) and shrubs (Salix rotundifolia), with abundant mosses (Calliergon richardsonii and Cinclidium subrotundum) and lichens (Peltigera sp.)11. At the KOPRI site (at a lower latitude), mosses (Sphagnum magellanicum, Sphagnum angustifolium, and Sphagnum fuscum), lichens (Cladonia mitis, Cladonia crispata, and Cladonia stellaris), and tundra tussock cottongrass (Eriophorum vaginatum) are abundant39. The active layer thickness is between 0.33 and 1.0 m, according to field data and radar-based estimates.Climatic limits imposed by temperature, water, and radiation were quantified following Nemani et al.42 at each site during the GS between 2001 and 2018 (Fig. S2) using the NARR data. In this study, we defined the GS as from May to Oct., early GS as between May and Jun., peak GS as between Jul. and Aug., and late GS as between Sep. and Oct. For a temperature limit scalar, the monthly mean temperature from −5  to 5 °C was linearly scaled between 100% (i.e., no growth) and 0% (i.e., no reduction in growing days). The monthly ratio of precipitation to potential evapotranspiration (PET by the Priestley-Taylor method43), ranging between 0 and 0.75, was linearly scaled from 100 to 0% as a water limit scalar. A radiation limit scalar was estimated as a 0.5% reduction in growing days for every 1% increase in monthly cloudiness above the 10% threshold (monthly cloudiness (n) was estimated44 as (R={R}_{0}(1-0.75{n}^{3.4})), where R and R0 are the monthly mean incoming radiation and clear-sky radiation45, respectively).The carbon flux response to climatic variations at each site was further analyzed using a forward stepwise multiple regression analysis11 between the NEE and meteorological variables (temperature, PAR, and VPD) using tower data during the GS. Interaction terms among the variables are also included to consider the convolved effects of the variables (Eq. (1)).$${Y}_{{{NEE}}}= {beta }_{0}+{beta }_{1}{X}_{{{T}}}+{beta }_{2},{X}_{{{VPD}}}+{beta }_{3}{X}_{{{PAR}}}+{beta }_{4}{X}_{{{T}}}* {X}_{{{VPD}}}+ldots$$$$qquad;;; {beta }_{5}{X}_{{{T}}}* {X}_{{{PAR}}}+{beta }_{6}{X}_{{{VPD}}}* {X}_{{{PAR}}}+{beta }_{7}{X}_{{{T}}}* {X}_{{{VPD}}}* {X}_{{{PAR}}}$$
(1)
where YNEE is the daily average NEE (µmol m−2 s−1), and XT, XVPD, and XPAR are daily average air temperature (°C), VPD (ha), and PAR (µmol Photon m−2 s−1), respectively. Regression coefficients (({beta }_{0},ldots ,{beta }_{7})), standard errors, significance (P-value), and R2 value of the final regression model are summarized in Table S2.MODIS: long-term trends of snowmelt and greenup timingsWe collected the gridded MODIS snow cover (MOD10A1.V00646 at a 500-m resolution every day) and phenology (MCD12Q2.V00634 at a 500-m resolution yearly) from the NASA Earthdata (https://earthdata.nasa.gov/). We estimated the snowmelt and snowpack timings at each site as the date when a logistic fit to the MODIS snow cover (quality flags of good and best) passed 0.1 each year. We rejected those snowmelt timings when the gaps in the daily MODIS snow cover were longer than 2 weeks around the time of the snowmelt event. The greenup and dormancy timings with a quality flag of best were taken from the MODIS phenology. Based on the spatial representativeness assessment (see Supplementary Note), we decided to use the snowmelt timing and greenup timing within a 1 × 1 pixel window.The significance of the long-term trends in greenup and snowmelt timings at each site was determined by Spearman’s rho and Mann-Kendall tests (Fig. S3). We further estimated the 95% confidence intervals of the trends from 3000 timing sets generated by bootstrap resampling from a normal distribution47 (mean equal to each greenup or snowmelt timing with three standard deviation set to 10 or 6.6 days, respectively, i.e., the root-mean-squared values between the ground data-based estimates and MODIS values in a 1 × 1 pixel window, Figs. S9 and S10).SGSI: snowmelt-growing season indexGrowing season index (GSI)48 is one of the novel phenology models49 and has been widely applied for the phenological representations of various ecosystems50,51. GSI is a product of three indices of climatic variables (Eq. (2), Fig. S5): daylength, VPD, and growing-degree-days (GDD)52. As a phenological measure for a given meteorological condition, we calculated the daily GSI for spring (from Jan. 1 to Jul. 31) and fall (from Aug. 1 to Dec. 31), respectively. For the spring-GSI, GDD is the degree sum when the daily mean temperature rises above −5 °C after Jan. 1. For the fall-GSI, GDD is the degree sum when the daily mean temperature falls below 20 °C after Aug. 1. We then revised the GSI by multiplying it by a snowmelt index (iS) and referred to this as the snowmelt-GSI (SGSI, Eq. (3), Fig. S5). This guarantees that vegetation greenup does not start unless snow is melted, even if the meteorological conditions are sufficient to trigger leaf-out. The iS was estimated to be 0 when the snow cover fraction (snowfac variable53 in ED2) was above 0.1 and 1 otherwise.$${{{GSI}}}={{iX}}_{1} times {{iX}}_{2} times {{iX}}_{3}$$
(2)
$${{{SGSI}}}={{{GSI}}} times {iS}$$
(3)
where iX (X1, X2, and X3 represent daylength, VPD, and GDD, respectively) is 0 (X ≤ Xmin), 1 (X ≥ Xmax), and (X − Xmin)/(Xmax − Xmin) otherwise. Xmax and Xmin are the maximum and minimum thresholds of each index, respectively. For the spring-GSI, Xmin was calculated as the minimum value among the values on the greenup day (from MCD12Q2.V006) for the study period of 2001–2018 at each study site, and Xmax was calculated as the minimum value among the values on the maturity day. For the fall-GSI, similarly, Xmin was the minimum value for the dormancy timings, and Xmax was the minimum value for the senescence timings. We incorporated GSI (or SGSI) into ED2 by multiplying it to the optimal value of leaf biomass on the day, where it operates as an upper limit of the leaf biomass.In this study, it was assumed that phenological stages are driven by meteorological conditions, not by other factors (e.g., no assumption of fixed phenological periods6,54,55). The development of a robust phenological model for the tundra ecosystem would be enabled by an increasing amount of ground-based phenology data (e.g., PhenoCam data37).Case study: flux data analysisThere were three sites where flux data is available for >5 years in Alaska; US-Atq site (flux data during 2004–2008 with delayed snowmelt in 2005), US-EML site (flux data during 2009–2017 with delayed snowmelt in 2017), and US-BZF (flux data during 2012–2018 with delayed snowmelt in 2017 and 2018). We first calculated two timings that are related to the meteorological conditions (0.1-GSI timing and half-max GSI timing51, Fig. S6) using the NARR data. The 0.1-GSI timing and half-max GSI timing were calculated on the day when the GSI passed 0.1 and the half-max value (i.e., 0.5), respectively, each year. To calculate two timings regarding the flux seasonal profile51 (source-sink transition timing and half-max productivity timing, Fig. S6), we used 30-min gap-filled FLUXNET201556 data (NEE and GPP; quality flags of measured or good) at the US-Atq site to calculate daily NEP (i.e., negative NEE) and daily GPP. At the US-EML and the US-BZF sites, we applied an open-source code called ONEFlux (Open Network-Enabled Flux processing pipeline for eddy-covariance data)57 using the ERA5 data (European Centre for Medium-Range Weather Forecast Reanalysis v558) which was downscaled with a quantile mapping method59. Using the gap-filled 30-min NEP and GPP data from the ONEFlux, we calculated the corresponding daily values and fitted a smoothing spline to the daily NEP and the daily GPP each year. The source-sink transition timing was defined as the day when the smoothing spline of the daily NEP passed zero in each year. The half-max productivity timing was set to the day when the smoothing spline of the daily GPP passed the half-max value in that year51.Further, we investigated whether the delayed snowmelt altered the relationships between meteorological conditions and the flux-threshold timings at each site based on (1) the correlation between the 0.1-GSI timing and the source-sink transition timing and (2) the correlation between the half-max GSI timing and the half-max productivity timing.ED2: model implementationWe used NARR data41 (0.3-degree resolution every 3 h; temperature, precipitation rate, pressure, v- and u-wind speed, downward longwave and shortwave radiation flux, and relative and specific humidity) for single-point ED2 implementation at each study site from 2001 to 2018. Vegetation structure (LAI, leaf and structural biomass, diameter at breast height, and population density) was initialized for each site by using the maximum annual LAI of cold-adapted shrubs and Arctic C3 grass from the Ent Global Vegetation Structure Dataset v1.0b (Ent-GVSD v1.0b) with the allometric equations in ED2. Ent-GVSD v1.0b provides plant functional types (from the MODIS land cover, MCD12C1.V00560) and maximum annual LAI values (from the MODIS LAI, MOD15A2.V00461,62) in subgrid cover fractions. We did not use canopy heights from Ent-GVSD v1.0b because of the absence of trees at the study sites. Soil texture (the ratio of sand:silt:clay) was set following the Harmonized World Soil Database v1.163 of the Food and Agriculture Organization of the United Nations (UN FAO). Soil carbon was initialized using the UN FAO Global Soil Organic Carbon Map64, and soil nitrogen was estimated using the soil C/N ratio of moist tundra (mean: 18.4)65.The prior distribution of each key variable was based on previous studies (Table S4), and 10,000 parameter sets were randomly generated from the prior distributions (the so-called Monte Carlo method). The best parameter set was selected based on statistical measures (r2 and root-mean-squared error) when compared to the data at the US-Atq site, i.e., NEP flux data for 2004–2006 and MODIS LAI data for 2003–2010 (MCD15A3H.V00666 at a 500-m resolution every 4 days) (Table S5). We then validated the performance of ED2 with this best parameter set by focusing on key ecosystem processes, such as NEP, ecosystem respiration, soil temperature, snowmelt timing, greenup timing, and the LAI at all sites (Table S5). The ED2 LAI was overestimated by 0.15–0.16 compared to the field-measured LAI values (Jul.–Aug. in 2006 at Barrow67 and Jun.–Aug. in 1996 at Toolik68).It is worth noting that the accuracy of the MODIS LAI has not been extensively evaluated at high latitudes because of limited ground measurements and few valid MODIS data points due to inadequate sun-sensor geometry, illumination conditions, and cloud contamination69,70. Furthermore, the heterogeneous landscapes of the region at the scale of remote sensing data (from hundreds of meters to a few kilometers) are also a major challenge that must be addressed before the data can be evaluated against ground measurements. According to the spatial representativeness assessment (see Supplementary Note, and Figs. S7 and S8), the landscapes around the flux towers generally have heterogeneity at a level similar to, or smaller than, the tower footprint size (200–300 m) during the early GS and peak GS in the MODIS 1 × 1 pixel window (i.e., 500 × 500 m2), but mostly higher than in the 3 × 3 pixel window (Table S6). This indicates that it is desirable to evaluate the MODIS 1 × 1 pixel LAI values against ground measurements, as both MODIS greenup and snowmelt timings agreed more with the ground data at the 1 × 1 pixel window scale than at the 3 × 3 pixel window scale (Figs. S9 and S10). A more thorough evaluation of both MODIS LAI data and ED2 LAI values is required in the coming years with the increase in ground data availability (e.g., National Ecological Observatory Network, NEON, LAI measurements).Correlation analysis: the effects of early or delayed snowmelt timingTo analyze the net and lagged effects of early or late snowmelt timing, it is necessary to constrain the contribution of interannual meteorological variation. Therefore, we compared only the years when meteorological conditions were similar, i.e., when the weekly mean GSI value was within one standard deviation of the weekly mean GSI during 2001–2018 (at the US-Hva site, weekly values during 1994–2018); meteorological conditions appeared similar for 8 or 9 years at each study site, except the US-BZF site, where the similarities were found for 10 years. We also limited the effect of greenup timing changes by excluding the years when greenup was earlier or later by one standard deviation of the mean of greenup timings during the study period. For the years satisfying the constraints, a least-squares linear regression was applied between the snowmelt timing change (deviation in snowmelt timing each year from the mean snowmelt timing) and the seasonal deviation (the difference of the seasonal mean from the mean value over the years) of each process from the ED2 results.To analyze the net and lagged effects of delayed snowmelt, we implemented the ED2 model in two schemes, (1) following the meteorologically-determined phenological index (i.e., the GSI, Eqs. (2) and (2) constraining leaf-out by snowmelt (i.e., the SGSI, Eq. (3)). For the years when unmelted snow delayed greenup, we took the difference between the modeled results (i.e., GSI and SGSI) for each process and applied a least-squares linear regression between the difference of each process and the delayed snowmelt days. More

Insect rearingA CLas negative colony of ACP was initially collected from CLas-free Murraya exotica L. growing in the ornamental landscape of South China Agricultural University (SCAU, Guangzhou, China) in May 2014. Then it was reared on potted M. exotica in a greenhouse at SCAU. M. exotica plants were pruned regularly to promote the growth flushes necessary to stimulate ACP oviposition. The ACP populations were periodically (at least once a month) tested to ensure the colony was CLas-free using nested quantitative PCR detection according to the method described by Coy et al.30.The parasitoid T. radiata used in the current study was initially collected from ACP hosts on M. exotica plants in the above-mentioned location during June 2015. Its population was maintained in rearing cages (60 × 60 × 60 cm) using a CLas-free ACP-M. exotica rearing system under laboratory conditions (26 ± 1 °C, RH 80 ± 10% with L:D = 14:10 photoperiods in insect incubators).Host plantsCLas-free and CLas-infected plants of Citrus reticulata Blanco cv. Shatangju were used in the current study. Both plant types were obtained from The Citrus Research Institute of Zhaoqing University (Guangdong, China). The CLas-infected plants were inoculated by shoot grafting. All plants were approximately 4-year old and 1.2−1.5 m in height, separated in nylon net greenhouses (70 mesh per inch2) at two different locations about 2.2 km apart in SCAU. Again, nested qPCR detection was performed periodically (at least once a month) to test for the presence or absence of CLas in the citrus plants according to the method described by Coy et al.30.Acquisition and persistence of CLas in Tamarixia radiata
When new shoots of CLas-infected C. reticulata plants were grown to 5–8 cm, 20 pairs of 1 week-old ACP adults were introduced into one nylon bag covering one fresh shoot to lay eggs for 48 h. When the progeny of ACP developed through to 4th or 5th instar nymph (CLas-donor ACP), which are the stages preferred by T. radiata parasitoids, 150 of the ACP nymphs were randomly selected and the remaining ones were removed. Following this, 10 pairs of 3-day old T. radiata adults, randomly selected from the population that has been tested to be CLas-free, were introduced into the nylon bag in order to parasitize the 4th or 5th instar ACP nymphs for 48 h before being recaptured. Then the potentially parasitized ACP nymphs together with the citrus plants were cultured in a plant growth chamber (Jiangnan Instrument Company, RXZ-500D, at 26 ± 1 °C, 60 ± 2% RH and 14:10 h L:D photoperiod of 3,000 lx illumination).When the progeny of T. radiata (considered F0 generation) developed to 3-day egg, 1st to 4th instar larvae, pupae, and adult stages respectively, they were identified and collected with the assistance of a stereomicroscope. DNA of each stage sample was extracted using the TIANamp Genomic DNA Kit (TIANGEN, Beijing, China) for CLas qPCR detection and titer quantification. Thirty eggs, 20 individuals of 1st or 2nd instar, 10 individuals of 3rd or 4th instar larvae or pupa, as well as three individuals of female or male adults were subsequently ground together to represent each life stage in qPCR, and each stage qPCR detection was repeated three times.The primers used for CLas qPCR detection were LJ900 primers, (F5′-GCCGTTTTAAC ACAAAAGATGAATATC-3′, R5′-ATAAATCAATTTGTTCTAGTTTAC GAC-3′), and 18S rRNA gene of T. radiata (F5′-AAACGGCTACCACATCCA-3′, R5′-ACCAGACT TGCCCTC CA-3′)31 was used as an internal control for DNA normalization and quantification. In order to normalize the qPCR values, each qPCR reaction was performed in three independent runs using SYBR Premix Ex Taq (Takara, Dalian, China) in Bio-Rad CFX Connect™ Real-Time PCR Detection System, with a protocol of initial denaturation at 95 °C for 3 min, followed by 40 cycles at 95 °C for 10 s, 60 °C for 20 s and 72 °C for 30 s.To monitor the CLas persistence in T. radiata, newly emerged female adults of T. radiata (considered F1 generation) were collected from the above experiment and fed with 20% honey water. After 1, 5, 10, and 15 days, 10 parasitoids were recaptured, subsequently ground for DNA extraction and CLas titer detection and quantification using qPCR. The protocol of DNA extraction and qPCR reaction was the same as above, and qPCR quantification was repeated three times for each treatment.Localization patterns of CLas in Tamarixia radiata
Localization patterns of CLas in different instars of T. radiataFluorescent in situ hybridization (FISH) was used to visualize the distribution of CLas in T. radiata exposed to CLas positive ACP, following the method of Gottlieb et al.32 with a slight modification. Eggs and different larval instars of T. radiata were collected and fixed in Carnoy’s solution (chloroform-ethanol-glacial acetic acid [6:3:1,vol/vol] formamide) overnight at 4 °C. After fixation, the samples were washed three times in 50% ethanol with 1× phosphate buffered saline (PBS) for 5 min. Then the samples were decolorized in 6% H2O2 in ethanol for 12 h, after which they were hybridized overnight in 1 ml hybridization buffer (20 mM Tris-HCl pH 8.0, 0.9 M NaCl, 0.01% sodium dodecyl sulfate, 30% formamide) containing 10 pmol of fluorescent probes/ml in a 37 °C water bath under dark conditions. The CLas probe used for FISH was 5′-Cy3-GCCTCGCGACTTCGCAACCCAT-3′. Finally, the stained T. radiata samples were washed three times in a washing buffer (0.3 M NaCl, 0.03 M sodium citrate, 0.01% sodium dodecyl sulfate, 10 min per time). After the samples were whole mounted and stained, the slides were observed and photographed using a Nikon eclipse Ti-U inverted microscope. For each stage sample, approximately 20 individuals were examined to confirm the results.Localization patterns of CLas in different organs of T. radiataDifferent organs (gut, fat body, ovary, poison sac, salivary glands, spermatheca, and chest muscle) were dissected from newly emerged adults of T. radiata in 1× phosphate buffered saline (PBS) under a stereomicroscope using a depression microscope slide and a fine anatomical needle. After a sufficient number of each tissue sample was collected (20 or more), the tissues were washed three times with 1 × PBS, followed by the fixation, decolorization, and hybridization procedures as outlined above, except that this time of decolorization was 2 h. After hybridization, nuclei in the different organs were counterstained with DAPI (0.1 mg/ml in 1 × PBS) for 10 min, then the samples were transferred to slides, mounted whole in hybridization buffer, and viewed using confocal microscopy (Nikon, Japan).Maternal transmission of CLas between Tamarixia generationsFive groups of experiments were used to clarify whether CLas can be transmitted vertically between different T. radiata generations. In the first group, 60 pairs of newly emerged T. radiata adults from the CLas-infected ACP colony (potential CLas-acquired parasitoid adults, F0 generation) were introduced into 60 nylon bags (one female per cage). Each bag covered one fresh citrus plant shoot with one marked CLas-free 4th instar nymph of ACP, the parasitoid females were given 24 h to oviposit, then transferred to another four groups successively to oviposit with intervals of 24 h before they were recaptured for CLas-PCR detection (58/60 and 56/60 T. radiata females and males respectively were CLas-infected). Only the progeny (F1 generation) in which parasitoid parents were both CLas-infected continued to be investigated.When the F1 progeny of CLas-infected parasitoid females developed to egg, larval, pupal, and adult stages respectively, they were collected and divided into two groups; in one group samples were used for the qPCR detection of the CLas titer, and the other group was used for the FISH visualization of CLas. The qPCR and FISH analysis protocols of CLas as well as the number of tested individuals were the same as previously outlined. Each stage was repeated three times.
CLas detection in T. radiata-inoculated ACPQuantitative PCR detection of CLasApproximately 60 newly emerged parasitoid adult females from CLas-infected ACP hosts (potential CLas-acquired parasitoid adults) were collected using an aspirator. They were first starved for 5 h, then released into finger tubes (diameter 6 mm × length 30 mm); one female per tube containing one 4th instar nymph of CLas-free ACP (this was treated as one experimental replicate). The probing behavior of the parasitoids was observed under a stereomicroscope, after which the parasitoids were recaptured for CLas PCR detection (similar to the above experiment, approximately 95% were CLas-infected). Only those 4th instar ACP nymphs, probed for egg-laying by a CLas-infected parasitoid but survived from the probing (the averaged proportion of such samples was 5.36 ± 0.47% and were 100% CLas infected), were transferred onto fresh CLas-free M. exotica shoots to complete their development (hereafter referred as “T. radiata-inoculated ACP”). The experiment was repeated in 32 parallel replicates (Supplementary Table 1), in which 103 T. radiata-inoculated ACP nymphs were finally obtained.Following the above, thirty T. radiata-inoculated ACP nymphs were collected when they developed into 5th instar nymphs (the stage when infection proliferation might have just begun since the infection was introduced at the 4th instar). In addition, thirty 8-day old adults that developed from the T. radiata-inoculated ACP nymphs were also collected. This was because the results in Wu et al.28 revealed that the proportion of CLas-infected ACP individuals exceeds 90% at the 12th day after infection acquisition, while ACP takes 4 days to develop into an adult from 5th instar stage. Their alimentary canals and salivary glands were dissected under a stereomicroscope using the methods of Ammar et al.33, and hemolymphs were collected with a 10 μl pipette tip using the method of Killiny et al.34. The DNA of the alimentary canals, salivary glands and hemolymphs were extracted using TIANamp Micro DNA Kit (Tiangen, Beijing, China), and the relative titers of CLas in each tissue of ACP nymphs and adults were detected by qPCR with of LJ900. The β-actin gene of ACP (F 5′-CCCTGGACTTTGAACAGGAA-3′; R 5′-CTCGTGGATACCGCAAGATT-3′) was selected as an internal control for data normalization and quantification35. For each sample, qPCR detection was repeated three times.FISH visualization of CLasThe alimentary canals and salivary glands of 5th instar nymphs and 8-day old adults of T. radiata-inoculated ACP were dissected as described above, and the distribution of CLas was visualized by FISH and confocal microscopy. The alimentary canals and salivary glands of CLas-infected ACP nymphs and adults (collected from CLas-infected citrus plants) were used as a positive control, and five to ten samples were detected by FISH for each tissue.
CLas transmission from T. radiata-inoculated ACP to citrus plantsAccording to the above experimental results, if the CLas could be detected in the salivary glands of the 8-day old ACP adults (T. radiata-inoculated ACP), 30 more of these adults were randomly selected to inoculate on fresh shoots of CLas-free citrus. ACP adults that acquired CLas from plants and CLas-free ACP adults were used as positive and negative controls respectively.After 20, 30, 40, and 50 days of feeding samples of the citrus leaves fed on by T. radiata-inoculated ACP (named as CLas-recipient citrus leaves), fed on by ACP that acquired CLas from plants (positive control), and fed on by CLas-free ACP (negative control) were cut (1 cm2). Their DNAs were extracted using DNAsecure Plant Kit (Tiangen, Beijing, China). The infections of CLas in these plants were detected by nested PCR based on the methods of Jagoueix et al.36 and Deng et al.37. The experiment was repeated in six plants for each of 20, 30, 40, and 50 days feeding duration, and the infection rates of CLas were calculated.Localization of CLas in citrus plants fed on by T. radiata-inoculated ACPIn order to further confirm the infection of CLas in the recipient citrus leaves, FISH was used to visualize the localization of CLas. According to the above experimental results, after being fed on for 50 days by the T. radiata-inoculated ACP adults, citrus leaf sections containing the midrib were cross-sliced in 30 µ sections using a cryostat (CM1950, Leica, Germany). The leaf samples were prepared for FISH vitalization according to the protocol described by Gottlieb et al.32. Citrus leaves from the plant that had been fed on by ACP adults that acquired CLas from plants and CLas-free ACP adults were used as positive and negative controls, respectively. Five to 10 leaf samples were detected by FISH for each treatment.Phylogenetic analysis of CLas bacteria in different ACP populations and citrus plantsTo assess the identity of the CLas bacteria in CLas donor ACP, CLas vectored parasitoids, T. radiata-inoculated ACP and recipient citrus leaves, the outer membrane protein gene (omp) of CLas was PCR amplified with the primers HP1asinv (5′-GATGATAGG TGCATAAAAGTACAGAAG-3′) and Lp1c (5′-AATACCCTTATGGGATACAAAAA-3′) following the procedure described in Bastianel et al.38. Then the PCR products were sent for sequencing after visualizing the expected bands on 1% agarose gels.All the DNA sequences of CLas omp gene were edited and aligned manually using Clustal X1.8339 in Mega 640. The best model and partitioning scheme were chosen using the Bayesian information criterion in PartitionFinder v.1.0.141. Phylogenetic analysis was undertaken using a maximum likelihood (ML) method with 1000 non-parametric bootstrap replications in RAxML42. Escherichia coli was used as an outgroup.Statistics and reproducibilityTaking 18S rRNA gene of T. radiata and the β-actin gene of ACP as housekeeping genes, the relative titers of CLas in different stages and different tissues of T. radiata and ACP were calculated using the method of 2[−ΔΔct 43. For the parallel experiments that had more than three replicates the differences were compared using analysis of variance (ANOVA) with SPSS 18.0 at a significance level α = 0.05; while for CLas titer, two-sample comparison between genders of Tamarixia adults analysis was performed using paired t-test. Fluorescent pictures were processed using Photoshop CS5 software.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Insect extractsThai edible insects (Fig. 1) were extracted and yield of each extract is shown in Fig. 2. Hexane extracts of most insects, except for P. succincta, provided the highest yield, followed by ethanolic extracts, and aqueous extracts, respectively. The reason might be due to a high amount of fat content of insects. Since these fat components are hydrophobic, they could be extracted well using nonpolar solvent, e.g. hexane. Semi-polar solvent like ethanol could also be used to extract hydrophobic compounds but with less extraction efficacy5. Several previous studies reported that fat was abundant in biomass of insects, ranging from 4.2 to 77.2%, which was accounted for about 26.8% on average dried insects6,7.Figure 1External appearances of Thai edible insects, including (a) rice grasshopper (Euconocephalus sp.), (b) bamboo caterpillar (O. fuscidentalis), (c) house cricket (A. domesticus), (d) silkworm pupae (B. mori), (e) Bombay locust (P. succincta), and (f) giant water bug (L. indicus).Full size imageFigure 2Yields of insect extracts, including B. mori (BM), O. fuscidentalis (OF), Euconocephalus sp. (EU), P. succincta (PS), A. domesticus (AD), and L. indicus (LI). The data are expressed as mean ± SD (n = 3). The Greek alphabet letters (α, β, γ, and δ) indicate significant differences among hexane extracts, the capital letters (A, B, C, and D) indicate significant differences among ethanolic extracts, and the small case letters (a, b, and c) indicate significant differences among aqueous extracts. The data were analyzed using One-Way ANOVA followed by post hoc Tukey test (p  More