Ecology

Timothy Searchinger and his colleagues raise concerns that the European Union’s plan to produce energy from biomass could compromise forest carbon stocks and biodiversity (Nature 612, 27–30; 2022). However, it is possible for improved forest management to reconcile increased bioenergy production by maintaining and restoring forest ecosystems.
Competing Interests
The authors declare no competing interests. More

Data collectionWe obtained trade data from two different sources: the United States Fish and Wildlife Service (USFWS) Law Enforcement Management Information System (LEMIS)31 and the International Union for Conservation of Nature (IUCN) Red List32. We used the former to obtain data on the live wildlife trade in general and the latter for data on the pet trade specifically. We then matched trade data with our previously compiled global scale datasets of life history traits and introductions in mammals, reptiles and amphibians25,26.We obtained data on the US live wildlife trade from LEMIS by a Freedom of Information Act Request on 12/08/2019. We requested summary data on all US imports and exports of wildlife across all available years (1999-2019) and all trade purposes, including information on species identities and shipment contents (e.g. live individuals, meat, skins, etc.). For each species, we summed the total number of recorded shipments of live individuals (including individuals that died in transit, and live eggs) as a measure of trade frequency. We classified species as in trade if there was at least one shipment of live individuals recorded in the LEMIS database, and as not traded otherwise. The LEMIS dataset is geographically limited to trade by the US, and therefore may not capture the full diversity of species involved in the wildlife trade. For example, the LEMIS database may be missing some species involved in the substantial trade in live wildlife between South–East Asian countries50. However, the US represents one of the most dominant players in the global market for live wildlife16, and by summing both imports and exports we capture demand for species in countries beyond the US to some extent. Supplementary Fig. 2 illustrates the frequency of trade between the US and countries represented in the US LEMIS dataset. LEMIS data should be considered a minimum estimate of the diversity of species involved in the wildlife trade since they mostly record only legal trade (although confiscated shipments are recorded), and shipments are sometimes not identified to the species level16,51,53,53. The LEMIS database also contains some mis-spelled and incorrectly identified species due to human input errors52. To minimise the effect of misidentified shipments on our species level classifications of US trade status, we discarded all LEMIS records that were not identified to the species level (i.e. those identified using genus, common or generic names only), and manually checked the LEMIS data for synonyms and alternate spellings when we could not automatically match any records in LEMIS with species in our life history datasets. Species classified as traded on the basis of a single recorded live shipment in LEMIS are most vulnerable to species level misclassification due to misidentified shipments. The vast majority of traded species have multiple shipments recorded in LEMIS (259/312 [83%] of traded mammals, 265/285 [93%] of traded reptiles and 72/75 [96%] of traded amphibians), reducing the potential impact of shipment level misidentification over the reliability of species level trade classifications. However, to investigate the robustness of our findings to possible errors in species identification in LEMIS, we re-ran our key analyses excluding species classified as traded on the basis of a single live shipment. We found qualitatively the same effects of life history traits on the probability of trade when removing these species as in our full sample (Supplementary Tables 25–27). Despite its limitations, LEMIS is an invaluable resource for identifying broad scale trends in the wildlife trade since few other countries maintain such detailed records, and it is the only large-scale international trade dataset that includes both CITES- and non-CITES-listed species16,41. Including non-CITES listed species in our analyses is important because CITES-listed species represent only a small minority of those in trade14 and are likely to be a biased sample in terms of life history traits, since species vulnerable to extinction typically have slower life histories40.We obtained separate data on the pet trade from the IUCN Red List. The IUCN has assessed the vast majority of mammal, reptile and amphibian species (91%, 79% and 86% respectively54). Here, we classified a species as involved in the pet trade if the IUCN species account included at least one clear description of involvement in the pet trade. Otherwise, we considered a species as not involved in the pet trade. Although LEMIS records the purpose of trade, it uses broad categories (e.g. ‘Commercial’, ‘Personal’, ‘Breeding in captivity’), none of which refers specifically to nor necessarily equates to trade for pets. Therefore, we sought this additional data on the pet trade from the IUCN Red List instead of following the approach of some previous studies which have used LEMIS data as a proxy for the pet trade (e.g. Refs. 15,19). In contrast, the IUCN Red List contains clear textual descriptions of use and trade for many species, allowing us to identify which species are traded specifically for pets32. The IUCN data has further complementary strengths compared with LEMIS in that it is global in scope and includes both legal and illegal trade. We obtained data from the IUCN Red List by manually searching the binomial name of each species in our samples and consulting the ‘Threats’ and ‘Use and Trade’ sections of the species accounts. We classified species as in the pet trade if the information clearly stated this was the case (e.g. “It has been recorded in the pet trade”, “This species appears in the international pet trade”). We discounted descriptions where the information was uncertain (e.g. the species is described as “probably” or “possibly” traded for pets). We did not count as pets those species that the IUCN categorises as used for “Pets/display animals, horticulture” but which are used only for zoos or captive display, such as beluga whales (Delphinapterus leucas). All species described as pets by the IUCN are ‘exotic’, i.e. those without a long history of domestication14, since the IUCN does not list domesticated species.We matched trade data with our previously published global scale datasets on life history traits and introductions25,26. Internationally traded species may or not be released in the wild outside their native range: some may remain in the confines of captivity (e.g. in zoos or kept by private owners). We defined a species as introduced if there was at least one reliable record of its release, by humans, into the wild outside of its native range, either accidentally or intentionally25,26. We included only species with complete data for the same life history traits as used in our prior analyses (mammals: body mass, gestation period, weaning age, neonatal body mass, litter size, litters per year, age at first reproduction and reproductive lifespan; reptiles: body mass, hatchling mass, clutch size, clutches per year, age of sexual maturity, reproductive lifespan and parity; amphibians: snout-vent length, egg size, clutch size, age of sexual maturity and reproductive lifespan) to facilitate direct comparisons with previous results and to allow us to account for covariation between life history traits55. Species with complete life history data represent 7.8%, 3.5% and 1.6% of the total estimated number of species of mammals, reptiles and amphibians respectively56,57,58. These samples are not random as they over-represent orders containing many species of interest and utility to humans (e.g. ungulates, primates, crocodilians) (Supplementary Tables 28–30). However, these biases are unlikely to undermine our results since we examine life history effects on trade and introduction within these samples. Trade and introduction data do not necessarily cover the same time periods: the US dataset covers only the years 1999-present and the IUCN descriptions also typically refer to recent trade. In contrast, our introduction dataset includes both historical and recent introductions25,26. Therefore, the goal of our analyses is not to test causal hypotheses on the direct relationship between trade and introduction but rather to investigate whether the same life history traits predispose species towards both trade and introduction across diverse taxa, locations and circumstances. When combining the datasets and phylogenies59,60,61,62,63, we resolved species name mis-matches by referring to taxonomic information from the IUCN Red List32, the Global Biodiversity Information Facility (GBIF33) and the Integrated Taxonomic Information System (ITIS64). Table 1 summarises final sample sizes and Supplementary Table 1 the degree of overlap between the trade datasets. Most species in the pet trade are also in the general live wildlife trade, but many more species are traded by the US for general purposes than are involved in the pet trade specifically.Finally, we obtained data for a proxy measure of species detectability in order to control for a potential confounding effect on relationships between life history traits and introduction: larger bodied and longer-lived species may be more likely to be recorded by human observers when introduced compared with smaller and shorter-lived species. We obtained data on species occurrence records, geographic range size and population density, assuming that highly detectable species will have a disproportionately large number of recorded observations than expected based on the size of their geographic ranges and average population densities, following similar approaches by e.g. Refs. 65,66. We obtained occurrence records from the Global Biodiversity Information Facility (GBIF33) via the R package rgbif67 selecting only records resulting from human observation. We obtained range sizes (in decimal degrees squared) from the IUCN Red List32 and processed them for analysis using functions from the rgdal package68, excluding areas of uncertain presence (i.e. limiting range to presence code 1, ‘extant’). We obtained population density estimates from the TetraDENSITY database (version 134), a global database of population density estimates for terrestrial vertebrates. Most species in the TetraDENSITY dataset are represented by estimates from multiple different studies (median = 3, range 1–408). We collapsed density estimates to the species level by taking the median value across studies, including all estimates regardless of sampling method to maximise sample size, and converting all units to individuals/km2 to ensure comparability.Statistical analysesTo investigate relationships between life history traits and trade, we run models treating US or pet trade as the outcome variable and life history traits as the predictors. For all analyses, all life history variables were included in the same models to account for covariation among life history traits55. For US trade, where data on trade frequency are available, we run models both in which trade is treated as a binary variable (traded vs. not traded) and as a count variable (frequency of live shipments, including zero values), while for the pet trade, we have no data on trade frequency and so we treat pet trade as a binary variable only. To investigate the effects of life history traits on introduction, we run models in which introduction is the outcome variable and life history traits are the predictors. In introduction models, we only include traded species (running separate models for the set of species in US trade and the set of species in the pet trade). This approach allows us to disentangle effects associated with trade and introduction and thus identify at which stage(s) life history biases emerge. We also run introduction models in which frequency of US trade is included as an additional predictor alongside life history traits, anticipating that highly traded species are more likely to be introduced. Finally, to investigate possible confounding effects of species detectability on relationships between life history traits and introduction, we investigate effects of number of observations, geographic range size and, where sample sizes allowed, population density on the probability of introduction. If highly detectable species are more likely to be recorded as introduced, we expect to find a positive effect of the number of observations (while accounting for geographic range size and population density) on the probability of introduction. If this effect confounds relationships between body mass/lifespan and introduction, the effect of these life history traits on the probability of introduction should disappear when detectability measures are included in the models alongside life history traits. All analyses were conducted using the R statistical programming environment (Version 4.2.069). Plots were coloured using palettes from the viridis package70.To estimate effects of predictor variables, we fit generalized linear mixed models (GLMMs) using Markov chain Monte-Carlo (MCMC) estimation, implemented in the MCMCglmm package35,36. For analyses with binary outcome variables (traded vs. not traded, introduced vs. not introduced) we fit probit models, while for analyses with US trade frequency as the outcome variable we fit hurdle models. Hurdle models estimate two latent variables: the probability that the outcome is zero (on the logit scale), and the probability of the outcome modelled as a Poisson distribution for non-zero values71. This method therefore allows us to estimate effects of life history traits on the probability and frequency of trade in the same model. While the binary component of a hurdle model estimates the probability that outcomes are zero, when reporting results we reverse the sign of coefficients from the binary model for ease of interpretation, so that effects can be interpreted as the probability that the outcome is not zero. Therefore, here predictors with consistent effects on the probability and frequency of trade in hurdle models will have the same sign (so that if, for example, litter size has a positive effect on both the probability and frequency of trade, both coefficients for litter size from the hurdle model will be positive).Datasets comprising biological measures from multiple related species violate the fundamental statistical assumption that observations are independent of one another, since closely related species are more phenotypically similar than expected by chance due to their shared evolutionary history72. To account for the non-independence of species due to shared ancestry, we included a phylogenetic random effect in all models, represented by a variance-covariance (VCV) matrix derived from the phylogeny. The off-diagonal elements of the VCV matrix contain the amount of shared evolutionary history for each pair of species35,37,38 based on the branch lengths of the phylogeny (here proportional to time)59,61,62,63,63. This approach allows us to estimate phylogenetic signal using the heritability (H2) parameter, which measures the proportion of total variance in the latent variable attributable to the phylogeny35,37,38. Heritability is interpreted in the same way as Pagel’s λ in phylogenetic generalized least squares regression35,37,38,72. Specifically, phylogenetic signal is constrained between 0, indicating no phylogenetic effect so that species can be treated as independent, and 1, indicating that similarity between species is directly proportional to their amount of shared evolutionary history35,38,72. As hurdle models estimate two latent variables, for each hurdle model we report two heritability estimates, one for the binary and one for the Poisson component. All continuous independent variables were log-10 transformed due to positively skewed distributions. Although GLMMs do not require normally distributed predictor variables, log-transforming positively skewed life history predictors in phylogenetic comparative analyses allows us to model life history evolution on proportional rather than absolute scales. This is important as it facilitates biologically meaningful comparisons between species across large scales of life history variation73. Further, log-transforming positively skewed predictors helps to meet assumptions of the underlying Brownian motion model of evolutionary change, which assumes that phenotypic change along branches of the phylogeny is normally distributed74.We calculated variance inflation factors (VIFs) using functions from the car R package75 to check for multicollinearity between predictor variables. Where any model reported a variance inflation factor of 5 or above, indicating potentially problematic levels of collinearity76, we re-ran the model removing the variable with the highest VIF iteratively until all the remaining variables had VIFs of More

