dataesr/scientific-paragraphs-categorization huggingface dataset:
Orthogroups reconstruction To cross expression data among species we reconstructed orthogroups including the nine species for which RNASeq data during symbiosis are available as well as 14 representative species of each NFN clades and Arabidopsis thaliana as outgroup (Table  S5 ). Orthogroups reconstruction was performed with OrthoFinder v2. 5. 2 (10. 1186/s13059-019-1832-y) with the ultra-sensitive Diamond mode (-S diamond_ultra_sens option). The inferred tree by OrthoFinder was then checked for consistent reconstruction with known species tree and OrthoFinder reran based on the species tree using alignment and phylogeny method with mafft v7. 313 (10. 1093/molbev/mst010) and fasttree v2. 1. 10 (10. 1093/molbev/msp077) respectively to infer orthogroups.
Statistical analysis The different scripts used to cross orthogroups and DEGs to estimate differentially expressed orthogroups (DEOGs), reconstruct the ancestral transcriptomic symbiotic response, identify symbiosome accommodation proteins and all statistical related analysis are freely available on the following dedicated GitHub repository https://github. com/CyrilLibourel/Universal_nodulation_transcriptomic_response. In order to estimate the null distribution of the number of DEOGs for each node, we first determined to which evolutionary node each orthogroup belongs using a maximum parsimony ancestral state reconstruction for discrete traits  24  (i. e. presence/absence of species genes in orthogroup). For each species, we randomly select genes accordingly to (i) the number of DEGs estimated in the species and (ii) the node that belongs the orthogroup containing DEGs.
Acknowledgements The authors thank  Oswaldo Valdes-Lopez ,  Pascal Gamas  and the  ENSA consortium  for discussions and feedback on the manuscript, and  B. Perret  for the Lupinus albus seeds. We are grateful to the genotoul bioinformatics platform  Toulouse Occitanie  ( Bioinfo Genotoul , https://doi. org/10. 15454/1. 5572369328961167E12 )for providing computing resources. This study was supported by the  Fédération de Recherche Agrobiosciences ,  Interactions et Biodiversité , the  French National Research Institute for Agriculture, Food and Environment (INRAE, Plant Health and Environment Division) , the  "Laboratoires d'Excellence (LABEX) "  TULIP  ( ANR-10-LABX-41 )" and by the  "École Universitaire de Recherche (EUR) "  TULIP-GS  ( ANR-18-EURE-0019 ). J. K. , C. L. and P-M. D are supported by the project  Engineering Nitrogen Symbiosis for Africa (ENSA)  currently funded through a grant to the  University of Cambridge  by the  Bill & Melinda Gates Foundation  ( OPP1172165 ) and the  UK Foreign, Commonwealth and Development Office as Engineering Nitrogen Symbiosis for Africa  ( OPP1172165 ).
This project has received funding from the  European Research Council (ERC)  under the European Union's  Horizon 2020 research and innovation program  (grant agreement No  101001675 -ORIGINS)  to  P-M. D. PR  has received funding from the  EU  in the framework of the  Marie-Curie   FP7 COFUND People Program , through the award of an  AgreenSkills+ fellowship  (under grant agreement  n°609398 ) and from the  European Union 's  Horizon 2020 research and innovation program  under the  Marie Skłodowska-Curie grant  agreement  N° 845838.
By testing the separate effects of herbivory and N-deposition, we found that increases in herbivory shifted forest plant communities towards species with higher nutrient demands as inferred from increases in the community-weighted mean of species indicator values for nutrients (CWM-N) (β = 0. 43, σ = 0. 15, Fig. 4a  and Supplementary Table  16 ). Similarly, greater increases in herbivory Color of points corresponds to the magnitude of change in site herbivore pressure between the baseline survey and resurvey (Supplementary Data file 1). Total nitrogen deposition (wet and dry, reduced and oxidized) is calculated using the EMEP database for the year 2000 and displayed across a color gradient of light to dark blue representing lowest to highest values at a spatial resolution of 10 km. were associated with a higher proportion of non-native species (β = 0.
36, σ = 0. 17; Fig. 4d ). , d, f, h ). Conditional effects of herbivory are depicted at the 10th (348 kg/ha; left) and 90th (1010 kg/ha; right) percentile of cumulative N-deposition in the data. There are n = 52 independent resurvey sites for all models. Lines and ribbons represent the posterior mean line and the 95% credible interval. Dashed lines represent statistically unclear relationships. Rugs in figure bottom in b, d, f, h depict the marginal distribution of the predictor. Cumulative N-deposition is calculated between the baseline and resurvey year per site. See Supplementary Fig. 3  and Supplementary Table  26 -29 for model outputs of the effects of N-deposition alone. See Supplementary Fig. 4  for interaction effects on species richness change and exchange ratio. Source data are provided as a Source Data file.
Nevertheless, our observational approach has limitations. Our approach cannot account for all potential confounding variables (or test for higherorder interactions among them). For example, although we did not find a statistically clear effect of site productivity here (Supplementary Tables  40 41 42 43 ), it would be useful to explore whether the interaction effects we observed would hold across larger productivity gradients. Furthermore, our dataset reflected the natural dominance of red, roe and fallow deer and wild boar species in Europe (Supplementary Data file 1). These species can have different feeding preferences, physiologies, and biomass requirements to the larger feeders such as bison and moose 72 that were only present in a few sites. Our results are likely to be driven by these dominant species, making it essential to understand whether the observed relationships would hold should populations of other herbivores increase.
Database We compiled baseline vegetation survey and resurvey data from 52 sites with associated herbivory data distributed across 13 European countries in the temperate deciduous forest biome (Fig. 1 ; see www. forestreplot. ugent. be and ref. 43 for inclusion criteria; Supplementary Tables  1, 2 ). These sites occur in historically continuously forested natural and seminatural forests that have not experienced any substantial change in land use (i. e. , no stand-replacing disturbance) either prior to the baseline survey or between the two surveys (sensu ref. 76). Site areas ranged from 5. 5 ha to 2. 5 × 10 6 ha (median: 2300 ha). At each site, researchers surveyed species in the herb, shrub and tree layer across 10 to 190 permanent or quasi-permanent plots per site (median: 50; total: 2928).
Time intervals between the baseline survey and resurvey ranged from 10 to 64 years (median: 47. 5 years). We accounted for changes in species taxonomy between surveys and sites by harmonizing species names following GBIF's backbone taxonomy  77, 78. This prevented double-counting species or inferring inflated estimates of turnover. Altogether, our dataset contains 1257 species across all sites and time periods. Note, because we include new resurveys here and herbivore densities were not available for all surveys in the forest-REplot database, our data comprise a different set of sites than previous forestREplot publications (e. g. , 50 and 15% overlap with refs. 43,56).
Explanatory variables Herbivore pressure. We quantified ungulate herbivore pressure at the level of a study site for the baseline survey and resurvey time period based on expert assessment from each site's dataset custodian in the forestREplot network  18. Custodians provided density estimates using the best available information alongside expert knowledge of the site. Many of these densities have been used in previously published analyses (see also refs. 8, 15, 37, 43, 79). Density estimates incorporate one or more of the following sources; interviews (e. g. , with local foresters, site managers, hunters, and national park administration), published and unpublished local data records, extrapolation of local/regional hunting statistics and/or direct animal count surveys. Herbivore densities were iteratively checked and revised twice by each dataset custodian and internally reviewed by the entire consortium.
8, 18, 82 (see Supplementary Data file 1 for all herbivore data and description of ordinal scale). We then also measured the equivalents of basal metabolic rates by multiplying the mean body mass of a species by its ordinal value, and summing across species per site. Mean body mass per species was taken from the Phylacine database  83. As this indicator correlated highly with the above herbivore pressure index (Pearson's ρ = 0. 82; Supplementary Fig. 6a ) and given an extreme outlier of one site (ρ = 0. 92 when this outlier was removed, Supplementary Fig. 6b ), we used the herbivore pressure index above in all analyses to avoid leverage. Temporal change in this index was calculated as the difference between the last resurvey and the baseline survey values per study site  18  , with change values ranging from -8 to 18 (Fig.
1  and Supplementary Data file 1). N-deposition. We quantified total cumulative N-deposition using the EMEP database (https://emep. int/mscw/mscw_moddata. html), using cumulative wet and dry deposition of oxidized and reduced nitrogen  84. We calculated the cumulative N-deposition between the baseline year and the year of the resurvey based on the methods described in ref. 8. First, we calculated N deposition between 1900 and the year of the baseline survey (N t1 ); second, we quantified the cumulative N deposition between 1900 and the resurvey (N t2 ); and third, we calculated the difference, N t2 -N t1 , to quantify cumulative N deposition between surveys. Therefore, cumulative N-deposition per site will be influenced by the rate of deposition per year, as well as the length of the intercensal intervals (i.
e. , sites with lower yearly rates but long intervals may have similar values to sites with high yearly rates but shorter intervals). The values of cumulative N-deposition ranged from 130 to 1296 kg ha -1 (Supplementary Table  2  and Supplementary Fig. 7b ). Site productivity. To control for the potentially confounding influence of productivity on vegetation responses to herbivory  85, 86  , we obtained local environmental data from the EuMedClim database  87  on both the potential evapotranspiration and annual precipitation for each site averaged across the baseline and resurvey years. We then calculated the annual precipitation to potential evapotranspiration ratio (AP:PET) as a productivity proxy metric  88. The EuMedClim database ranges from 1901 to 2014, so any sites that had been resurveyed since 2014 were given the 2014 value as their resurvey value.
Species richness change and exchange ratio. Species richness change was calculated as the difference in the number of herb layer species at each site between the resurvey and baseline survey  8, 90. Herb layer species turnover was estimated using the richness-based species exchange ratio  91  calculated at the site level as E = (S imm + S ext )/ S tot , where S imm is the number of species gained at resurvey, S ext is the number of species lost at resurvey and S tot is the total number of unique species at baseline and resurvey. The range of change values for species richness and exchange ratios were -98 to 90 and 0. 24 to 0. 65, respectively. Species of conservation concern. We identified species of conservation concern using two criteria: (1) Species listed as threatened in national Red Lists based on a recent database synthesizing national Red Lists across Europe  92.
1%, resurvey site mean = 3. 8%). Non-native and nutrient-demanding species. We identified the nonnative species present at each site using the Global Register of Introduced and Invasive Species (GRIIS; http://www. griis. org) and its designations of which species are non-native in each country. We then calculated the percentage of non-native species per site per survey period and the difference over time (baseline site mean = 3. 0%, resurvey site mean = 4. 8%). We estimated shifts in species' N-demands using ecological indicator values (EIVs) compiled from ref. 94 (sci. muni. cz/botany/juice/ELLENB. TXT), filling data gaps with values from ref. 95. Coverage of N-numbers was 92% of species (1156 out of 1257 species). For each study site and survey period, we quantified the community-weighted mean N-number (CWM-N), weighted by baseline and resurvey occupancy per species per study site (i.
Acknowledgements J. S. was supported by the project:  TERRANOVA , the  European Landscape Learning Initiative , which has received funding from the  European Horizon   2020 research and innovation program  under the  Marie Sklodowska-Curie  grant agreement No  813904. H. M. P, I. R. S and N. F were funded by iDiv via the  German Research Foundation  ( DFG FZT 118 ). D. T. received funding from the  EU Horizon 2020   INFRAIA  project  eLTER-PLUS (871128 ). F. M. was supported by the  Slovak Research and Development Agency  (project  APVV-19-0319 ) and also project  VEGA1/0624/ 21. T. S. was supported by  National Research, Development and Innovation Fund  ( K-135252 ). M. C. , R. H. , M. K. , M. M. and O. V. were supported by the  Czech Academy of Sciences  (project  RVO   67985939 ).
Data availability The data generated in this study have been made available in the Supplementary Information. Source data for figures are provided with this paper. The underlying species composition data were available from forestreplot. ugent. be, but restrictions apply to the availability of these data, which were used under license for the current study and so are not publicly available. These data were, however, available from the authors upon request and with the permission of the forestREplot consortium. Source data are provided with this paper.
Data-sets We created a dataset consisting of the variables included in this study from 1995 -2018 in an 8daily time resolution and a 3 km spatial resolution over the Euro-CORDEX 36 domain. The time period was limited by the necessary overlap of different data products, and the temporal resolution by the lowest one among the data products. The spatial extent and resolution were chosen to match upcoming land-surface model outputs. Therefore, where needed, the variables were interpolated in time to daily mean across 8 days and in space to the 3 km Euro-CORDEX grid. We used GLASS12B02 V40 (http://www. glass. umd. edu/introduction. html)  32  GPP data, which was created with the revised EC-LUE model 33,  49  which, importantly, includes long-term effects of CO 2 , temperature (T) and VPD on GPP.
The ET data used in this work is from GLASS11B02 V41, calculated by a multimodel ensemble approach merging five process-based ET data-sets  32, 50. Importantly, we compared GPP and ET from FLUXNET2015 (https://fluxnet. org/data/fluxnet2015dataset/)  51  and ICOS Drought-2018 (https://www. icos-cp. eu/data-products/YVR0-4898 ) 52 station data-sets to corresponding grid cells of GLASS, ERA5-Land (https://cds. climate. copernicus. eu/cdsapp#!/dataset/reanalysis-era5-land)  53  and GLEAM3. 5a (https://www. gleam. eu/)  54  and opted to use GLASS ET for the analyses based on the outcome (Supplementary Figure  S6 , S11 and S12, Tables  S3  and S4 ). Moreover, the T r data is from GLEAM3. 5a and was also compared to upscaled sap flux data from the SAPFLUXNET 0. 1. 5 dataset (http://sapfluxnet. creaf.
cat/)  55  (Supplementary Fig. S6 , Table  S5 ). To determine SSI droughts we acquired SM data from ERA5-Land for two soil layers, 0-7 cm and 7-28 cm. P data to identify hydro-climates (with mean annual P, Supplementary Table  S1 ) and SPI droughts, as well as other meteorological variables used here, temperature (T), relative humidity (RH), shortwave incoming radiation (R s ) and wind speed (W), were taken from the COSMO-REA6 data-set (https://reanalysis. meteo. uni-bonn. de/?COSMO-REA6)  48.
Survey procedure and analysis All the participants have been invited to access an online survey between March 1 st and March 15 th , 2020. Invitations to participate to the survey were sent by e-mail (one invitation followed by 2 reminders). The survey was available through the software Microsoft Forms (License Office 365 A1 for Angers University). The design and validation of the survey was performed by all authors who were also 4 pedagogical referents in our institution. The survey is reported in Table  1. Five-item Likert scales were used for questions 1 to 20. Questions 1 to 17 assessed students' and teachers' opinions (Likert items: "strongly agree", "agree", "neutral", "disagree" and "strongly disagree") and questions 18 to 20 their satisfaction (Likert items: "very satisfied", satisfied, "neutral", "unsatisfied", "very unsatisfied").
Statistical analysis Statistical analysis was performed with the SPSS 15. 0 Software ® (IBM Corp. , Armonk, NY, USA) and Systat statistical software v13 (Systat Software, Inc. , San José, CA, USA). All data were expressed as means ± standard deviation. Qualitative and quantitative variables were compared using Chi-square and Mann-Whitney tests. Differences between SCT and PCC were searched for each subject and compared using a Wilcoxon test. Paired analysis testing was performed for each student. The Spearman rank correlation test was used to assess the correlation. Statistical significance was defined as a p < 0. 05.
The snow heat flux was computed using observations collected on a 1 m pole equipped with 18 T-type thermocouples deployed from ground level up to a height of 0. 82 m (see inset in Fig. 1 ). The thermocouples were unequally spaced, with separations ranging from 2. 5 to 10 cm. To avoid disturbing the snowpack in cases where the spacing was small, some thermocouples were deployed with a 90 • deviation from the main axis. Four TP08 heated needle probes (Hukseflux, the Netherlands) were deployed at heights of 7, 27, 47, and 67 cm above ground. These probes recorded temperature and effective thermal conductivity according to the measurement principle and data treatment detailed in  Morin et al. (2010)  and  Domine et al. (2015 Domine et al. ( , 2016)). In short, the needle is heated at a constant power (0.
In addition to the automatic measurements, snow field surveys were conducted in April 2018 and March 2019. In both cases, three snow pits were dug at different locations in the vicinity of the experimental setup to capture the spatial variability of the snow cover. Snow pit data from previous campaigns  (Domine et al. , 2015)  were also included in the analysis. In each snow pit, the stratigraphy was analyzed, and profiles of density, temperature, and thermal conductivity were collected. Thermal conductivity was measured with a portable instrument featuring a TP02 heated needle, while density was measured with a 100 cm 3 box cutter  (Conger and McClung, 2009)  and a field scale. Monitoring and maintenance of the instruments was done during field trips in winter and every summer. No field trips to the site were possible in 2020 due to the COVID-19 pandemic.
Turbulent heat flux A detailed explanation of the procedure for obtaining the turbulent heat fluxes from raw eddy covariance data is provided by  Lackner et al. (2021). In short, turbulence data were processed using EddyPro ® (version 7. 0. 3; LI-COR Biosciences, USA), a software package that computes fluxes from raw 10 Hz data, while accounting for several corrections, including the application of a double rotation on the raw data to align the coordinate system with the current snow surface. EddyPro also includes a thorough QA/QC procedure. A program called PyFluxPro  (Isaac et al. , 2017)  was also used to remove spikes and erroneous data that persisted despite the EddyPro ® processing. Data gaps were filled using ERA5 https://doi. org/10. 5194/tc-16-127-2022 The Cryosphere, 16, 127-142, 2022 (  Hersbach et al.
Forcing variables and model setup The ISBA-Crocus model was run in offline mode, meaning it was not coupled to an atmospheric model. Instead, the model was driven by local observations of the following meteorological variables: air temperature, specific humidity, wind speed, incoming shortwave and longwave radiation, atmospheric pressure, and (solid and liquid) precipitation rates. These observations have been collected at the study site since 2012, except for atmospheric pressure, for which data are available from June 2017 onwards. Data from ERA5 were used to fill this gap by correcting for any bias based on the period when both ERA5 and the observations were available. Note that the sensitivity of the model to atmospheric pressure is known to be very low. Shortwave and longwave downwelling radiation between mid-December 2019 and mid-February 2020 had to be replaced by ERA5 data due to a malfunction of the radiometer.
In this study, we noticed only small differences between these three options (< 1 W m -2 at the daily scale), whereas major differences were reported from Col de Porte, France, as well as at almost all the sites in the ESM-SnowMIP dataset  (Ménard et al. , 2019). This is probably be-cause our site is windy, which helps to keep Ri b low and away from the very stable conditions that are more difficult to handle. Nevertheless, our results objectively confirm the statement of  Menard et al. (2021)  that turbulent fluxes are by far the most uncertain component of the energy balance simulation in snow cover modeling. We therefore highly recommend future complementary evaluations of turbulent fluxes of snow cover models, as was done in this paper for various climate and environments. Meanwhile, further studies are underway at the site targeted in this study to evaluate optimal turbulent heat flux parametrizations, particularly for sub-daily timescales.
Comparison to previous modeling attempts Studies comparing the full simulated and observed energy balance of an Arctic snowpack are sparse. Westermann et al. (2016)  used CryoGrid 3 (https://github. com/CryoGrid, last access: 15 August 2021) to simulate the energy budget of a polygon tundra site in Siberia. They reported that the model satisfactorily reproduced the winter energy balance but also found that it underestimated all heat fluxes. Helgason and Pomeroy (2012)  used the SNTHERM model in the central Canadian Prairies. They found that the latent heat flux and dU/dt were well simulated, while the sensible heat flux was overestimated. However, they concluded that the sensible heat flux was undermeasured and that the modeled flux was more applicable.
Acknowledgements. We are grateful to  Denis Sarrazin  for technical support of the CEN weather stations and the community of  Umiujaq  for granting permission to conduct our research. Financial support. This study was funded by  Sentinel North , a  Canada First Research Excellence Fund , under the theme 1 project titled "Complex systems: structure, function and interrelationships in the North". This research has also been supported by the ( Discovery Grant ), the  French Polar Institute (IPEV)  (grant no. 1042 ), and the  Horizon 2020 program  of the  European Commission  (INTAROS; grant no. 727890 ). CNRM/CEN is part of  Labex OSUG@2020  ( ANR-10-LABX-0056 ). Review statement. This paper was edited by  Chris Derksen  and reviewed by three anonymous referees.
Code availability. The source files of SURFEX code are provided at the Git repository (http://git. umr-cnrm. fr/git/Surfex_Git2. git, last access: 20 June 2021, Masson et al. , 2013)  with several code management tools (history management, bug fixes, documentation, interface for technical support, etc. ). Registration is required; a description of the procedure is described at https://opensource. cnrm-game-meteo. fr/projects/ snowtools/wiki/Procedure_for_new_users (last access: 25 July 2021). EddyPro ® is available on the LI-COR website (https://www. licor. com/env/support/EddyPro/software. html, Fratini and Mauder, 2014)  together with the corresponding manual. PyFluxPro is provided in the Git repository of OzFlux (https: //github. com/OzFlux/PyFluxPro; last access: 25 July 2021). The different levels are described by Isaac et al.