dataesr/scientific-paragraphs-categorization huggingface dataset:
p function. All rate estimates were calculated with three ε values (ε = 0/0. 5/0. 9), knowing that the extinction rate in swallowtails is usually low  34  (supported by the results of this study). Second, we relied on the state-dependent speciation and extinction (SSE) model, in which speciation and extinction rates are associated with phenotypic evolution of a trait along a phylogeny  143. In particular, we used the MuSSE 144 implemented in the R-package diversitree 0. 9-10 151 , which allows multiple character states to be studied. Larval host-plant data were taken from previous works  5, 18, 30, 34, 112, 113, 152. The following ten host-plant character states and corresponding ratios of sampled species in the tree of all known species for each character (sampling fractions) were used: 1 = Aristolochiaceae (110/152), 2 = Annonaceae (69/138), 3 = Lauraceae (33/39), 4 = Apiaceae (9/10), 5 = Rutaceae (119/163), 6 = Crassulaceae (19/19), 7 = Papaveraceae (44/44), 8 = Fabaceae (1/1), 9 = Zygophyllaceae (2/2), and 10 = Magnoliaceae (2/2).
After a burn-in of 500 steps, we estimated posterior density distribution for speciation, extinction and transition rates. There have been concerns about the power of SSE models to infer diversification dynamics from a distribution of species traits  [153] [154] [155]  , hence other birth-death models were used to corroborate the results obtained with SSE models. Third, to provide an independent assessment of the relationship between diversification rates and host specificity, we used the ML approach of Morlon et al. 145  implemented in the R-package RPANDA 1. 3  156. This is a birth-death method in which speciation and/or extinction rates may change continuously through time. This method has the advantage of not assuming a constant extinction rate over time (unlike BAMM  146  ), and allows clades to have declining diversity since extinction can exceed speciation, meaning that diversification rates can be negative  145.
We accounted for non-random incomplete taxon sampling using the implemented analytical correction; we set a sampling fraction per genus based on the known species diversity of each genus. Mixing and convergence among runs (ESS > 200 after 15% burn-in) were assessed with the R-package BAMMtools 2. 1  159  to estimate (1) the mean global rates of diversification through time, (2) the estimated number of rate shifts evaluating alternative diversification models comparing priors and posterior probabilities and (3) the clade-specific rates through time when a distinct macroevolutionary regime is identified. Fifth, BAMM has been criticized for incorrectly modelling rate shifts on extinct lineages, that is, unobserved (extinct or non-sampled) lineages inherit the ancestral diversification process and cannot experience subsequent diversification-rate shifts  158, 160.
To solve this, we used a Bayesian approach implemented in RevBayes 1. 0. 10 161 that models rate shifts consistently on extinct lineages by using the SSE framework  147, 158. Although there is no information of rate shifts for unobserved/ extinct lineages in a phylogeny including extant species only, these types of events must be accounted for in computing the likelihood. The number of rate categories is fixed in the analysis but RevBayes allows any number to be specified, thus allowing direct comparison of different macroevolutionary regimes. Finally, we evaluated the impact of abrupt changes in diversification using the Bayesian episodic birth-death model of CoMET 148 implemented in the R-package TESS 2. 1  162. These models allow detection of discrete changes in speciation and extinction rates concurrently affecting all lineages in a tree, and estimate changes in diversification rates at discrete points in time, but can also infer mass extinction events (sampling events in which the extant diversity is reduced by a fraction  163  ).
2. 0 171. Finally, CDS alignments were strongly cleaned from misaligned sequences (gene by gene) using HMMCleaner 1. 8  172. A last cleaning step was performed using trimAl 1. 2rev59  173  , which is designed to trim alignments for large-scale phylogenomic analyses. The resulting dataset comprised 6621 genes in at least four sampled species (median of 32% of missing data), which was used to reconstruct a robust phylogenomic tree of Papilionidae 46 (Supplementary Fig. 18 ). We used this genomic dataset of 45 species representing all genera in which the resulting genus-level swallowtail phylogenomic tree  46  accurately represents the evolutionary associations with host plants as estimated using the ancestral-state analyses applied to the species-level phylogeny  34  (Fig. 1  and Supplementary Figs.
4  and 5 ). We thus transferred the inference of ancestral host-plant shifts on the phylogenomic tree and selected the branches representing a host-plant shift and/or a shift of climate preference (in general from tropical to temperate conditions; Supplementary Fig. 10 ). We also selected branches with no change as negative controls  73. As a result, 14 branches are selected to measure the impact of a host-plant shift and 14 branches are selected as controls (Supplementary Fig. 18 ). Within these 14 branches with an ecological change, nine branches represent host-plant shifts only, and five branches correspond to shifts in both host plant and environment (from tropical to temperate conditions). To test the impact of these different changes on the genomes, two datasets were created, Dataset 1 and 2.
Given the low quality of the genomes of Allancastria cerisyi and Parnassius imperator, these two genomes were discarded for the downstream analyses. We first selected the genes from the 6621-gene dataset for each focal branch using three criteria: (1) the dataset is composed only of orthologous protein-coding genes (OrthoFinder 2. 2  171  ), (2) the species needed to accurately define the branch were available (i. e. crown node of the clade) and (3) for each branch, one species per tribe was available, and therefore include a different number of genes per branch. Thus, for Dataset 1, only the genes containing sequences for the species needed to generate all focal branches were selected. This stringent selection leads to Dataset 1, comprising only 520 genes but the same genes for all branches (no missing genes).
For Dataset 2, the genes were selected for each branch independently (i. e. for a given branch, a gene was selected if the sequence needed to generate that branch was present). This second selection leads to 1439 genes per branch on average among a total of 1533 genes, which were selected at least once for one branch. The genomic dataset is available in Figshare (see Data availability). We studied the ratio (ω) of dN/dS to find genes under positive selection  76, 174. The dN/dS ratio is traditionally used to estimate selective pressure from proteincoding sequences. If host-plant shifts have no effect on the selection of a given gene, we expect a dN/dS = 1 and the selective regime is considered neutral. However, if host-plant shifts result in positive selection on coding genes, the ratio increases such that dN/dS > 1.
5). The fit for including positive selection is tested using a likelihood ratio test comparing the null model with the alternative model with one degree of freedom  76, 176. If the alternative model is better suited to host-shift branches, it is more likely the gene was under positive selection during the host-plant shifts. For each gene and for each branch, both the null and alternative models using CodeML were implemented in PAML 4  177  with a fixed topology (as inferred with the phylogenomic dataset  46  ) and the nucleotide alignment of each gene. To test the robustness of the estimations, we used a false discovery rate test to control false positives  178. Finally, for each branch, we reported the number of genes under positive selection (i. e. for which the alternative model including the site class with dN/dS > 1 have a better likelihood) on the total gene number.
Acknowledgements This project has received funding from the  Marie Curie International Outgoing Fellow  under the European Union's  Seventh Framework Programme  (project  BIOMME , agreement no. 627684 ), a  PICS grant  from the  CNRS  (project  PASTA ), an ' Investissement d'Avenir ' grant from the  Agence Nationale de la Recherche  (project  CASMA, CEBA, ref. ANR-10-LABX-25-01 ) and the  European Research Council (ERC)  under the European Union's  Horizon 2020 research and innovation programme  (project  GAIA , agreement no. 851188 ) to F. L. C. ; a  Natural Sciences and Engineering Research Council of Canada (NSERC)   Discovery Grant  ( RGPIN-2018-04920 ) to F. A. H. S. ; and a  German Research Foundation  grant ( WA 2461/9-1 ) to S. W. We are grateful to  Sophie Dang ,  Troy Locke  and  Corey Davis  at the  Molecular Biology Service Unit of the University of Alberta  for their help, assistance, and advice on next-generation sequencing.
The availability of phosphorus (P) and nitrogen (N) constrains the ability of ecosystems to use resources such as light, water and carbon. In turn, nutrients impact the distribution of productivity, ecosystem carbon turnovers and their net exchange of CO 2 with the atmosphere in response to variation of environmental conditions in both space and time. In this study, we evaluated the performance of the global version of the land surface model ORCHIDEE-CNP (v1. 2), which explicitly simulates N and P biogeochemistry in terrestrial ecosystems coupled with carbon, water and energy transfers. We used data from remote sensing, ground-based measurement networks and ecological databases. Components of the N and P cycle at different levels of aggregation (from local to global) are in good agreement with datadriven estimates.
Second, models ignoring nutrients were often calibrated on available carbon data so that a new model with nutrients inevitably needs a parameter recalibration to reach similar performances as the same model without nutrients. Third, for evaluating a large-scale model resolving both nutrient and carbon biogeochemistry, one needs specific nutrient-related datasets, which are more scarce than classical biomass, productivity and soil carbon data used for benchmarking carbon-only models. The evaluation for N and P, together with carbon cycling in global LSMs, remains very limited  (Wang et al. , 2010; Goll et al. , 2012) , but recent advances in ground-based measurements, ecological datasets and process understanding have made a better evaluation of C, N and P models feasible. The available nutrient datasets have allowed for meta-analyses of site-level nutrient fertilization experiments (e.
g. ,  Yuan and Chen, 2015; Wright, 2019) , data-driven assimilation schemes to constrain nutrient budgets  (Wang et al. , 2018) , new knowledge about the critical P processes of sorption  (Helfenstein et al. , 2018 (Helfenstein et al. , , 2020) )  and phosphatase-mediated mineralization  (Sun et al. , 2020) , global datasets of leaf nutrient content  (Butler et al. , 2017) , and empirical constraints on the CO 2 fertilization effect on land carbon storage  (Terrer et al. , 2019; Liu et al. , 2019). In addition to direct comparison with nutrient datasets, it is also possible to diagnose emerging model responses in terms of ecosystem resource use efficiencies (RUE) and confront them with observations for identifying how ecosystems adjust and optimize nutrient, water, light and carbon resource availabilities  (Fernández-Martínez et al.
Model description ORCHIDEE-CNP simulates the exchange of greenhouse gases (i. e. , carbon dioxide, nitrous oxide), water and energy at the land surface and features a detailed representation of the root uptake of dissolved N and P, the allocation of N and P among tissues, and the N and P turnover in litter and soil organic matter  (Goll et al. , 2017a (Goll et al. , , 2018) )  (Fig. 1 ). In this study, we present the first global application of the model and an evaluation against global carbon and nutrient datasets. ORCHIDEE-CNP simulates the cycles of C, N and P, which are described in detail elsewhere  (Krinner et al. , 2005; Zaehle and Friend, 2010; Goll et al. , 2014 Goll et al. , , 2017a Goll et al. , , 2018)). We give a brief overview here. P enters the ecosystem through release from minerals into the soil solution, whereas N is biologically fixed from an ample reservoir of dinitrogen.
Meteorological data The model was forced by CRU-JRA-55 meteorological data provided at a spatial resolution of 0. 5 • × 0. 5 • and upscaled to a resolution of 2   (Harris et al. , 2014)  and the Japanese 55-year Reanalysis (JRA-55) data  (Kobayashi et al. , 2015) , covering the period 1901 to 2017. This climate dataset was provided by the TRENDY v6 model intercomparison project  (Le Quéré et al. , 2018).
Land cover The historic land cover change maps were based on the European Space Agency Climate Change Initiative (ESA-CCI) land cover data  (Bontemps et al. , 2013). To be used by global vegetation models like ORCHIDEE-CNP, ESA-CCI land cover data were aggregated to 2 • × 2 • and grouped into PFTs using the reclassification method from  Poulter et al. (2011 Poulter et al. ( , 2015)). The fraction of cropland and pasture in the PFT map was further constrained by the cropland area and the sum of pasture and rangeland area for the year 2010, respectively, in the History Database of the Global Environment land use dataset  (HYDE 3. 2; Klein Goldewijk et al. , 2017a, b) , which were also aggregated to 2 • × 2 •. The above processes produced a reference ESA-CCI-based PFT map for the year 2010. The land use changes derived from Land-Use Harmonization (LUH) v2 (http://luh.
Soil and lithology datasets ORCHIDEE-CNP v1. 2 is forced by information on soil texture, pH, bulk density and soil types  (Goll et al. , 2017a). We used a global gridded map of three soil texture classes from  Zobler (1986)  to derive soil-texture-specific parameters for soil water capacity, hydraulic conductivity and thermal conductivity. We used global gridded data on bulk density from the Harmonized World Soil Database (HWSD; FAO/IIASA/ISRIC/ISSCAS/JRC, 2012) and soil pH from the International Geosphere-Biosphere Programme Data Information System Soil Data  (Global Soil Data Task Group, 2000). Soil pH forcing maps are needed to simulate the dynamics of NH 3 and NH + 4 in soil in ORCHIDEE  (Zaehle and Friend, 2010). We used a global gridded map with the dominant soil orders (following the USDA Soil Taxonomy) at 1 • × 1 • resolution to derive soil-order-specific soil phosphorus sorption parameters  (Goll et al.
Nutrient management For croplands, we used yearly gridded mineral N and P fertilizer application data from  Lu and Tian (2017)  available for the period 1960 to 2017. This dataset is based on nationallevel data on crop-specific fertilizer application amounts from the International Fertilizer Industry Association (IFA) and the FAO. N and P mineral fertilization between 1900 and 1959 was linearly extrapolated assuming that fertilizer applications for 1900 are zero and that there were no N and P fertilizers applied before 1900. For pasture, we used global gridded datasets of mineral N fertilizer application rates from  Lu and Tian (2017) , developed by combining country-level statistics (FAO) and land use datasets (HYDE 3. 2)  (Xu et al. , 2019). For both cropland and pasture, N and P in mineral fertilizer were assumed to go directly into soil mineral pools, where all mineral N fertilizer was assumed to be in the form of ammonium nitrate, with half of N as ammonium (NH + 4 ) and half as nitrate (NO - 3 ).
Manure applications are also included as a model forcing, given their significant input contribution to agricultural soils. For cropland, we used gridded annual manure N application data for the period 1860-2014 from  Zhang et al. (2017)  compiled and downscaled based on country-specific annual live-stock population data from FAOSTAT. For the period before 1860, we assumed N and P deposition rates of the year 1860. For pasture, we used global gridded datasets of N manure application rates from  Lu and Tian (2017). The application of manure P in cropland and pasture was derived from manure N assuming a manure P : N ratio of 0. 2. This ratio is a weighted value by the amount of manure N applied to soil and derived from ruminants (14. 4 Tg N yr -1 ) and monogastric animals (10. 1 Tg N yr -1 ) from FAOSTAT for the year 2000 with P : N ratios of 0.
Evaluation We evaluated the performance of ORCHIDEE-CNP v1. 2 based on four major aspects (Fig. 1 ). Firstly, we evaluated the global C, N and P flows and storages. In the absence of robust spatially resolved estimates of N and P fluxes, we used the data-driven reconstruction of steady-state C, N and P fluxes on the biome level from the Global Observationbased Land-ecosystems Utilization Model of Carbon, Nitrogen and Phosphorus (GOLUM-CNP) v1. 0  (Wang et al. , 2018)  (Table  1 ). Secondly, we evaluated plant resource use efficiencies (RUEs) of light, water, C, N and P on global and biome scales. RUEs reflect how ecosystems adjust and adapt to the availability of nutrient, water, light, and carbon resources  (Fernández-Martínez et al. , 2014; Hodapp et al. , 2019). For this, we used estimates from site measurements and observation-based gridded datasets.
Thirdly, we evaluated the response of gross primary productivity (GPP) to elevated CO 2 to assess the response of plant productivity to changing resource availability (i. e. , CO 2 ) and historical perturbation C fluxes. For this, we used observationbased estimates  (Ehlers et al. , 2015; Campbell et al. , 2017). Fourthly, we evaluated large-scale patterns of vegetation and soil N : P ratios as well as the N and P openness and turnover rates on the ecosystem level to assess spatial variation in nutrient limitation and the underlying drivers. For this, we used estimates from GOLUM-CNP, site measurements and observation-based gridded datasets  (Kerkhoff et al. , 2005; McGroddy et al. , 2004; Reich and Oleksyn, 2004; Tipping et al. , 2016; Butler et al. , 2017; Wang et al. , 2018). Finally, we showed the implications of ORCHIDEE-CNP for C cycling by evaluating the spatiotemporal patterns of terrestrial C fluxes and pools of the two versions of ORCHIDEE.
Ecosystem productivity Different data-driven maps of NPP and GPP based on remote sensing and climate datasets were used (Table  1 ), thereby accounting for the uncertainty of each product and for the uncertainty from the spread between different products. Uncertainties of each NPP and GPP product were derived according to the original publications. We used a 20 % uncertainty of gridded GPP from the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Breathing Earth System Simulator (BESS) (Sect. S1C in the Supplement;  Turner et al. , 2006; Jiang and Ryu, 2016)  at a 2 • scale. This is a coarse extrapolation of uncertainty reported at the grid cell scale, since none of these products reported spatial error covariance information, allowing us to upscale this uncertainty at 2 • resolution.
Further, for some products, uncertainty was defined as the bias against local measurements  (Turner et al. , 2006)  and for others by using different climate input fields (Table  1 ). For multi-tree ensemble (MTE) GPP (Table  1 ), we used the spread (1σ standard deviation) from an ensemble of 25 members produced by different machine-learning methods  (Jung et al. , 2009). For MODIS NPP (Table  1 ), we used a 19 % uncertainty as assessed by  Turner et al. (2006). For BETHY NPP we do not have an uncertainty  (Tum et al. , 2016). For Global Inventory Modeling and Mapping Studies (GIMMS) NPP (Table  1 ), we used the variance of three sets of products (Table  1 ) based on different climate datasets  (Smith et al. , 2016). Two statistical indices were used to summarize the performance of ORCHIDEE and ORCHIDEE-CNP with respect to the interannual and seasonal variability of GPP and the interannual variability of net biome productivity (NBP) (Sect.
4. 6): the coefficient of determination (R 2 ) and relative mean square deviation (rMSE). The rMSE is defined as rMSE = n j =1 X model,j -X ref,j 2 n j =1 X ref,j -Xref,j 2. (1) X model and X ref are values from models (i. e. , ORCHIDEE and ORCHIDEE-CNP) and referenced datasets (i. e. , MTE and BESS; Sect. S1C), respectively, and Xref,j is the mean value across all years (for interannual variability evaluation) or all months (for seasonality evaluation).
The impact of such perturbations on the nutrient flows are analyzed in detail in Sects. S2 and S3 in the Supplement. In terms of C and nutrient storages, ORCHIDEE-CNP simulated comparable soil C, N and P storage (soil organic matter and litter) but higher vegetation C, N and P than GOLUM-CNP. Detailed comparisons for the spatial pattern of soil organic carbon (SOC) and forest aboveground C against observation-based datasets can be found in Figs. S2  and S3.
Resource use efficiencies Here we evaluate the resource use efficiencies of GPP for light (L), water (W), C, N and P defined by  where GPP is the annual gross primary productivity (g C m -2 yr -1 ), fAPAR the fraction of absorbed photosynthetically active radiation (%), PAR the annual photosynthetically active radiation (W m -2 yr -1 ), ET the annual evapotranspiration (mm m -2 yr -1 ), and F N and F P the total N uptake (g N m -2 yr -1 ) and P uptake by plants (g P m -2 yr -1 ), respectively. We calculated fAPAR in ORCHIDEE-CNP and ORCHIDEE as a function of leaf area index (LAI): fA-PAR = 1 -exp(-0. 5 • LAI)  (Ito et al. , 2004). LUE = GPP f APAR × PAR , (2) WUE = GPP ET , (3) CUE = GPP GPP , (4) NUE = GPP F N , (5) Compared to observed LUE (Sect. S1E in the Supplement), ORCHIDEE-CNP modeled median values at the biome level are generally lower but still within the ranges of uncertainties of observation-based datasets (Sect.
Compared to the three sets of atmospheric inversions (CAMS, JENA and CTracker), ORCHIDEE(-CNP) performs slightly worse than the mean of predictions from 16 land surface models from TRENDY ensembles (v6) (Fig. 10c ). ORCHIDEE-CNP shows a worse performance in interannual variability of NBP than ORCHIDEE when compared against inversion datasets at global scale and for the Northern Hemisphere. However, ORCHIDEE-CNP improved the performance of the interannual variability of NBP against inversion datasets relative to ORCHIDEE for tropical regions (higher R 2 and lower rMSE), with closer or even better fitness against inversion datasets than the mean value of TRENDY ensemble models (Fig. 10c ).
Uncertainties in evaluation datasets hamper a more detailed evaluation of the variations of WUE among biome types. We found that the inclusion of nutrient cycles improved the spatial variability in simulated CUE, but general biases remain (Fig. 3c ), and uncertainties in observation-based estimates are large. Improvements are mainly found in temperate biomes (TEDF, TECF and TEG), indicating that the allocation of GPP to respiration and biomass growth, which is controlled by nutrient availability, works reasonably well. ORCHIDEE-CNP underestimates CUE for tropical biomes (TRF and TRG) more strongly than ORCHIDEE, despite substantially reduced biases in NPP and GPP (Fig. S4d ). However, we should be cautious in drawing conclusions considering the large uncertainty in MODIS CUE  (He et al. , 2018).