Study siteThe field experiment was conducted at Namco Station (30°47’N, 90°58’E, altitude 4730 m) of the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS), which is located in the alpine steppes of TP in China. The experiment was permitted by ITPCAS, complied with local and national guidelines and regulations. From 2006 to 2017, the mean annual temperature (MAT) and mean annual precipitation (MAP) was about − 0.6 °C and 406 mm, respectively. Monthly mean temperature varied from − 10.8 °C in January to 9.1 °C in July and most of the precipitation occurred from May to October37,38. During our six-year observations (2010, 2011, 2012, 2013, 2015 and 2017), climate change during the growing season from May to September varied differently, with the annual precipitation ranged from 255.9 mm to 493.8 mm and the MAT from 6.7 to 7.4 °C. Androsace tapete, Kobresia pygmaea, Stipa purpurea and Leontopodium pusillum were the dominant plant species at the alpine steppe.Experimental design and treatmentsThe long-term experiment began in May, 2010. Three homogenous plots were randomly arranged as replicates at the alpine steppe and six subplots (~ 13 m2) were distributed in each plot by a cycle, with a 2 m buffer zone between each adjacent subplot (Appendix S1: Fig. S1). In this experiment, six treatments of N fertilization rate (0, 1, 2, 4, 8, and 16 g N m−2 yr−1) were clockwise applied in each subplot. The subplots of 0 g N m−2 yr−1 were control group. We sprayed NH4NO3 solution on the first day of each month in the growing season (from May to September) each year. After fertilizing, we rinsed plant residual fertilizer with a little deionized water (no more than 2 mm rainfall). For the control groups, we added equivalent amount of water. The experiment was conducted from 2010 to 2017 (it should be pointed out that there was no fertilization in 2014 and 2016).Sampling and measurementsThe samples were collected with the training and permission of ITPCAS and involved plants that are common species and not endangered or protected. The identification of the plants was done by referring to a book of Chen and Yang39. Pictures of the corresponding specimens can be seen on the website of ITPCAS (http://itpcas.cas.cn/kxcb/kxtp/nmc_normal_plant/).Vegetation samples were collected in August in 2011 and repeated at the same time in 2012, 2013, 2015 and 2017. We established one 50 × 50 cm quadrat in each subplot, clipped aboveground biomass (AGB) and sorted species by families. The biomass was used to measure ANPP (g m−2 yr−1). Following aboveground portion collected, we used three soil cores (5 cm diameter) to collect the belowground roots at 0–30 cm depth and mixed into one sample, which were used to assess belowground net primary productivity (BNPP, g m−2 yr−1). The roots were cleaned with running water to remove sand and stones.Both plant and root samples were dried at 75 °C for 48 h and then ground into powder (particle size ~ 5 μm) by a laboratory mixer mill (MM400, Retsch). To determine N and C content of plants, we weighed the samples into tin capsules and measured with the elemental analyzer (MAT253, Finnigan MAT GmbH, Germany).Estimation of the critical N rate (Ncr), N retention fraction (NRF), N retention capacity and N-induced C gainAccording to the N saturation hypothesis, plant productivity increases gradually during N addition, reaches a maximum at the Ncr, and eventually declines16,17. We considered the Ncr to be the rate where ANPP no longer remarkably changed with N addition (Fig. 1).We defined plant N retention fraction (NRF, %; Eq. 1) as the aboveground N storage caused by unit N addition rate, and N retention capacity (g N m−2 yr−1; Eq. 2) was the increment of N storage due to exogenous N addition compared to the control40. The equations are as following:$$N;retention;fraction = frac{{ANPP_{tr} times N;content_{tr} – ANPP_{ck} times N;content_{ck} }}{N;rate}$$
(1)
$$N;retention;capacity = ANPP_{tr} times N;content_{tr} – ANPP_{ck} times N;content_{ck}$$
(2)
where ANPPtr and N contenttr (%) refer to those in the treatment (tr) groups, and ANPPck and N contentck refer to those in the control (ck) groups. These expressions are also used in the following equations (Eqs. 3–5).The N-induced C gain (g C m−2 yr−1; Eq. 3) was estimated by the increment of C storage owing to exogenous N addition compared to the control40. Maximum N retention capacity (MNRC, Eq. 4) and maximum N-induced C gain (Eq. 5) mean the maximum N and C storage increment in plant caused by exogenous N input at Ncr, respectively. The formulas are as following:$$N{text{-}}induced;C;gain = ANPP_{tr} times C;content_{tr} – ANPP_{ck} times C;content_{ck}$$
(3)
$$MNRC = ANPP_{max } times N;content_{max } – ANPP_{ck} times N;content_{ck}$$
(4)
$$Maximum;N{text{-}}induced;C;gain = ANPP_{max } times C;content_{max } – ANPP_{ck} times C;content_{ck}$$
(5)
where ANPPmax, N contentmax and C contentmax refer to the value of ANPP, N content and C content at Ncr, respectively.Data synthesisTo evaluate N limitation and saturation on the TP more accurately, we searched papers from the Web of Science (https://www.webofscience.com) and the China National Knowledge Infrastructure (https://www.cnki.net). The keywords used by article searching were: (a) N addition, N deposition or N fertilization, (b) grassland, steppe or meadow. Article selection was based on the following conditions. First, the experimental site must be conducted in a grassland ecosystem. Second, the experiment had at least three N addition levels and a control group. Third, if the experiment lasted for many years, we analyzed data with multi-year average. Based on the above, we collected 89 independent experimental cases. Among these, 27 cases were located on the TP alpine grasslands, 25 in the Inner Mongolia (IM) grasslands and 37 in other terrestrial grasslands (detailed information sees Appendix S2: Table S1).We extracted ANPP data and N addition rate of these cases and estimated Ncr and ANPPmax (Appendix S2: Fig. S2). We then calculated NRF, N retention and C gain of each group of data for further analysis (Appendix S2: Table S2). Most of the 89 cases did not have data on N and C content. To facilitate the calculation, we summarized N and C content from 40 articles in the neighboring areas of the cases and divided the N and C content into seven intervals according to the N addition rate (Appendix S2: Table S3 and Fig. S3). The unit of N addition rate was unified to “g N m−2 yr−1”. All the original data were obtained directly from texts and tables of published papers. If the data were displayed only in graphs, Getdata 2.20 was used to digitize the numerical data. For the estimation of N retention and C gain of the TP at current N deposition rates and future at Ncr, we fitted the exponential relationship to the data from 27 cases on the TP, and then substituted N rates into the fitted equations (Eq. 6):$$y = a times left[ {1 – exp left( { – bx} right)} right].$$
(6)
We also included MAT, MAP, soil C:N ratio, fencing management (fencing or grazing) and grassland type (meadow, steppe and desert steppe) of the experiment sites for exploring the drivers affecting N limitation (Appendix S2: Table S1). When climatic data were missing from the article, MAT and MAP were obtained from the WorldClim (http://www.worldclim.org).Species were usually divided into four functional groups (grasses, sedges, legumes and forbs) to study the response of species composition to N addition in previous study41. We synthesized 13 TP experimental cases (including our field experiment) from the data synthesis and each case included at least three functional groups (detailed references see Appendix S2).Statistical analysisThere were 42 species in our field experiment. We divided them by family into eleven groups: Asteraceae (forbs), Poaceae (grasses), Leguminosae (legumes), Rosaceae (forbs), Boraginaceae (forbs), Caryophyllaceae (forbs), Cyperaceae (sedges), Labiatae (forbs), Primulaceae (forbs), Scrophulariaceae (forbs) and Others. Due to species in the group of Others contributed only 1.22% of AGB, we analyzed AGB and foliar stoichiometry among other ten families (Appendix S1: Table S1). In Namco steppe, forbs, grasses, sedges and legumes accounted for 78.0%, 7.4%, 8.2% and 5.2% of the AGB respectively (Appendix S1: Table S1 and Fig. S2). Such a large number of forbs suggested that our experiment was conducted on a severely degraded grassland.For our field data, two-way ANOVAs were used to analyze the effects of year, N fertilization rate and their interactions on species AGB. One-way ANOVAs were used to test the response of ANPP, BNPP, root:shoot ratio, species foliar C content, N content and C:N ratio to N addition rate. Duncan’s new multiple range test was used to compare the fertilization influences at each rate in these ANOVAs. Prior to the above ANOVAs, we performed homogeneity of variance test and transformed the data logarithmically when necessary. Simple regression was used to estimate the relevance among ANPP, NRF, N retention capacity and C gain with N addition rates.Structural equation modeling (SEM) was used to explore complex relationships among multiple variables. To quantify the contribution of drivers such as climate and soil to Ncr, ANPP, NRF and MNRC, we constructed SEM based on existing ecological knowledge and the possible relationships between variables. We considered environmental factors (MAT, MAP and soil C:N) and ANPPck as explanatory variables, and Ncr, NRF and MNRC as response variables. We included the ANPPck in the SEM rather than the ANPPmax because we wonder whether there was a relationship between ANPP in the absence of exogenous N input and the ecosystem N retention in the presence of N saturation. This has important implications for assessing N input. Before constructing the SEM, we excluded collinearity between the factors. In addition, Student’s t-test and one-way ANOVAs were performed to explain the effect of fencing management and grassland type on above response variables, respectively. The SEM was constructed using the R package “piecewiseSEM”42. Fisher’s C was used to assess the goodness-of-model fit, and AIC was for model comparison.Given the influence of extreme values in the data synthesis, we calculated the geometric mean of Ncr, NRF, N retention and N-induced C gain. All statistical analyses were performed with SPSS 26.0 and RStudio (Version 1.2.1335) based on R version 3.6.2 (R Core Team, 2019). More

Nadell, C. D., Drescher, K. & Foster, K. R. Spatial structure, cooperation and competition in biofilms. Nat. Rev. Microbiol. 14, 589–600 (2016).Article 
CAS 

Google Scholar 
Palmer, J. D. & Foster, K. R. Bacterial species rarely work together. Science 376, 581–582 (2022).Article 
ADS 
CAS 

Google Scholar 
Pande, S. & Kost, C. Bacterial unculturability and the formation of intercellular metabolic networks. Trends Microbiol. 25, 349–361 (2017).Article 
CAS 

Google Scholar 
Nadell, C. D., Xavier, J. B. & Foster, K. R. The sociobiology of biofilms. FEMS Microbiol. Rev. 33, 206–224 (2009).Article 
CAS 

Google Scholar 
Fritts, R. K., McCully, A. L. & McKinlay, J. B. Extracellular metabolism sets the table for microbial cross-feeding. Microbiol. Mol. Biol. Rev. 85, 135 (2021).Article 

Google Scholar 
D’Souza, G. et al. Ecology and evolution of metabolic cross-feeding interactions in bacteria. Nat. Prod. Rep. 35, 455–488 (2018).Article 

Google Scholar 
Libby, E., Hébert-Dufresne, L., Hosseini, S.-R. & Wagner, A. Syntrophy emerges spontaneously in complex metabolic systems. PLoS Comput. Biol. 15, e1007169 (2019).Article 

Google Scholar 
Staley, J. T. & Konopka, A. Measurement of in situ activities of nonphotosynthetic microorganisms in aquatic and terrestrial habitats. Annu. Rev. Microbiol. 39, 321–346 (1985).Article 
CAS 

Google Scholar 
Zachar, I. Closing the energetics gap. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01839-3 (2022).Article 

Google Scholar 
Zachar, I. & Boza, G. Endosymbiosis before eukaryotes: mitochondrial establishment in protoeukaryotes. Cell. Mol. Life Sci. 77, 3503–3523. https://doi.org/10.1007/s00018-020-03462-6 (2020).Article 
CAS 

Google Scholar 
Zachar, I. & Szathmáry, E. Breath-giving cooperation: critical review of origin of mitochondria hypotheses. Biol. Direct 12, 19. https://doi.org/10.1186/s13062-017-0190-5 (2017).Article 
CAS 

Google Scholar 
Booth, A. & Doolittle, W. F. Eukaryogenesis, how special really?. Proc. Natl. Acad. Sci. 112, 10278–10285 (2015).Article 
ADS 
CAS 

Google Scholar 
Morris, B. E. L., Henneberger, R., Huber, H. & Moissl-Eichinger, C. Microbial syntrophy: Interaction for the common good. FEMS Microbiol. Rev. 37, 384–406 (2013).Article 
CAS 

Google Scholar 
Szathmáry, E. On the propagation of a conceptual error concerning hypercycles and cooperation. J. Syst. Chem. 4, 2208 (2013).Article 

Google Scholar 
Seth, E. C. & Taga, M. E. Nutrient cross-feeding in the microbial world. Front. Microbiol. 5, 350 (2014).Article 

Google Scholar 
Piccardi, P., Vessman, B. & Mitri, S. Toxicity drives facilitation between 4 bacterial species. Proc. Natl. Acad. Sci. 116, 15979–15984 (2019).Article 
ADS 
CAS 

Google Scholar 
Yurtsev, E. A., Conwill, A. & Gore, J. Oscillatory dynamics in a bacterial cross-protection mutualism. Proc. Natl. Acad. Sci. 113, 6236–6241 (2016).Article 
ADS 
CAS 

Google Scholar 
Kehe, J. et al. Positive interactions are common among culturable bacteria. Sci. Adv. 7, 45 (2021).Article 

Google Scholar 
Momeni, B., Xie, L. & Shou, W. Lotka-Volterra pairwise modeling fails to capture diverse pairwise microbial interactions. Elife 6, 25051 (2017).Article 

Google Scholar 
Zengler, K. & Zaramela, L. S. The social network of microorganisms: How auxotrophies shape complex communities. Nat. Rev. Microbiol. 16, 383–390 (2018).Article 
CAS 

Google Scholar 
Koschwanez, J. H., Foster, K. R. & Murray, A. W. Sucrose utilization in budding yeast as a model for the origin of undifferentiated multicellularity. PLoS Biol. 9, e1001122 (2011).Article 
CAS 

Google Scholar 
Ciofu, O., Beveridge, T. J., Kadurugamuwa, J., Walther-Rasmussen, J. & Høiby, N. Chromosomal beta-lactamase is packaged into membrane vesicles and secreted from Pseudomonas aeruginosa. J. Antimicrob. Chemother. 45, 9–13 (2000).Article 
CAS 

Google Scholar 
Xenophontos, C., Harpole, W. S., Küsel, K. & Clark, A. T. Cheating promotes coexistence in a two-species one-substrate culture model. Front. Ecol. Evol. 9, 78006 (2022).Article 

Google Scholar 
West, S. A., Diggle, S. P., Buckling, A., Gardner, A. & Griffin, A. S. The social lives of microbes. Annu. Rev. Ecol. Evol. Syst. 38, 53–77 (2007).Article 

Google Scholar 
Flemming, H.-C. & Wingender, J. The biofilm matrix. Nat. Rev. Microbiol. 8, 623–633 (2010).Article 
CAS 

Google Scholar 
Kümmerli, R. & Brown, S. P. Molecular and regulatory properties of a public good shape the evolution of cooperation. Proc. Natl. Acad. Sci. 107, 18921–18926 (2010).Article 
ADS 

Google Scholar 
Griffin, A. S., West, S. A. & Buckling, A. Cooperation and competition in pathogenic bacteria. Nature 430, 1024–1027 (2004).Article 
ADS 
CAS 

Google Scholar 
Kramer, J., Özkaya, Ö. & Kümmerli, R. Bacterial siderophores in community and host interactions. Nat. Rev. Microbiol. 18, 152–163 (2019).Article 

Google Scholar 
van der Meij, A., Worsley, S. F., Hutchings, M. I. & van Wezel, G. P. Chemical ecology of antibiotic production by actinomycetes. FEMS Microbiol. Rev. 41, 392–416 (2017).Article 

Google Scholar 
Kümmerli, R., Schiessl, K. T., Waldvogel, T., McNeill, K. & Ackermann, M. Habitat structure and the evolution of diffusible siderophores in bacteria. Ecol. Lett. 17, 1536–1544 (2014).Article 

Google Scholar 
Jautzus, T., van Gestel, J. & Kovács, Á. T. Complex extracellular biology drives surface competition in lessigreaterBacillus subtilisless/igreater. Ecol. Lett. 16, 2320–2328. https://doi.org/10.1101/2022.02.28.482363 (2022).Article 
CAS 

Google Scholar 
Sachs, J. L., Mueller, U. G., Wilcox, T. P. & Bull, J. J. The evolution of cooperation. Q. Rev. Biol. 79, 135–160 (2004).Article 

Google Scholar 
Hillesland, K. L. & Stahl, D. A. Rapid evolution of stability and productivity at the origin of a microbial mutualism. Proc. Natl. Acad. Sci. 107, 2124–2129 (2010).Article 
ADS 
CAS 

Google Scholar 
Bruno, J. F., Stachowicz, J. J. & Bertness, M. D. Inclusion of facilitation into ecological theory. Trends Ecol. Evol. 18, 119–125 (2003).Article 

Google Scholar 
Gore, J., Youk, H. & van Oudenaarden, A. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256 (2009).Article 
ADS 
CAS 

Google Scholar 
Sorg, R. A. et al. Collective resistance in microbial communities by intracellular antibiotic deactivation. PLoS Biol. 14, e2000631 (2016).Article 

Google Scholar 
Karray, F. et al. Extracellular hydrolytic enzymes produced by halophilic bacteria and archaea isolated from hypersaline lake. Mol. Biol. Rep. 45, 1297–1309 (2018).Article 
CAS 

Google Scholar 
Datta, M. S., Sliwerska, E., Gore, J., Polz, M. F. & Cordero, O. X. Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat. Commun. 7, 11965 (2016).Article 
ADS 
CAS 

Google Scholar 
Tarnita, C. E. The ecology and evolution of social behavior in microbes. J. Exp. Biol. 220, 18–24 (2017).Article 

Google Scholar 
Özkaya, Ö., Xavier, K. B., Dionisio, F. & Balbontn, R. Maintenance of microbial cooperation mediated by public goods in single- and multiple-trait scenarios. J. Bacteriol. 199, 22 (2017).Article 

Google Scholar 
Yang, D.-D. et al. Fitness and productivity increase with ecotypic diversity among Escherichia coli strains that coevolved in a simple, constant environment. Appl. Environ. Microbiol. 86, 8 (2020).Article 

Google Scholar 
Pande, S. et al. Fitness and stability of obligate cross-feeding interactions that emerge upon gene loss in bacteria. ISME J. 8, 953–962 (2013).Article 

Google Scholar 
Zhou, K., Qiao, K., Edgar, S. & Stephanopoulos, G. Distributing a metabolic pathway among a microbial consortium enhances production of natural products. Nat. Biotechnol. 33, 377–383 (2015).Article 
CAS 

Google Scholar 
Harcombe, W. R., Chacón, J. M., Adamowicz, E. M., Chubiz, L. M. & Marx, C. J. Evolution of bidirectional costly mutualism from byproduct consumption. Proc. Natl. Acad. Sci. 115, 12000–12004 (2018).Article 
ADS 
CAS 

Google Scholar 
Summers, Z. M. et al. Direct exchange of electrons within aggregates of an evolved syntrophic coculture of anaerobic bacteria. Science 330, 1413–1415 (2010).Article 
ADS 
CAS 

Google Scholar 
Maddamsetti, R., Lenski, R. E. & Barrick, J. E. Adaptation, clonal interference, and frequency-dependent interactions in a long-term evolution experiment with Escherichia coli. Genetics 200, 619–631 (2015).Article 
CAS 

Google Scholar 
Gerrish, P. J. & Lenski, R. E. The fate of competing beneficial mutations in an asexual population. Genetica 102, 127–144 (1998).Article 

Google Scholar 
Popat, R. et al. Quorum-sensing and cheating in bacterial biofilms. Proc. R. Soc. B 279, 4765–4771 (2012).Article 
CAS 

Google Scholar 
Rainey, P. B. & Rainey, K. Evolution of cooperation and conflict in experimental bacterial populations. Nature 425, 72–74 (2003).Article 
ADS 
CAS 

Google Scholar 
Hardin, G. Tragedy of the commons. Science 162, 1243 (1968).Article 
ADS 
CAS 

Google Scholar 
West, S. A., Cooper, G. A., Ghoul, M. B. & Ten Griffin, A. S. recent insights for our understanding of cooperation. Nat. Ecol. Evol. 5, 419–430 (2021).Article 

Google Scholar 
MacArthur, R. Species packing and competitive equilibrium for many species. Theor. Popul. Biol. 1, 1–11 (1970).Article 
CAS 

Google Scholar 
Oliveira, N. M., Niehus, R. & Foster, K. R. Evolutionary limits to cooperation in microbial communities. Proc. Natl. Acad. Sci. 111, 17941–17946 (2014).Article 
ADS 
CAS 

Google Scholar 
Tilman, D. Resource Competition and Community Structure. Monographs in Population Biology, Vol. 17 (Princeton University Press, 1982).
Google Scholar 
Ferenci, T. Trade-off mechanisms shaping the diversity of bacteria. Trends Microbiol. 24, 209–223 (2016).Article 
CAS 

Google Scholar 
Rozen, D. E., Philippe, N., de Visser, J. A., Lenski, R. E. & Schneider, D. Death and cannibalism in a seasonal environment facilitate bacterial coexistence. Ecol. Lett. 12, 34–44 (2009).Article 

Google Scholar 
Brännström, Å., Johansson, J. & von Festenberg, N. The Hitchhiker’s Guide to Adaptive Dynamics. Games 4, 304–328 (2013).Article 
MATH 

Google Scholar 
Ramin, K. I. & Allison, S. D. Bacterial tradeoffs in growth rate and extracellular enzymes. Front. Microbiol. 10, 2956 (2019).Article 

Google Scholar 
Imachi, H. et al. Isolation of an archaeon at the prokaryote–eukaryote interface. Nature 577, 519–525 (2020).Article 
ADS 
CAS 

Google Scholar 
Wintermute, E. H. & Silver, P. A. Emergent cooperation in microbial metabolism. Mol. Syst. Biol. 6, 407 (2010).Article 

Google Scholar 
Libby, E., Kempes, C. & Okie, J. Metabolic compatibility and the rarity of prokaryote endosymbioses. BioRxiv https://doi.org/10.1101/2022.04.14.488272 (2022).Article 

Google Scholar 
Pauli, B., Oña, L., Hermann, M. & Kost, C. Obligate mutualistic cooperation limits evolvability. Nat. Commun. 13, 27630 (2022).Article 

Google Scholar 
Oña, L. & Kost, C. Cooperation increases robustness to ecological disturbance in microbial cross-feeding networks. Ecol. Lett. 25, 1410–1420 (2022).Article 

Google Scholar 
Machado, D. et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat. Ecol. Evol. 5, 195–203 (2021).Article 

Google Scholar 
Mee, M. T., Collins, J. J., Church, G. M. & Wang, H. H. Syntrophic exchange in synthetic microbial communities. Proc. Natl. Acad. Sci. 111, E2149–E2156 (2014).Article 
ADS 
CAS 

Google Scholar 
Goldford, J. E. et al. Emergent simplicity in microbial community assembly. Science 361, 469–474 (2018).Article 
ADS 
CAS 

Google Scholar 
McCutcheon, J. P. The genomics and cell biology of host-beneficial intracellular infections. Annu. Rev. Cell Dev. Biol. 37, 115–142 (2021).Article 
CAS 

Google Scholar 
Sousa, F. L., Neukirchen, S., Allen, J. F., Lane, N. & Martin, W. F. Lokiarchaeon is hydrogen dependent. Nat. Microbiol. 1, 5 (2016).Article 

Google Scholar 
Spang, A. et al. Proposal of the reverse flow model for the origin of the eukaryotic cell based on comparative analyses of Asgard archaeal metabolism. Nat. Microbiol. 4, 1138–1148 (2019).Article 
CAS 

Google Scholar 
Martin, W. & Müller, M. The hydrogen hypothesis for the first eukaryote. Nature 392, 37–41 (1998).Article 
ADS 
CAS 

Google Scholar 
López-García, P. & Moreira, D. The Syntrophy hypothesis for the origin of eukaryotes revisited. Nat. Microbiol. 5, 655–667 (2020).Article 

Google Scholar 
Mills, D. B. et al. Eukaryogenesis and oxygen in Earth history. Nat. Ecol. Evol. 6, 520–532 (2022).Article 

Google Scholar 
Liu, Y. et al. Expanded diversity of Asgard archaea and their relationships with eukaryotes. Nature 593, 553–557 (2021).Article 
ADS 
CAS 

Google Scholar 
Zachar, I., Szilágyi, A., Számadó, S. & Szathmáry, E. Farming the mitochondrial ancestor as a model of endosymbiotic establishment by natural selection. Proc. Natl. Acad. Sci. USA. 115, E1504–E1510. https://doi.org/10.1073/pnas.1718707115 (2018).Article 
ADS 
CAS 

Google Scholar 
Cavalier-Smith, T. & Chao, E.E.-Y. Multidomain ribosomal protein trees and the planctobacterial origin of neomura (eukaryotes, archaebacteria). Protoplasma https://doi.org/10.1007/s00709-019-01442-7 (2020).Article 

Google Scholar 
Searcy, D. G. Nutritional syntrophies and consortia as models for the origin of mitochondria. Symb. Mech. Model Syst. 1, 163–183. https://doi.org/10.1007/0-306-48173-1_10 (2002).Article 

Google Scholar 
Müller, N., Timmers, P., Plugge, C. M., Stams, A. J. M. & Schink, B. Syntrophy in methanogenic degradation. Endosymb. Methanog. Archaea 1, 153–192. https://doi.org/10.1007/978-3-319-98836-8_9 (2018).Article 

Google Scholar 
Searcy, D. G. Metabolic integration during the evolutionary origin of mitochondria. Cell Res. 13, 229–238 (2003).Article 
CAS 

Google Scholar 
Flemming, H.-C. & Wuertz, S. Bacteria and archaea on Earth and their abundance in biofilms. Nat. Rev. Microbiol. 17, 247–260 (2019).Article 
CAS 

Google Scholar 
Spang, A. et al. Asgard archaea are the closest prokaryotic relatives of eukaryotes. PLoS Genet. 14, e1007080 (2018).Article 

Google Scholar 
Burns, J. A., Pittis, A. A. & Kim, E. Gene-based predictive models of trophic modes suggest Asgard archaea are not phagocytotic. Nat. Ecol. Evol. 2, 697–704 (2018).Article 

Google Scholar 
Seitz, K. W. et al. Asgard archaea capable of anaerobic hydrocarbon cycling. Nat. Commun. 10, 1 (2019).Article 
ADS 
CAS 

Google Scholar 
Jimenez, P. & Scheuring, I. Density-dependent private benefit leads to bacterial mutualism. Evolution 75, 1619–1635. https://doi.org/10.1111/evo.14241 (2021).Article 

Google Scholar 
Preussger, D., Giri, S., Muhsal, L. K., Oña, L. & Kost, C. Reciprocal fitness feedbacks promote the evolution of mutualistic cooperation. Curr. Biol. 30, 3580-3590.e7 (2020).Article 
CAS 

Google Scholar 
Monaco, H. et al. Spatial-temporal dynamics of a microbial cooperative behavior resistant to cheating. Nat. Commun. 13, 3580 (2022).Article 
ADS 

Google Scholar 
Yanni, D., Márquez-Zacarias, P., Yunker, P. J. & Ratcliff, W. C. Drivers of spatial structure in social microbial communities. Curr. Biol. 29, 545–550 (2019).Article 

Google Scholar  More

D’Costa, V. M. et al. Antibiotic resistance is ancient. Nature 477, 457–461. https://doi.org/10.1038/nature10388 (2011).Article 
ADS 

Google Scholar 
Cytryn, E. The soil resistome: The anthropogenic, the native, and the unknown. Soil Biol. Biochem. 63, 18–23. https://doi.org/10.1016/j.soilbio.2013.03.017 (2013).Article 

Google Scholar 
Holmes, A. H. et al. Understanding the mechanisms and drivers of antimicrobial resistance. Lancet 387, 176–187. https://doi.org/10.1016/S0140-6736(15)00473-0 (2016).Article 

Google Scholar 
Regulation (EC) No 1831/2003 of the European parliament and of the council of 22 September 2003 on additives for use in animal nutrition.European Commission. Communication from the commission to the European parliament, the council, the European economic and social committee and the committee of the regions: A farm to fork strategy for a fair, healthy and environmentally-friendly food system COM/2020/381 Final, (2020).Kumar, K. C., Gupta, S. C., Chander, Y. & Singh, A. K. Antibiotic use in agriculture and its impact on the terrestrial environment. Adv. Agron. 87, 1–54. https://doi.org/10.1016/S0065-2113(05)87001-4 (2005).Article 

Google Scholar 
Chee-Sanford, J. C. et al. Fate and transport of antibiotic residues and antibiotic resistance genes following land application of manure waste. J. Environ. Qual. 38, 1086–1108. https://doi.org/10.2134/jeq2008.0128 (2009).Article 

Google Scholar 
Heuer, H., Schmitt, H. & Smalla, K. Antibiotic resistance gene spread due to manure application on agricultural fields. Curr. Opin. Microbiol. 14, 236–243. https://doi.org/10.1016/j.mib.2011.04.009 (2011).Article 

Google Scholar 
Epelde, L. et al. Characterization of composted organic amendments for agricultural use. Front. Sustain. Food Syst. 2, 44. https://doi.org/10.3389/fsufs.2018.00044 (2018).Article 

Google Scholar 
Youngquist, C. P., Mitchell, S. M. & Cogger, C. G. Fate of antibiotics and antibiotic resistance during digestion and composting: A review. J. Environ. Qual. 45, 537–545. https://doi.org/10.2134/jeq2015.05.0256 (2016).Article 

Google Scholar 
Ma, X., Xue, X., González-Mejía, A., Garland, J. & Cashdollar, J. Sustainable water systems for the city of tomorrow: A conceptual framework. Sustainability 7, 12071–12105. https://doi.org/10.3390/su70912071 (2015).Article 

Google Scholar 
Wang, Y. et al. Degradation of antibiotic resistance genes and mobile gene elements in dairy manure anerobic digestion. PLoS ONE 16, e0254836. https://doi.org/10.1371/journal.pone.0254836 (2021).Article 

Google Scholar 
Thanomsub, B. et al. Effects of ozone treatment on cell growth and ultrastructural changes in bacteria. J. Gen. Appl. Microbiol. 48, 193–199. https://doi.org/10.2323/jgam.48.193 (2002).Article 

Google Scholar 
Sousa, J. M. et al. Ozonation and UV254nm radiation for the removal of microorganisms and antibiotic resistance genes from urban wastewater. J. Hazard. Mater. 323, 434–441. https://doi.org/10.1016/j.jhazmat.2016.03.096 (2017).Article 

Google Scholar 
Park, J. H., Choppala, G. K., Bolan, N. S., Chung, J. W. & Chuasavathi, T. Biochar reduces the bioavailability and phytotoxicity of heavy metals. Plant Soil 348, 439–451. https://doi.org/10.1007/s11104-011-0948-y (2011).Article 

Google Scholar 
Jeffery, S. et al. The way forward in biochar research: targeting trade-offs between the potential wins. GCB Bioenergy 7, 1–13. https://doi.org/10.1111/gcbb.12132 (2015).Article 

Google Scholar 
Krasucka, P. et al. Engineered biochar: A sustainable solution for the removal of antibiotics from water. Chem. Eng. J. 405, 126926. https://doi.org/10.1016/j.cej.2020.126926 (2021).Article 

Google Scholar 
Ken, D. S. & Sinha, A. Recent developments in surface modification of Nano zero-valent iron (nZVI): remediation, toxicity and environmental impacts. Environ. Nanotechnol. Monit. Manag. 14, 100344. https://doi.org/10.1016/j.enmm.2020.100344 (2020).Article 

Google Scholar 
Zhao, X. et al. An overview of preparation and applications of stabilized zero-valent iron nanoparticles for soil and groundwater remediation. Water Res. 100, 245–266. https://doi.org/10.1016/j.watres.2016.05.019 (2016).Article 

Google Scholar 
Diao, M. & Yao, M. Use of zero-valent iron nanoparticles in inactivating microbes. Water Res. 43, 5243–5251. https://doi.org/10.1016/j.watres.2009.08.051 (2009).Article 

Google Scholar 
Shi, C. J., Wei, J., Jin, Y., Kniel, K. E. & Chiu, P. C. Removal of viruses and bacteriophages from drinking water using zero-valent iron. Sep. Purif. Technol. 84, 72–78. https://doi.org/10.1016/j.seppur.2011.06.036 (2012).Article 

Google Scholar 
Anza, M., Salazar, O., Epelde, L., Alkorta, I. & Garbisu, C. The application of nanoscale zero-valent iron promotes soil remediation while negatively affecting soil microbial biomass and activity. Front. Environ. Sci. 7, 19. https://doi.org/10.3389/fenvs.2019.00019 (2019).Article 

Google Scholar 
FAOSTAT. Mushrooms and truffles, production quantity (tons). https://www.tilasto.com/en/topic/geography-and-agriculture/crop/mushrooms-and-truffles/mushrooms-and-truffles-production-quantity/spain, (2020).Polat, E., Uzun, H., Topçuo, B., Önal, K. & Onus, A. N. Effects of spent mushroom compost on quality and productivity of cucumber (Cucumis sativus L.) grown in greenhouses. Afr. J. Biotechnol. 8, 176–180 (2009).
Google Scholar 
Fazaeli, H. & Masoodi, A. R. T. Spent wheat straw compost of Agaricus bisporus mushroom as ruminant feed. Asian-Australas. J. Anim. Sci. 19, 845–851. https://doi.org/10.5713/ajas.2006.845 (2006).Article 

Google Scholar 
Phan, C. W. & Sabaratnam, V. Potential uses of spent mushroom substrate and its associated lignocellulosic enzymes. Appl. Microbiol. Biotechnol. 96, 863–873. https://doi.org/10.1007/s00253-012-4446-9 (2012).Article 

Google Scholar 
Lau, K. L., Tsang, Y. Y. & Chiu, S. W. Use of spent mushroom compost to bioremediate PAH-contaminated samples. Chemosphere 52, 1539–1546. https://doi.org/10.1016/S0045-6535(03)00493-4 (2003).Article 
ADS 

Google Scholar 
Mayans, B. et al. An assessment of Pleurotus ostreatus to remove sulfonamides, and its role as a biofilter based on its own spent mushroom substrate. Environ. Sci. Pollut. Res. Int. 28, 7032–7042. https://doi.org/10.1007/s11356-020-11078-3 (2021).Article 

Google Scholar 
Congilosi, J. L. & Aga, D. S. Review on the fate of antimicrobials, antimicrobial resistance genes, and other micropollutants in manure during enhanced anaerobic digestion and composting. J. Hazard. Mater. 405, 123634. https://doi.org/10.1016/j.jhazmat.2020.123634 (2021).Article 

Google Scholar 
Oliver, J. P. et al. Invited review: fate of antibiotic residues, antibiotic-resistant bacteria, and antibiotic resistance genes in US dairy manure management systems. J. Dairy Sci. 103, 1051–1071. https://doi.org/10.3168/jds.2019-16778 (2020).Article 

Google Scholar 
Beneragama, N. et al. Survival of multidrug-resistant bacteria in thermophilic and mesophilic anaerobic co-digestion of dairy manure and waste milk. Anim. Sci. J. 84, 426–433. https://doi.org/10.1111/asj.12017 (2013).Article 

Google Scholar 
Sun, W., Qian, X., Gu, J., Wang, X. J. & Duan, M. L. Mechanism and effect of temperature on variations in antibiotic resistance genes during anaerobic digestion of dairy manure. Sci. Rep. 6, 30237. https://doi.org/10.1038/srep30237 (2016).Article 
ADS 

Google Scholar 
Sun, W., Gu, J., Wang, X., Qian, X. & Peng, H. Solid-state anaerobic digestion facilitates the removal of antibiotic resistance genes and mobile genetic elements from cattle manure. Bioresour. Technol. 274, 287–295. https://doi.org/10.1016/j.biortech.2018.09.013 (2019).Article 

Google Scholar 
Zou, Y., Xiao, Y., Wang, H., Fang, T. & Dong, P. New insight into fates of sulfonamide and tetracycline resistance genes and resistant bacteria during anaerobic digestion of manure at thermophilic and mesophilic temperatures. J. Hazard. Mater. 384, 121433. https://doi.org/10.1016/j.jhazmat.2019.121433 (2020).Article 

Google Scholar 
Agga, G. E., Kasumba, J., Loughrin, J. H. & Conte, E. D. Anaerobic digestion of tetracycline spiked livestock manure and poultry litter increased the abundances of antibiotic and heavy metal resistance genes. Front Microbiol. 11, 614424. https://doi.org/10.3389/fmicb.2020.614424 (2020).Article 

Google Scholar 
Jauregi, L., Epelde, L., González, A., Lavín, J. L. & Garbisu, C. Reduction of the resistome risk from cow slurry and manure microbiomes to soil and vegetable microbiomes. Environ. Microbiol. 23, 7643–7660. https://doi.org/10.1111/1462-2920.15842 (2021).Article 

Google Scholar 
Zhang, Z. et al. Assessment of global health risk of antibiotic resistance genes. Nat Commun 13, 1553. https://doi.org/10.1038/s41467-022-29283-8 (2022).Article 
ADS 

Google Scholar 
He, Y. et al. Antibiotic resistance genes from livestock waste: occurrence, dissemination, and treatment. npj Clean Water 3, 4. https://doi.org/10.1038/s41545-020-0051-0 (2020).Article 

Google Scholar 
Cui, E., Wu, Y., Zuo, Y. & Chen, H. Effect of different biochars on antibiotic resistance genes and bacterial community during chicken manure composting. Bioresour. Technol. 203, 11–17. https://doi.org/10.1016/j.biortech.2015.12.030 (2016).Article 

Google Scholar 
Fu, Y., Zhang, A., Guo, T., Zhu, Y. & Shao, Y. Biochar and hyperthermophiles as additives accelerate the removal of antibiotic resistance genes and mobile genetic elements during composting. Materials (Basel) 14, 5428. https://doi.org/10.3390/ma14185428 (2021).Article 
ADS 

Google Scholar 
Forsberg, K. J. et al. Bacterial phylogeny structures soil resistomes across habitats. Nature 509, 612–616. https://doi.org/10.1038/nature13377 (2014).Article 
ADS 

Google Scholar 
Li, H. et al. Effects of bamboo charcoal on antibiotic resistance genes during chicken manure composting. Ecotoxicol. Environ. Saf. 140, 1–6. https://doi.org/10.1016/j.ecoenv.2017.01.007 (2017).Article 
ADS 

Google Scholar 
Bondarenko, O., Ivask, A., Käkinen, A. & Kahru, A. Sub-toxic effects of CuO nanoparticles on bacteria: Kinetics, role of Cu ions and possible mechanisms of action. Environ. Pollut. 169, 81–89. https://doi.org/10.1016/j.envpol.2012.05.009 (2012).Article 

Google Scholar 
Wang, X. et al. Bacterial exposure to ZnO nanoparticles facilitates horizontal transfer of antibiotic resistance genes. NanoImpact 10, 61–67. https://doi.org/10.1016/j.impact.2017.11.006 (2018).Article 
ADS 

Google Scholar 
Qiu, X., Zhou, G. & Wang, H. Nanoscale zero-valent iron inhibits the horizontal gene transfer of antibiotic resistance genes in chicken manure compost. J. Hazard. Mater. 422, 126883. https://doi.org/10.1016/j.jhazmat.2021.126883 (2022).Article 

Google Scholar 
Zeng, T., Wilson, C. J. & Mitch, W. A. Effect of chemical oxidation on the sorption tendency of dissolved organic matter to a model hydrophobic surface. Environ. Sci. Technol. 48, 5118–5126. https://doi.org/10.1021/es405257b (2014).Article 
ADS 

Google Scholar 
Pak, G. et al. Comparison of antibiotic resistance removal efficiencies using ozone disinfection under different pH and suspended solids and humic substance concentrations. Environ. Sci. Technol. 50, 7590–7600. https://doi.org/10.1021/acs.est.6b01340 (2016).Article 
ADS 

Google Scholar 
Zhuang, Y. et al. Inactivation of antibiotic resistance genes in municipal wastewater by chlorination, ultraviolet, and ozonation disinfection. Environ. Sci. Pollut. Res. Int. 22, 7037–7044. https://doi.org/10.1007/s11356-014-3919-z (2015).Article 

Google Scholar 
Park, S., Rana, A., Sung, W. & Munir, M. Competitiveness of quantitative polymerase chain reaction (qPCR) and droplet digital polymerase chain reaction (ddPCR) technologies, with a particular focus on detection of antibiotic resistance genes (ARGs). Appl. Microbiol. 1, 426–444. https://doi.org/10.3390/applmicrobiol1030028 (2021).Article 

Google Scholar 
European Medicines Agency. European surveillance of veterinary antimicrobial consumption, (2020). Sales of Veterinary Antimicrobial Agents in 31 European Countries in 2018 (EMA/24309/2020).Heuer, H. et al. The complete sequences of plasmids pB2 and pB3 provide evidence for a recent ancestor of the IncP-1β group without any accessory genes. Microbiology (Reading) 150, 3591–3599. https://doi.org/10.1099/mic.0.27304-0 (2004).Article 

Google Scholar 
World Health Organization. Critically Important Antimicrobials for Human Medicine, 6th Revision (WHO, Geneva, Switzerland, 2019).Zhu, Y. G. et al. Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc. Natl Acad. Sci. U. S. A. 110, 3435–3440. https://doi.org/10.1073/pnas.1222743110 (2013).Article 
ADS 

Google Scholar 
Guo, T. et al. Increased occurrence of heavy metals, antibiotics and resistance genes in surface soil after long-term application of manure. Sci. Total Environ. 635, 995–1003. https://doi.org/10.1016/j.scitotenv.2018.04.194 (2018).Article 
ADS 

Google Scholar 
Nõlvak, H. et al. Inorganic and organic fertilizers impact the abundance and proportion of antibiotic resistance and integron-integrase genes in agricultural grassland soil. Sci. Total Environ. 562, 678–689. https://doi.org/10.1016/j.scitotenv.2016.04.035 (2016).Article 
ADS 

Google Scholar 
Chen, Q. L. et al. Effect of biochar amendment on the alleviation of antibiotic resistance in soil and phyllosphere of Brassica chinensis L.. Soil Biol. Biochem. 119, 74–82. https://doi.org/10.1016/j.soilbio.2018.01.015 (2018).Article 

Google Scholar 
Zhu, B., Chen, Q., Chen, S. & Zhu, Y. G. Does organically produced lettuce harbor higher abundance of antibiotic resistance genes than conventionally produced?. Environ. Int. 98, 152–159. https://doi.org/10.1016/j.envint.2016.11.001 (2017) .Article 

Google Scholar 
Métodos, M. A. P. A. Oficiales de análisis de suelos y Aguas Para riego. Minist. Agric. Pesca Aliment. Métodos Oficiales Anal. III (1994).Muziasari, W. I. et al. Aquaculture changes the profile of antibiotic resistance and mobile genetic element associated genes in Baltic Sea sediments. FEMS Microbiol. Ecol. 92, fiw052. https://doi.org/10.1093/femsec/fiw052 (2016).Article 

Google Scholar 
Muurinen, J. et al. Influence of manure application on the environmental resistome under Finnish agricultural practice with restricted antibiotic use. Environ. Sci. Technol. 51, 5989–5999. https://doi.org/10.1021/acs.est.7b00551 (2017).Article 
ADS 

Google Scholar 
Muziasari, W. I. et al. The resistome of farmed fish feces contributes to the enrichment of antibiotic resistance genes in sediments below Baltic Sea fish farms. Front. Microbiol. 7, 2137. https://doi.org/10.3389/fmicb.2016.02137 (2017).Article 

Google Scholar 
Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408. https://doi.org/10.1006/meth.2001.1262 (2001).Article 

Google Scholar 
Ovreås, L., Forney, L., Daae, F. L. & Torsvik, V. Distribution of bacterioplankton in meromictic Lake Saelenvannet, as determined by denaturing gradient gel electrophoresis of PCR-amplified gene fragments coding for 16S rRNA. Appl. Environ. Microbiol. 63, 3367–3373. https://doi.org/10.1128/aem.63.9.3367-3373.1997 (1997) .Article 
ADS 

Google Scholar 
Caporaso, J. G. et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624. https://doi.org/10.1038/ismej.2012.8 (2012).Article 

Google Scholar 
Lanzén, A. et al. Multi-targeted metagenetic analysis of the influence of climate and environmental parameters on soil microbial communities along an elevational gradient. Sci. Rep. 6, 28257. https://doi.org/10.1038/srep28257 (2016).Article 
ADS 

Google Scholar 
Pinna, N. K., Dutta, A., Monzoorul, H. M. & Mande, S. S. Can targeting non-contiguous V-regions with paired-end sequencing improve 16S rRNA-based taxonomic resolution of microbiomes?: An in silico evaluation. Front. Genet. 10, 653. https://doi.org/10.3389/fgene.2019.00653 (2019).Article 

Google Scholar 
Andrews, S. FastQC: A quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. https://doi.org/10.14806/ej.17.1.200 (2011).Article 

Google Scholar 
Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. https://doi.org/10.1038/s41587-019-0209-9 (2019).Article 

Google Scholar 
Amir, A. et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems 2, e00191-e216. https://doi.org/10.1128/mSystems.00191-16 (2017).Article 

Google Scholar 
Yang, Y., Li, B., Zou, S., Fang, H. H. P. & Zhang, T. Fate of antibiotic resistance genes in sewage treatment plant revealed by metagenomic approach. Water Res. 62, 97–106. https://doi.org/10.1016/j.watres.2014.05.019 (2014).Article 

Google Scholar 
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).Book 
MATH 

Google Scholar 
de Mendiburu, F. Agricolae: Statistical procedures for agricultural research. R package version 1.3-3. https://CRAN.R-project.org/package=agricolae (2020).Paradis, E. & Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528. https://doi.org/10.1093/bioinformatics/bty633 (2019).Article 

Google Scholar 
Oksanen, J. et al. Vegan: Community ecology package. R Package Version 2.3-1. (2015). More

Stefansson, S. O., Björnsson, B. T., Ebbesson, L. O. E. & McCormick, S. D. Smoltification. In Fish Larval Physiology (eds Finn, R. N. & Kapoor, B. G.) 639–681 (CRC Press, 2020).Chapter 

Google Scholar 
Kaleka, A. S., Kaur, N. & Bali, G. K. Larval development and molting. In Edible Insects (ed. Mikkola, H.) 17 (IntechOpen, 2019).
Google Scholar 
Butler, L. K. & Rohwer, V. G. Feathers and molt. in Ornithology: Foundation, Analysis, and Application (eds Morrison, M. L. et al.) 242–270 (JHU Press, 2018).
Google Scholar 
Swaddle, J. P., Witter, M. S., Cuthill, I. C., Budden, A. & McCowen, P. Plumage condition affects flight performance in common starlings: Implications for developmental homeostasis, abrasion and moult. J. Avian Biol. 27, 103–111 (1996).Article 

Google Scholar 
Norris, D. R., Marra, P. P., Montgomerie, R., Kyser, T. K. & Ratcliffe, L. M. Reproductive effort, molting latitude, and feather color in a migratory songbird. Science 306, 2249–2250 (2004).Article 
ADS 
CAS 

Google Scholar 
Delhey, K., Peters, A. & Kempenaers, B. Cosmetic coloration in birds: Occurrence, function, and evolution. Am. Nat. 169, S145–S158 (2007).Article 

Google Scholar 
Tomotani, B. M. & Muijres, F. T. A songbird compensates for wing molt during escape flights by reducing the molt gap and increasing angle of attack. J. Exp. Biol. 222, 195396 (2019).Article 

Google Scholar 
Galván, I., Negro, J. J., Rodriguez, A. & Carrascal, L. M. On showy dwarfs and sober giants: Body size as a constraint for the evolution of bird plumage colouration. Acta Ornithol. 48, 65–80 (2013).Article 

Google Scholar 
Speakman, J. R. & Król, E. Maximal heat dissipation capacity and hyperthermia risk: Neglected key factors in the ecology of endotherms. J. Anim. Ecol. 79, 726–746 (2010).
Google Scholar 
Wolf, B. O. & Walsberg, G. E. The role of the plumage in heat transfer processes of birds. Am. Zool. 40, 575–584 (2000).
Google Scholar 
Berthold, P. & Querner, U. Genetic basis of moult, wing length, and body weight in a migratory bird species, Sylvia atricapilla. Experientia 38, 801–802 (1982).Article 

Google Scholar 
Gwinner, E., Neusser, V., Engl, D., Schmidl, D. & Bals, L. Haltung, Zucht und Eiaufzucht afrikanischer und europäischer Schwarzkehlchen Saxicola torquata. Gefied. Welt 111, 118–120 (1987).
Google Scholar 
Berthold, P. & Querner, U. Microevolutionary aspects of bird migration based on experimental results. Isr. J. Ecol. Evol. 41, 377–385 (1995).
Google Scholar 
Helm, B. & Gwinner, E. Timing of postjuvenal molt in African (Saxicola torquata axillaris) and European (Saxicola torquata rubicola) stonechats: Effects of genetic and environmental factors. Auk 116, 589–603 (1999).Article 

Google Scholar 
Helm, B. & Gwinner, E. Timing of molt as a buffer in the avian annual cycle. Acta Zool. Sin. 52, 703–706 (2006).
Google Scholar 
Rohwer, S., Ricklefs, R. E., Rohwer, V. G. & Copple, M. M. Allometry of the duration of flight feather molt in birds. PLoS Biol. 7, e1000132 (2009).Article 

Google Scholar 
Jenni, L. & Winkler, R. The Biology of Moult in Birds (Bloomsbury Publishing, 2020).
Google Scholar 
Tonra, C. M. & Reudink, M. W. Expanding the traditional definition of molt-migration. Auk Ornithol. Adv. 135, 1123–1132 (2018).
Google Scholar 
Rohwer, S., Butler, L. K., Froehlich, D. R., Greenberg, R. & Marra, P. P. Ecology and demography of east–west differences in molt scheduling of Neotropical migrant passerines. Birds Two Worlds Ecol. Evol. Migr. (R. Greenb. PP Marra, Eds.). Johns Hopkins Univ. Press. Balt. Maryl., 87–105 (2005).Bensch, S., Åkesson, S. & Irwin, D. E. The use of AFLP to find an informative SNP: Genetic differences across a migratory divide in willow warblers. Mol. Ecol. 11, 2359–2366 (2002).Article 
CAS 

Google Scholar 
Ruegg, K. Genetic, morphological, and ecological characterization of a hybrid zone that spans a migratory divide. Evol. Int. J. Org. Evol. 62, 452–466 (2008).Article 

Google Scholar 
Delmore, K. E., Fox, J. W. & Irwin, D. E. Dramatic intraspecific differences in migratory routes, stopover sites and wintering areas, revealed using light-level geolocators. Proc. R. Soc. B Biol. Sci. 279, 4582–4589 (2012).Article 

Google Scholar 
Delmore, K. E. et al. Individual variability and versatility in an eco-evolutionary model of avian migration. Proc. R. Soc. B 287, 20201339 (2020).Article 

Google Scholar 
Procházka, P. et al. Across a migratory divide: divergent migration directions and non-breeding grounds of Eurasian reed warblers revealed by geolocators and stable isotopes. J. Avian Biol. 49, 012516 (2018).Article 

Google Scholar 
Bensch, S., Grahn, M., Müller, N., Gay, L. & Åkesson, S. Genetic, morphological, and feather isotope variation of migratory willow warblers show gradual divergence in a ring. Mol. Ecol. 18, 3087–3096 (2009).Article 

Google Scholar 
Rohwer, S. & Irwin, D. E. Molt, orientation, and avian speciation. Auk 128, 419–425 (2011).Article 

Google Scholar 
Pageau, C., Sonnleitner, J., Tonra, C. M., Shaikh, M. & Reudink, M. W. Evolution of winter molting strategies in European and North American migratory passerines. Ecol. Evol. 11, 13247–13258 (2021).Article 

Google Scholar 
Butler, L. K., Rohwer, S. & Rogers, M. Prebasic molt and molt-related movements in Ash-throated Flycatchers. Condor 108, 647–660 (2006).Article 

Google Scholar 
Barry, J. H., Butler, L. K., Rohwer, S. & Rohwer, V. G. Documenting molt-migration in Western Kingbird (Tyrannus verticalis) using two measures of collecting effort. Auk 126, 260–267 (2009).Article 

Google Scholar 
Hobson, K. A. & Wassenaar, L. I. Linking breeding and wintering grounds of neotropical migrant songbirds using stable hydrogen isotopic analysis of feathers. Oecologia 109, 142–148 (1996).Article 
ADS 
CAS 

Google Scholar 
Hobson, K. A. & Wassenaar, L. I. Tracking Animal Migration with Stable Isotopes (Academic Press, 2018).
Google Scholar 
Rubenstein, D. R. & Hobson, K. A. From birds to butterflies: Animal movement patterns and stable isotopes. Trends Ecol. Evol. 19, 256–263 (2004).Article 

Google Scholar 
Bearhop, S. et al. Assortative mating as a mechanism for rapid evolution of a migratory divide. Science 310, 502–504 (2005).Article 
ADS 
CAS 

Google Scholar 
Eppig, J. T. et al. The mouse genome database (MGD): Comprehensive resource for genetics and genomics of the laboratory mouse. Nucleic Acids Res. 40, D881–D886 (2012).Article 
CAS 

Google Scholar 
Contina, A., Bridge, E. S. & Kelly, J. F. Exploring novel candidate genes from the mouse genome informatics database: Potential implications for avian migration research. Integr. Zool. 11, 240 (2016).Article 

Google Scholar 
Yang, J. et al. Common SNPs explain a large proportion of the heritability for human height. Nat. Genet. 42, 565–569 (2010).Article 
CAS 

Google Scholar 
Thompson, C. W. Is the Painted Bunting actually two species? Problems determining species limits between allopatric populations. Condor 93, 987–1000 (1991).Article 

Google Scholar 
Contina, A., Bridge, E. S., Seavy, N. E., Duckles, J. M. & Kelly, J. F. Using geologgers to investigate bimodal isotope patterns in Painted Buntings (Passerina ciris). Auk 130, 265 (2013).Article 

Google Scholar 
Besozzi, E., Chew, B., Allen, D. C. & Contina, A. Stable isotope analysis of an aberrant Painted Bunting (Passerina ciris) feather suggests post-molt movements. Wilson J. Ornithol. 133, 151 (2021).Article 

Google Scholar 
Sharp, A. et al. Spatial and Temporal Scale-Dependence of the Strength of Migratory Connectivity in a North American Passerine. https://assets.researchsquare.com/files/rs-1483049/v1/72236b63-952d-4870-89e7-461056b8625b.pdf?c=1648893558 (2022).Pyle, P. et al. Temporal, spatial, and annual variation in the occurrence of molt-migrant passerines in the Mexican monsoon region. Condor 111, 583–590 (2009).Article 

Google Scholar 
Bridge, E. S., Fudickar, A. M., Kelly, J. F., Contina, A. & Rohwer, S. Causes of bimodal stable isotope signatures in the feathers of a molt-migrant songbird. Can. J. Zool. 89, 951 (2011).Article 
CAS 

Google Scholar 
Seutin, G., White, B. N. & Boag, P. T. Preservation of avian blood and tissue samples for DNA analyses. Can. J. Zool. 69, 82–90 (1991).Article 
CAS 

Google Scholar 
Ali, O. A. et al. RAD capture rapture: Flexible and efficient sequence-based genotyping. Genetics 202, 389–400 (2016).Article 
CAS 

Google Scholar 
Contina, A. et al. Characterization of SNP markers for the Painted Bunting (Passerina ciris) and their relevance in population differentiation and genome evolution studies. Conserv. Genet. Resour. 11, 5–10 (2019).Article 
ADS 

Google Scholar 
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: An analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140 (2013).Article 

Google Scholar 
Parker, P., Li, B., Li, H. & Wang, J. The genome of Darwin’s Finch (Geospiza fortis). Gigascience 10, 100040 (2012).
Google Scholar 
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).Article 
CAS 

Google Scholar 
Van der Auwera, G. A. et al. From FastQ data to high-confidence variant calls: The genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinform. 43, 1–33 (2013).
Google Scholar 
McKenna, A. et al. The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).Article 
CAS 

Google Scholar 
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).Article 
CAS 

Google Scholar 
Anderson, E. genoscapeRtools: Tools for Building Migratory Bird Genoscapes (2019).Purcell, S. et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).Article 
CAS 

Google Scholar 
Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).Article 
CAS 

Google Scholar 
Alexander, D. H. & Lange, K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinform. 12, 246 (2011).Article 

Google Scholar 
Francis, R. M. pophelper: An R package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 17, 27–32 (2017).Article 
CAS 

Google Scholar 
Chew, B., Kelly, J. & Contina, A. Stable isotopes in avian research: a step by step protocol to feather sample preparation for stable isotope analysis of carbon (δ13C), nitrogen (δ15N), and hydrogen (δ2H). Version 1.1. https://doi.org/10.17504/protocols.io.z2uf8ew (2019).Wassenaar, L. I. & Hobson, K. A. Comparative equilibration and online technique for determination of non-exchangeable hydrogen of keratins for use in animal migration studies. Isotopes Environ. Health Stud. 39(3), 211–217 (2003).Article 
CAS 

Google Scholar 
Bowen, G. J., Wassenaar, L. I. & Hobson, K. A. Global application of stable hydrogen and oxygen isotopes to wildlife forensics. Oecologia 143, 337–348 (2005).Article 
ADS 

Google Scholar 
R Core Team: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021).Wassenaar, L. I. & Hobson, K. A. Stable-hydrogen isotope heterogeneity in keratinous materials: Mass spectrometry and migratory wildlife tissue subsampling strategies. Rapid Commun. Mass Spectrom. 20, 2505–2510 (2006).Article 
ADS 
CAS 

Google Scholar 
Zhou, X., Carbonetto, P. & Stephens, M. Polygenic modeling with Bayesian sparse linear mixed models. PLoS Genet. 9, e1003264 (2013).Article 
CAS 

Google Scholar 
Guan, Y. & Stephens, M. Bayesian variable selection regression for genome-wide association studies and other large-scale problems. Ann. Appl. Stat. 5, 455 (2011).Article 
MathSciNet 
MATH 

Google Scholar 
Marchini, J., Cardon, L. R., Phillips, M. S. & Donnelly, P. The effects of human population structure on large genetic association studies. Nat. Genet. 36, 512–517 (2004).Article 
CAS 

Google Scholar 
Browning, S. R. & Browning, B. L. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am. J. Hum. Genet. 81, 1084–1097 (2007).Article 
CAS 

Google Scholar 
Chaves, J. A. et al. Genomic variation at the tips of the adaptive radiation of Darwin’s finches. Mol. Ecol. 25, 5282–5295 (2016).Article 
CAS 

Google Scholar 
Zhang, Y.-W. et al. mrMLM v4.0.2: An R platform for multi-locus genome-wide association studies. Genom. Proteom. Bioinform. 18, 481–487 (2020).Article 

Google Scholar 
Grabherr, M. G. et al. Genome-wide synteny through highly sensitive sequence alignment: Satsuma. Bioinformatics 26, 1145–1151 (2010).Article 
CAS 

Google Scholar 
Quinlan, A. R. & Hall, I. M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).Article 
CAS 

Google Scholar 
Ellis, N., Smith, S. J. & Pitcher, C. R. Gradient forests: Calculating importance gradients on physical predictors. Ecology 93, 156–168 (2012).Article 

Google Scholar 
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

Google Scholar 
Anderson, E. C. snps2assays: Prepare SNP Assay Orders from ddRAD or RAD Loci (2015).Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).Article 
CAS 

Google Scholar 
Patterson, N., Price, A. L. & Reich, D. Population structure and eigenanalysis. PLoS Genet. 2, e190 (2006).Article 

Google Scholar 
Price, A. L. et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909 (2006).Article 
CAS 

Google Scholar 
Ruegg, K. et al. Ecological genomics predicts climate vulnerability in an endangered southwestern songbird. Ecol. Lett. 21, 1085–1096 (2018).Article 

Google Scholar 
Bay, R. A. et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science 359, 83–86 (2018).Article 
ADS 
CAS 

Google Scholar 
Hedenström, A. Adaptations to migration in birds: Behavioural strategies, morphology and scaling effects. Philos. Trans. R. Soc. B Biol. Sci. 363, 287–299 (2008).Article 

Google Scholar 
Buehler, D. M. & Piersma, T. Travelling on a budget: Predictions and ecological evidence for bottlenecks in the annual cycle of long-distance migrants. Philos. Trans. R. Soc. B Biol. Sci. 363, 247–266 (2008).Article 

Google Scholar 
Schieltz, P. C. & Murphy, M. E. The contribution of insulation changes to the energy cost of avian molt. Can. J. Zool. 75, 396–400 (1997).Article 

Google Scholar 
Carling, M. D. & Thomassen, H. A. The role of environmental heterogeneity in maintaining reproductive isolation between hybridizing Passerina (Aves: Cardinalidae) buntings. Int. J. Ecol. 2012, 1–11 (2012).Article 

Google Scholar 
Irwin, D. E. Incipient ring speciation revealed by a migratory divide. Mol. Ecol. 18, 2923–2925 (2009).Article 

Google Scholar 
Thomas, D. W., Blondel, J., Perret, P., Lambrechts, M. M. & Speakman, J. R. Energetic and fitness costs of mismatching resource supply and demand in seasonally breeding birds. Science 291, 2598–2600 (2001).Article 
ADS 
CAS 

Google Scholar 
Rohwer, V. G., Rohwer, S. & Ortiz-Ramirez, M. F. Molt biology of resident and migrant birds of the monsoon region of west Mexico. Ornitol. Neotrop. 20, 565–584 (2009).
Google Scholar 
Bensch, S., Andersson, T. & Åkesson, S. Morphological and molecular variation across a migratory divide in willow warblers, Phylloscopus trochilus. Evolution 53, 1925–1935 (1999).Article 

Google Scholar 
Turbek, S. P., Scordato, E. S. C. & Safran, R. J. The role of seasonal migration in population divergence and reproductive isolation. Trends Ecol. Evol. 33, 164–175 (2018).Article 

Google Scholar 
Scordato, E. S. C. et al. Migratory divides coincide with reproductive barriers across replicated avian hybrid zones above the Tibetan Plateau. Ecol. Lett. 23, 231–241 (2020).Article 

Google Scholar 
Battey, C. J. et al. A migratory divide in the Painted Bunting (Passerina ciris). Am. Nat. 191, 259–268 (2018).Article 
CAS 

Google Scholar 
Contina, A. et al. Genetic structure of the Painted Bunting and its implications for conservation of migratory populations. Ibis 161, 372 (2019).Article 

Google Scholar 
Butler, L. K. The grass is always greener: Do monsoon rains matter for molt of the Vermilion Flycatcher (Pyrocephalus rubinus)? Auk 130, 297–307 (2013).Article 

Google Scholar 
Turbek, S. P. et al. A migratory divide spanning two continents is associated with genomic and ecological divergence. Evolution 76, 722 (2022).Article 

Google Scholar 
Dietz, M. W., Daan, S. & Masman, D. Energy requirements for molt in the kestrel Falco tinnunculus. Physiol. Zool. 65, 1217–1235 (1992).Article 

Google Scholar 
Vézina, F., Gustowska, A., Jalvingh, K. M., Chastel, O. & Piersma, T. Hormonal correlates and thermoregulatory consequences of molting on metabolic rate in a northerly wintering shorebird. Physiol. Biochem. Zool. 82, 129–142 (2009).Article 

Google Scholar 
Bazzi, G. et al. Candidate genes have sex-specific effects on timing of spring migration and moult speed in a long-distance migratory bird. Curr. Zool. 63, 479–486 (2017).CAS 

Google Scholar 
Busby, L. et al. Sonic hedgehog specifies flight feather positional information in avian wings. Development 147, 188821 (2020).Article 

Google Scholar 
Eichberger, T. et al. GLI2-specific transcriptional activation of the bone morphogenetic protein/Activin antagonist Follistatin in human epidermal cells. J. Biol. Chem. 283, 12426–12437 (2008).Article 
CAS 

Google Scholar 
Matzuk, M. M. et al. Multiple defects and perinatal death in mice deficient in follistatin. Nature 374, 360–363 (1995).Article 
ADS 
CAS 

Google Scholar 
Patel, K., Makarenkova, H. & Jung, H.-S. The role of long range, local and direct signalling molecules during chick feather bud development involving the BMPs, follistatin and the Eph receptor tyrosine kinase Eph-A4. Mech. Dev. 86, 51–62 (1999).Article 
CAS 

Google Scholar 
Nakamura, M. et al. Control of pelage hair follicle development and cycling by complex interactions between follistatin and activin. FASEB J. 17, 1–22 (2003).Article 
MathSciNet 

Google Scholar 
Pays, L., Charvet, I., Hemming, F. J. & Saxod, R. Close link between cutaneous nerve pattern development and feather morphogenesis demonstrated by experimental production of neo-apteria and ectopic feathers: Implication of chondroitin sulphate proteoglycans and other matrix molecules. Anat. Embryol. 195, 457–466 (1997).Article 
CAS 

Google Scholar 
Pyle, P., Saracco, J. F. & DeSante, D. F. Evidence of widespread movements from breeding to molting grounds by North American landbirds. Auk Ornithol. Adv. 135, 506–520 (2018).
Google Scholar 
De Mita, S. et al. Detecting selection along environmental gradients: Analysis of eight methods and their effectiveness for outbreeding and selfing populations. Mol. Ecol. 22, 1383–1399 (2013).Article 

Google Scholar 
Lotterhos, K. E. & Whitlock, M. C. Evaluation of demographic history and neutral parameterization on the performance of FST outlier tests. Mol. Ecol. 23, 2178–2192 (2014).Article 

Google Scholar 
Frichot, E., Schoville, S. D., de Villemereuil, P., Gaggiotti, O. E. & François, O. Detecting adaptive evolution based on association with ecological gradients: Orientation matters!. Heredity (Edinb.) 115, 22–28 (2015).Article 
CAS 

Google Scholar 
Trivedi, A. K. et al. Temperature alters the hypothalamic transcription of photoperiod responsive genes in induction of seasonal response in migratory redheaded buntings. Mol. Cell. Endocrinol. 493, 110454 (2019).Article 
CAS 

Google Scholar  More

Morphometric measurements are subjective and used to assess the carcass development and quantitatively measure the muscular distribution in the carcass with estimates of its conformation. In the present study there were not significative differences observed for these parameters or for carcass compactness index (CCI), inferring that the use of cactus pear silage as well as intermittent water supply combined or alone did not alter animal growth and/or carcass conformation, maintaining the muscle pattern achieved by the control diet (usual) and demonstrating body and carcass uniformity. Since animals used in this study were homogeneous and had similar age and body performance, as indicated by the carcass morphometric measurements and by the difference between the empty carcass and hot carcass weights, which resulted in the sum of head + limb with an average of 8.2 ± 0.13 kg between treatments, giving an idea that the animals were similar in chronological age, since the allometric growth of the body occurs from the extremities to the interior of the body.The significant difference between treatments with inclusion of cactus pear silage for hot carcass yield (HCY) and cold carcass yield (CCY) may be related to the weight of the full gastrointestinal tract, which showed higher values for animals fed with a higher proportion of Tifton 85 grass hay in the diet (0% CPS). Increasing the NDF content of the diet reduces the passage rate of digesta, and the emptying of the gastrointestinal tract (GT) that cause a distension of the rumen-reticulum and increase the weight of the gastrointestinal tract, resulting in lower HCY and consequently lower CCY. While the diets with inclusion of CPS increase NFC content, such as pectin, which have higher rates of rumen degradability and, higher rates of passage7,8,9.Measurements and evaluations carried out on the carcass, such as the carcass compactness index and loin eye area (LEA), are parameters that quantitatively measure the muscle distribution in the carcass, an edible part of greater financial return, which indicates the conformation of these animals3, while the body condition score (BCS) and the measure C, which are highly correlated, measure the distribution of fat on the carcass, giving an idea of the carcass finish, in which the higher these variables, the greater the proportion of fat that allows for less water loss due to carcass cooling10. These variables in the present study were also not influenced by the levels of cactus pear silage and water restrictions, presenting an overall mean of 0.17 kg/cm, 7.68 cm, 2.42 points and 0.7 mm respectively, and consequently did not influence the losses due to cooling, which presented an average loss of 1.48%.The main cuts of the goat carcass are the neck, leg, shoulder, loin, and rib. Their economic values differ, and their proportions become an important index to evaluate the carcass quality9. The cuts of greatest importance and commercial values are the leg and the loin, called noble cuts because they present greater yield and muscle tenderness, being interesting that they present a good proportion in the carcass, for providing greater edible tissue content, mainly muscle.Carcasses with similar weight tend to have equivalent proportions of cuts, as they exhibit isogonic growth. As the cold carcass weight (CCW) and the conformation of the animals were similar, with similar morphometric measurements, they had a direct relationship in the absence of an effect on commercial cuts.The commercial value of the carcass, whether through carcass yield and/or the proportions of the cuts, is also linked to tissue composition, thus the dissection of the leg represents an estimate of measuring the tissue composition of the carcass, in which is sought a greater proportion of muscle, intermediate proportion of fat and less bone in carcasses11. In this way, diets with cactus pear silage and the different levels of intermittent water supply resulted in the constancy in the amount of muscle, fat, and bone in legs of goats. The similarity in muscle proportion is related to the lack of effects on slaughter weight and CCW, as the weight of muscles is highly correlated to carcass weight. The average muscle yield was above 60% in all treatments, confirming that the animals showed good efficiency to the diets and adapted well to the water supply levels. Although the diets with cactus silage had high amounts of metabolizable energy (ME) and no difference in DM intake, the energy input was similar that not influencing carcass weights and carcass compactness index. That is, it did not influence muscle deposition in the carcass, probably due to synchronicity of energy and protein.As for the weight and proportion of bone tissue, it is believed that because this is a tissue with early development in relation to muscle and fat2, diets in the final stages of growth (average of 8 months) would hardly change their participation in the tissue composition, where the relationship of this tissue with the others is usually only increased when there are changes in the proportion of muscle and/or fat.Water restriction, as long as it is moderate and acute, mainly affects the loss of body water and not tissues, which does not cause deleterious effects on animal productivity and growth.The muscle:fat ratio indicates the state of leg fattening, while the muscle:bone ratio estimates the carcass muscularity, both being attributes of quality3. The similarity previously reported in the weight of fat, bone and muscle corroborates that these relationships also do not have differences. The same occurs for the leg muscularity index (LMI), due to the weight of the five muscles used to determine the index and the length of the femur which had been similar between the animals.Nevertheless, when considering fat as a percentage of participation in leg weight, it is possible to observe that the intermittency in water supply in both intervals (24 and 48 h) reduced the proportion of fat in the leg. Although in this research, the water supply levels did not affect the daily intake of dry matter from animals, with average intake of 650.67 g/kg DM, ranging from 599 to 682 g/kg DM between treatments7, during days of water deprivation, fat mobilization for energy availability may occur, possibly offsetting water stress and influencing not only feed intake, on these days of deprivation but also affecting energy metabolism, which results in the mobilization of energy reserves2.When the physicochemical composition of the meat was evaluated, it was observed that the diets and water supply levels probably did not affect the reserves of muscle glycogen during the pre-slaughter management as can be seen through pHinitial and pHfinal. The pHinitial right after slaughter should be close to neutrality, as well as in the live animal, indicating that the animal did not suffer from stress during the pre-slaughter period. The pHfinal, on the other hand, is expected to show a considerable variation, between 5.55 and 6.2 for goat meat; and due be inversely proportional to the concentration of muscle glycogen at the time of slaughter, that is, a more intense expenditure of glycogen stores results in less lactic acid production and higher pHfinal10,12,13. In this research, the pHfinal had an average of 5.74, a pH higher than the isoelectric point of muscle proteins (5.2–5.3). This result is favorable, since it is above the neutral charge and presenting an excessive negative charge that provides the repulsion of filaments, which allows water molecules to bind and improve the organoleptic characteristics of the meat, through succulence and texture of meat13 evaluated by cooking loss, moisture, and shear force, principally. The cooking loss (CL), moisture and shear force (SF) were within the values recommended (20–35% CL, moisture above 70% and SF up to 44.13 Newton (N) for goat meat) to classify the meat as soft and tender14. Statistically, interactions were found between the supply of silage and intermittent water supply, in which goats on a diet without cactus pear silage and without intermittent water supply showed higher values of cooking losses and shear force.Higher concentrations of collagen content and/or greater activities of calpastatin (which inhibit the action of calpains), as well as larger fascicles and greater number of fibers present in each muscle fascicle, as was visually observed in the meat of the animals in this research, can lead to reductions in meat tenderness15. Because goat carcasses are generally small, with low marbling degree and a thin layer of subcutaneous fat, there is rapid heat dissipation at the beginning of the post-mortem period, which can lead to cold shortening, muscle hardening, and less tender meats16.pHfinal of the meat has a high correlation with color parameters (L*—lightness, a*—redness, b*—yellowness and Chroma), as the pHfinal can affect the reaction of myoglobin to oxymyoglobin. The b* index in meat, on the other hand, may be related to the concentration of fat and/or the presence of carotenoids in the diet which can be affected by forage preservation processes, such as silage and hay, which significantly reduces by up to 80% carotenoids levels13. It is believed that the carotenoid concentrations in the diet of this study were similar between treatments and consequently in values of b* of meat. Values of a* and Chroma directly depend on the content and state of the heme pigments in the muscle, due to the chemical state of iron (Fe), playing an important role in meat color10. These parameters showed no significant difference between treatments, however, higher values of a* and Chroma in meat are desired, as a result of the increase in oxymyoglobin and decrease in metmyoglobin that provides the meat’s “bloom”. According to Dawson et al.17, the minimum critical value for meat luminosity (L*) is 34. Lower values of L are related to elevating pHfinal, which results in the high concentration of metmyoglobin, making the meat darker, which causes rejection by consumers for associating dark meat to as old meat.The meat’s presentation and more precisely its color is an important factor that can influence a consumer’s purchase decision, as it gives us the idea of freshness and meat’ quality. The L* and a* color parameters are the most representative for these characteristics18. Although in our research it did not have a significant effect on the color parameters, we can indicate that the meat obtained in this research would be well accepted by consumers, because Hopkins19 suggests that consumers will consider meat color acceptable when the L* value is equal to or exceeds 34, and a* value below 19 or equal to or exceeds 9.5 according to Khliji et al.18. In the present study, all values for L* remained above this aforementioned threshold and the values of a* remained within these values which suggests that meats from all diets and water supply levels had an acceptable color for consumers.When evaluating the chemical composition of meat, no significant differences were observed between treatments, except for the ash content, that remained above the average values found in the literature, which is 0.99–1.10%16. It is believed that because cactus pear is a rich source of Ca, Mg, K and with increasing level of cactus pear silage in the diet31, these minerals were consumed in larger amounts, which could have resulted in a higher proportion of minerals in the meat of animals that received 42% cactus pear silage.The lipid fatty acid profile in meat has a major impact on sensory properties and nutritional quality, influencing acceptance and health for consumers20,21. Intermittent water supply, cactus pear silage, and interaction between water supply and cactus pear silage did not influence most fatty acids present in the Longissimus lumborum muscle of the animals under study, except only a few saturated fatty acids e.g. docosanoic acid (C22:0), tricosanoic acid (C23:0), BCFA, anteiso-tridecanoic acid (C13:0 anteiso) and anteiso-pentadecanoic acid (C15:0 anteiso).Biohydrogenation of ruminal bacteria results in a circumstantial variety of fatty acids (FA), which will be absorbed in the intestine and later incorporated into the meat of goats. In addition to the diet and the biohydrogenation, the meat lipid profile can vary due to de novo synthesis, desaturation, duration of the feeding period and differences in pathways of various FA by the animal organism22.A high concentration of saturated fatty acids present in meat is not desirable, as there is evidence that saturated fatty acids, mainly C16:0, as well as myristic (C14:0) and lauric (C12:0) increase the blood cholesterol and low-density lipoproteins (LDL) concentration, due to interferences with hepatic LDL receptors23, however, in the studied treatments, there were no significant differences for these fatty acids. On the other hand, C18:0 has no impact on cholesterol levels, due to being poorly digested and easily desaturated to C18:1 by Δ9-desaturase24, present in the cell endoplasmic reticulum. This fatty acid is not harmful to health and is considered the only desirable SFA. As the levels of C18:0 in diets tend to be minimal, their main origin is the biohydrogenation of PUFA and de novo syntheses in diets with a high energy pattern25.In addition to carrying out the biohydrogenation process, ruminal bacteria synthesize a series of FA, mainly those of odd and branched chain, that comprise mainly the lipids of the bacterial membrane26,27, to maintain membrane fluidity. Linear odd-chains fatty acids are formed when propionyl-CoA, instead of acetyl-CoA, is used as a de novo synthesis initiator25. On the other hand, iso and anteiso FA are synthesized by the precursors branched-chain amino acids (valine, leucine, and isoleucine) and their corresponding branched- short-chain carboxylic acids (isobutyric, isovaleric and 2-methyl butyric acids)28.There is an increasing interest to study odd-and branched-chain fatty acids (OBCFAs) from animal products, mainly in milk due to its higher concentration compared to meat. Researchers reported that several OBCFAs have potential health benefits in humans29 as improved gut health30 and presenting anti-cancer activity31, as well as improve the sensory characteristics of the meat, providing a greater sensation of tenderness and juiciness, because BCFA content are associated with a less consistent fat in meat from lambs due to its lower melting point and its chain structure32.The FAs profile in the ruminal bacteria is largely composed by OBCFAs (C15:0; anteiso C15:0; iso C15:0; C17:0; iso C17:0; C17:1 and anteiso C17:0) in the bacteria membrane lipids24. Thus, the higher concentration of OBCFAs might be the result of the difference in the rumen bacterial populations induced by variation in the dietary carbohydrate, that is, a higher concentration of cellulolytic bacteria in relation to amylolytic bacteria, due to the high neutral detergent fiber (NDF) content in the diet with 0% cactus forage silage. It is also known that amylolytic bacteria produce more linear odd chain and anteiso FAs than iso FAs, whereas cellulolytic bacteria produce more iso FAs28,32. As the Tifton 85 grass hay-based diet had the highest neutral detergent fiber corrected for ash and protein (NDFap) and starch content (highest % of ground corn), the meat of those animals had higher concentrations of anteiso C15:0 and anteiso C13:0 compared to animals fed diets with the inclusion of cactus pear silage, also influencing the total sum of branched chain fatty acids.Although levels of intermittent water supply have generated punctual changes in tricosanoic acid (C23:0) SFA, the same was not observed for MUFA and PUFA, due to changes in the rumen environment, promoted by water restrictions, which were not sufficient to circumstantially modify biohydrogenation, resulting in similarities in concentrations of unsaturated fatty acids in goat meat.The animals subjected to 24 h of intermittent water supply (IWS) presented the highest concentration of C23:0 in relation to other treatments, which is interesting because it is involved in the synthesis of ceramide and reduces the risk of diabetes in humans33.The cactus pear has high non-fibrous carbohydrate (NFC) content (mainly pectin), having 59.5% high and medium rumen degradation carbohydrates which provide a higher production rate and removal of short-chain fatty acids and changes in rumen bacterial populations34. The inclusion of CPS resulted in a higher passage rate of digesta, affected biohydrogenation, and resulted in the escape of intermediate fatty acids isomers that are absorbed in the small intestine. Consequently, there was changing composition of fatty acids in the muscle of these animals, with a significant effect being observed only in the cis-13 C18:1. Furthermore, diets with high proportions of cactus pear silage (CPS), such as 42% CPS diet, can decrease ruminal pH and affect the final stages of biohydrogenation, resulting in the escape of intermediate fatty acids isomers, that are absorbed in the small intestine, which can explain the similarity of the C20:1 in 42% CPS diet from the Tifton hay-based diet, with differences between goat meat from 21% CPS diet and Tifton hay-based diet.Oleic acid (c9-C18:1) was the MUFA with the highest participation in the lipid profile of goat meat, which is interesting because it has a hypocholesterolemic effect, being a desirable fatty acid (DFA) for not reducing the serum high density lipoproteins (HDL) levels and thus prevent cardiovascular disease by reducing LDL levels35. The high concentrations of c9-C18:1 in ruminant meat come from the food intake, the effect of biohydrogenation, and mainly of the high activity of Δ9-desaturase, necessary for animal biosynthesis through desaturation of C18:0 to c9-C18:127. This fatty acid in the lipid profile of red meat varies between 30 and 43%36, confirming that the meat in the present study had a good concentration of this fatty acid.Much of unsaturated fatty acids, which have 18 carbons or 16 carbons, are largely converted to C18:0 and C16:0 through biohydrogenation, and when this process is not 100% completed, in addition to the PUFA that pass through this process intact, some product intermediates are formed, reaching the duodenum and are absorbed by the animal, in which significant amounts of cis and trans-monounsaturated, such as vaccenic fatty acid (t11-C18:1), reach the duodenum and are absorbed, later composing the muscle tissue22.The literature indicates that the precursor of conjugated linoleic acid (CLA) in the meat of animals is trans vaccenic acid (t11-C18:1), so the enzyme ∆9-desaturase, besides acting in the conversion of stearic into oleic fatty acid, also converts the trans-vaccenic acid to its corresponding CLA isomer, c9t11-C18:236. This pathway is more expressive in the mammary gland, and as the concentration of vaccenic acid (t11-C18:1) was not different, the concentration of CLA was not affected by the supply of silage and intermittent water supply, in the same way, that there are also no differences in the activity of ∆9-desaturase. Nevertheless, it is worth noting that in the human adipose tissue there is also the presence of ∆9-desaturase, and therefore, increased intake of vaccenic fatty acid could have the same beneficial effects associated with the intake of CLA, where the dietary vaccenic fatty acid shows 19–30% conversion rate37.Tifton hay is a natural source of n-3 fatty acids, mainly C18:3 n-3 with up to 20% participation in the lipid profile2, allowing a certain part of these PUFAs to be absorbed and increased in the tissue muscle, with 10 to 30% PUFAs in the diet generally escaping from biohydrogenation.Linoleic fatty acid (c9c12 C18:2) and α-linolenic acid (C18:3 n-3) are essential fatty acids for humans, that serve as precursors of the n-3 and n-6 pathways, distinct families, but synthesized by some of the same enzymes (∆4-desaturase, ∆5-desaturase, and ∆6-desaturase)25. Arachidonic fatty acid (C20:4 n-6) comes from elongation and desaturation of linoleic acid, where its concentrations, even close to that of its precursor, may indicate that there was a high activity of ∆6-desaturase (desaturation to γ-linolenic), elongase (elongation of γ-linolenic to dihomo-gamma-linolenic) and ∆5-desaturase. This fatty acid was influenced by the diets, presenting lower concentrations in the meat of animals fed the 42% cactus pear silage when compared to the Tifton hay diet (0% cactus pear silage).A higher concentration of long-chain PUFA n-3, docosahexaenoic (C22:6 n-3), was observed in the muscle of animals fed on Tifton hay. This was probably due to the high concentration of C18:3 n-3, precursor of C22:6 n-3, that the hay presents in relation to the cactus pear silage.The ratios and proportions of fatty acids are used to determine nutritional and nutraceutical values of the product or diet, and mainly, to indicate the cholesterolemic potential4. It is interesting that the n-6/n-3 ratio is low due to the pro-inflammatory properties of n-6; it is recommended to decrease its intake to assist in disease prevention38, while n-3 fatty acids are anti-inflammatory, antithrombotic, antiarrhythmic and reduce blood lipids, with vasodilating properties, being interesting that they present a higher proportion24. n-6 fatty acids tend to have a higher percentage in meat, and this directly influences the formation of n-3 isomers, since linoleic acid, when in excess, can reduce the synthesis of linolenic acid metabolites. The percentage of FA in one group can interfere with the metabolism of the other, reducing its incorporation into tissue lipids and altering its general biological effects38. Therefore, it is not recommended that the n-6/n-3 ratio be kept above 5 or 639, demonstrating that the averages of the current research remained acceptable.In relation to atherogenicity index (AI) and thrombogenicity index (TI), Ulbricht and Southgate39 proposed that sheep meat should have values of up to 1.0 and 1.58, respectively, and the lower the values for these indices in the lipid fraction, the greater the prevention of early stages of cardiovascular diseases. In the present study, the general averages observed were 0.29 for the AI, and 0.81 for the TI, although there were no significant differences, all treatments are within the recommended range, despite having been used as comparative standard to sheep, due to the absence of the proposed standard for goat meat.The h:H ratio did not differ for diets and water supply levels, but had an average of 1.90, below the reference value for meat products, which is 2.0. Values above 2.0 are recommended and favorable40, as it indicates a higher proportion of hypocholesterolemic fatty acids, that are beneficial to human health.The ∆9-desaturase enzyme that acts on both the mammary gland and adipose tissue, responsible for the transformation of SFA into unsaturated fatty acids (UFA), as well as in the endogenous conversion of CLA37 did not differ between treatments. On the other hand, the elongase showed less activity. Probably there was a greater “de novo” synthesis which resulted in a greater accumulation of palmitic fatty acid, and a reduction in the activity of the elongase enzyme.The crossbred goats demonstrated to present efficient mechanisms for adapting to water restrictions, especially when receiving feed with higher water content, such as cactus pear silage, being able to replace Tifton hay with 42% cactus pear silage in the diet for goats in confinement without negatively affecting the carcass traits and meat quality. Because, although these animals have shown some differences in the indices of tenderness and juiciness of their meats, however, all presented values of juiciness and tenderness compatible with meat extremely appreciated by the consumer market, and even goat meat showing some fatty acids with different concentrations induced by the supply of silage and water intermittence, the final lipid profile was appropriate to the health of consumers, observed by the absence of differences in the total concentrations of PUFA and in the main nutraceutical parameters (DFA, n-6/n-3; h:H; AI and TI).These results are relevant, indicating that goat feedlots in regions with low water availability may adopt strategies of lesser demand for drinking water and considerable concentrations of cactus pear silage in the diet, can reduce production costs without considerably affecting the product to be marketed, and therefore, provide higher profitability of the system. More

Disturbance detectionWe used a well-established spectral-temporal segmentation method, Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr), to detect disturbances within the Google Earth Engine (GEE) cloud-computing platform57,58. The core of the LandTrendr is to extract a set of disturbance-related metrics by breaking pixel-level annual time-series spectral trajectories into linear features using Landsat observations. The LandTrendr has been widely used for change detection in various forest settings, and details about the algorithms can be found in previous publications57. Here we briefly described the key steps in generating the year and type of disturbances in China’s forests using the LandTrendr within the GEE platform. The overall analytic flow can be found in Supplementary Fig. 10.First, we generated annual spectrally consistent time-series data by using all available, good quality (cloud cover ≤ 20) Tier 1 Landsat 5 (Thematic Mapper), Landsat 7 (Enhanced Thematic Mapper Plus), and Landsat 8 (Operational Land Imager) images acquired during the peak growing seasons (June 1—September 30) from 1986 to 2020. The peak growing seasons were selected to exclude compounding influences from ice, snow, and soil, and to maximize the spectral changes after forest disturbances. To tackle the spectral inconsistency among Landsat sensors, we harmonized spectral values via linear transformations according to band-respective coefficients presented in59. Clouds, cloud shadows, snow, and water were masked out using the Fmask algorithm60. The annual band composites at 30-meter spatial resolution during 1986–2020 were computed using the Medoid method61.Secondly, we ran the LandTrendr using five spectral indices, including two spectral bands (shortwave infrared I and II that were B5 and B7), tasseled cap wetness (TCW), normalized burn ratio (NBR), and normalized difference vegetation index. These five indices were effective indictors to represent vegetation greenness and structures, and were commonly used for detecting changes in forest disturbance and recovery62. For each spectral index, the LandTrendr produced a set of parameters to describe a possible disturbance event at the pixel level, including spectral values at pre-disturbance level (preval), magnitude of change (mag), duration (dur) and rate of change (rate), and the signal-to-noise ratio (dsnr) (n = 5). Using these five spectral indices, we generated a stack of disturbance-related parameter layers (n = 25, 5 spectral indices × 5 parameters), which were later used to detect and classify disturbances using machine learning models derived from reference data (described below).Disturbance classificationReference dataHigh-quality consistent reference data is key to train and classify disturbance types. To do so, we generated a total of 31225 reference points using a hierarchical approach. We first generated a large number of potential disturbance points using forest loss data from 2001 to 20203. Then we separated fire disturbances from non-fire disturbances by overlaying MODIS burned area (BA) with potential disturbance points following the procedure used by63. Specifically, fire disturbances were determined if the MODIS BA data coincided with the Landsat-derived forest loss for the fire year and 2 years postfire (i.e., t + 0, t + 1, t + 2) to account for delayed post-fire tree mortality. Following this step, we derived points as potential disturbances that consisted of fires and non-fire disturbances (including forest conversion to other land use types and silvicultural practices at various intensities). We also generated roughly the same number of points that experienced no disturbances (e.g., persistent forests), which were determined by selecting pixels with very few changes in spectral indices. These reference points, including fire, non-fire disturbances, and persistent forests, were then used to sample the time-series spectral data from 1986 to 2020. Finally, time-series spectral data from each reference point were visually checked to make sure they accurately represented disturbance events. This process resulted into a total of 31225 reference data points, including 2356 fire disturbance points, 13,242 non-fire disturbance points, and 15,627 no disturbance points (persistent forests) (Supplementary Fig. 2).Random forest classificationWe used machine learning modeling to classify each pixel into fire disturbance, non-fire disturbance, or no disturbance. The reference data points were used to sample the LandTrendr-derived disturbance-related parameter layers described above, which resulted into a dataset consisting of disturbance types. We divided the dataset into 70% of training data, and 30% as validation data. Using the training data, a Random Forest (RF) model was trained to classify each reference point into fire, non-fire disturbance, or no disturbance. Our RF approach showed that short-wave infrared (SWIR)-based moisture indices (e.g., B7, TCW) were strong predictors for detecting forest disturbances (Supplementary Fig. 11) likely because of their sensitivity to vegetation water content and canopy structure64. Finally, we applied the trained RF model to the full classification stack to consistently map the disturbance types from 1986 to 2020 across China’s forests, assuming that the spectral trajectories derived from reference data period 2001–2020 can be extrapolated to the whole mapping period 1986–2020. However, note that our approach was meant to detect relatively acuate and discrete disturbances that caused canopy opening, rather than subtle changes of forest structure or composition resulted from low intensive silvicultural practices and chronic disturbances.Year of disturbanceWe used the LandTrendr to determine the year of disturbance as the onset of magnitude of spectral change. Since we ran LandTrendr on five spectral indices, there were five possible years of disturbance for each pixel. Thus, we determined the year of disturbance using the median value from at least three different indices (i.e., NDVI, NBR, TCW, B5, B7). In this way, we only kept pixels that were detected as disturbances using at least three indices, thus reducing commission errors. The year with the greatest spectral changes generated by the LandTrendr often had an accuracy within 3 years11. A confidence level was also assigned to each disturbed pixel based on numbers of indices which showed possible disturbance events. Specially, low, medium, and high confidence were assigned if the disturbance was detected by three, four, or five spectral indices, respectively.ValidationsWe validated the disturbance map at the pixel and national levels. At the pixel level, we validated the final map using the validation sub-sample described in the previous section. We derived a confusion matrix to report user’s and producer’s accuracy (Supplementary Table 1) as the main accuracy assessment metrics. At the national level, we compared forest disturbance detected in this study to available existing dataset. Specifically, we compared the area of forest fire disturbance between our study and the national fire records during 2003–2009 (Supplementary Fig. 5). We compared the disturbance rates between our study and Landsat-derived global forest cover changes from 2001 to 20193 (Supplementary Fig. 4).Post-processingWe applied a series of spatial filters to minimize the unrealistic outliers from two potential sources of uncertainty, including speckle in time-series spectral trajectories or misregistration among images. This may lead to individual pixel or small patches including only a few pixels, which were (a) detected as disturbances, thus increasing the commission errors, or (b) not detected as disturbances, while their surrounding pixels were mostly disturbed, thereby increasing the omission errors. To address the issue (a), we removed all single-pixel disturbance patches through setting the minimum mapping unit as two 30 × 30 m2 pixels (0.18 ha). To address the issue (b), we applied a 3 by 3 moving window to fill holes through assigning the year of disturbance based on the years in the surrounding pixels. Finally, we smoothed the year of disturbance by assigning the center pixel using majority rules from surrounding pixels within the 3 by 3 windows, thus accounting for artefacts associated with uncertainties in the correct identification of the disturbance year.Characterizing disturbance regimes and their trendsWe characterized the disturbance regime using five indicators within each 0.5° grid cell (n = 1946) across China’s forests based on annual forest disturbance maps generated from the previous step. Within each grid cell, we calculated (1) total annually disturbed forest area (km2 yr−1), (2) percentage of forest disturbed annually (% yr−1), as annual disturbed forest area divided by the total forested area, (3) disturbance size (ha), as the number of disturbed pixels for each individual patch using an eight-neighbor rule, (4) disturbance frequency (# of patches per 1000 km2 forested area each year), as the number of disturbance patches per year divided by the total forested area, (5) disturbance severity (ΔNDVI = NDVIt−1 − NDVIt+1), as magnitude of NDVI change 1 year before and 1 year after disturbance, obtained from the LandTrendr analysis. We used (1) and (2) to characterize the disturbance rate, and (3)–(5) to describe the patch characteristics. The (2) and (4) were normalized by forest area within each grid cell, thus making them comparable among grid cells. For (3)–(5), we only calculated the patch size >0.45 ha (five 30 × 30-m2 pixels), because patches  TC2000), and the expansion of forested area from 1986 to 2000 (e.g., TC1986  20% following Liu et al., (2019). We should note that our study area did not include the newly afforested area after 2000. All analyses were performed within the forest mask, thus excluding the potential confounding factors from other land cover types. The description of TC1986 and TC2000 can be found in3,32. More