dataesr/scientific-paragraphs-categorization huggingface dataset:
1). Subsequently, chimera detection was carried on in QIIME with BLAST  24  and OTUs were taxonomically annotated employing the Greengenes 13. 8 database  25. At this point, two samples did not satisfy the quality filters and were discarded. Thus, for the remaining 285 samples, a dataset containing 1,294 OTUs was obtained after filtering out singletons and OTUs representing less than 0. 005% of the total number of annotated reads  23. From this dataset, 33 OTUs had missing taxonomic ranks and were discarded. Finally, 1,261 OTUs in the 285 samples were considered for further analysis (Supplementary Table  S1 ). The 1,261 OTUs were grouped in 18 phyla and 101 genera through the tax_glom method within the phyloseq package  26  in R (www. r-project. org). Besides, genera that belonged to a higher taxonomy rank but lacked the genus information were merged and marked as unspecified (g__unsp).
The analyses of α and β-diversities in the 285 samples were carried on with the vegan R package  27  , and the non-metric multidimensional scaling (NMDS) plot was performed using phyloseq 26 and ggplot2  28. For the α-diversity study, the Shannon index was employed, whereas the β-diversity study was represented using the Whittaker index. Additionally, the dissimilarity between pairs of samples was estimated with the NMDS method using the Bray-Curtis dissimilarity  29. Host DNA extraction and SNP genotyping. Pig genomic DNA was extracted from the Longissimus dorsi muscle of all the 288 samples using the standard phenol-chloroform method  30. The DNA concentration and purity was measured with a ND-1000 spectrophotometer (NanoDrop) afterwards. A total of 288 pigs were genotyped with the GeneSeek Genomic Profiler (GGP) Porcine HD v1 (70 K) array (Illumina, San Diego, CA, USA) using the Infinium HD Assay Ultra protocol (Illumina).
Genotypes were obtained with the GenomeStudio software (2011. 1 version, Illumina) and filtered with the PLINK software  31  (1. 90b5 version). Further analyses were conducted using only SNPs that mapped in the Sscrofa11. 1 assembly, with a minor allele frequency (MAF) >5% and missing genotypes <5%, retaining a total of 45,508 SNPs.
Acknowledgements This work was funded by the  Fondo Europeo de Desarrollo Regional (FEDER)  and the  Spanish Ministerio de Economía y Competitividad  ( MINECO AGL2014-56369-C2  and  AGL2017-82641-R ). D. Crespo-Piazuelo was funded by a " Formació i Contractació de Personal Investigador Novell " ( FI-DGR) Ph. D grant  from the  Generalitat de Catalunya  ( ECO/1788/2014 ) and by the  PiGutNet COST Action  (www. pigutnet. eu)for a  Short Term Scientific Mission  at the  GABI laboratory (INRA, France)  under the supervision of J. Estellé. Contract of L. Migura-Garcia was supported by  INIA  and the  European Social Fund. L. Criado-Mesas  was funded with a  FPI grant  from the  AGL2014-56369-C2  project. M. Revilla  was also funded by a  FI-DGR  ( ECO/1639/2013 ). M. Ballester was financially supported by a  "Ramón y Cajal " contract ( RYC-2013-12573 ) from the  Spanish Ministry of Economy and Competitiveness.
Phenotypic Innate Lymphoid Cell Studies NK cells, invariant NKT cells (iNKT) and γδ T cell subsets were characterized from fresh whole blood through multiparametric flow cytometry (Gallios flow cytometer and Kaluza Version 3. 0 analysis software, Beckman-Coulter, Kendon, FL, USA). NK cells were identified as CD3 -CD56 +/-CD16 +/-within the gated lymphocyte population previously defined in the FSC/SSC dot plot. The same strategy was used to define the iNKT cells, which were identified using the monoclonal antibody 6B11, specific to the conserved CDR3 region of the canonical hallmark T cell receptor (TCR)-invariant chain Vα24Jα18  (40) (41) (42)  and defined as CD3 + CD56 +/-Vα24Jα18 + lymphocytes, as detailed in the Supplementary Figure  1A. In its turn, the γδ T cell subsets were identified as CD3 + CD56 + Panγδ + lymphocytes coexpressing or not the Vδ2 chain.
Briefly, the slope indicates the growth/degrowth rate of an indicator during the "window of linearity, " which is the highest linearity region of the curve. The maximum and minimum points of the second derivate curve were used as boundaries for this region. The statistical analysis was performed using the Stata software version 11. 0 (STATA Corp. College Station, TX, USA), and the GraphPad Prism version 6. 0 software (GraphPad Software, La Jolla, CA, USA) was used to generate the graphics. Sigmoidal fitting was performed using R version 3. 0. 3 (2014-03-06-"Warm Puppy").
Disseminated TB was diagnosed in 15. 1% of the TB/HIV patients, whereas no case was detected in the TB group. Among the 33 TB/HIV recruited participants, only 4 (12. 1%) developed symptoms consistent with IRIS after beginning cART, which was in agreement with the frequency of IRIS observed in our previous studies  (28, 45). The IRIS onset was observed from 1-3 weeks after cART initiation and was clinically diagnosed as defined above. In general, the patients presented worsening of the signs/symptoms related to TB, like fever and enlargement of the right or the left cervical/submandibular lymph nodes, which were either self-resolving, or the patients were treated with corticoid-based therapy, such as Prednisone (Supplementary Table  1 ). At D0, the IRIS and non-IRIS TB/HIV patients presented similar CD4 + T cell absolute counts (data not shown).
Analysis of NK Subsets in IRIS Patients The impacts of the NK subsets on the inflammatory process were assessed in the TB/HIV IRIS patients despite the low number of confirmed IRIS cases. Prior to TB treatment, the non-IRIS patients presented a lower frequency of NK cells than the HC (p = 0. 009), and a trend for a lower frequency compared to that of the IRIS patients (p = 0. 055). Conversely, the IRIS patients and HCs had more similar NK cell frequencies (Figure  4A ). Interestingly, three NK receptors (NKG2C, NKp80, and CD158a, Figures  4B-D ) were decreased in the non-IRIS patients but were expressed at similar levels between the IRIS patients and HC. At IRIS onset, no significant difference was found in the NK cell subsets between the patients with and without IRIS (data not shown). Despite the limited number of IRIS cases, these observations suggest a more preserved NK cell profile among the IRIS patients, similar to the HC.
Potential for NK Cell Degranulation and Cytokine Production The potential for NK cell degranulation and cytokine production was also investigated among the groups. NK cell degranulation was addressed through CD107a expression; the results are shown in Figure  5. A significantly higher percentage of NK cells expressing CD107a was observed in the TB/HIV patients at D0 than among the HIV, TB, and HC groups (Figure  6A ). However, when we analyzed the degranulation capacity of NK cells according to the development of IRIS, no difference was observed between the IRIS patients and the other groups. Moreover, no difference was observed in CD107a-expressing NK cells in the TB/HIV patients regardless of the TB clinical presentation (data not shown). We also analyzed CD107a expression on NK cells under stimulated and non-stimulated conditions.
ACKNOWLEDGMENTS We thank the patients who participated in the study and all staff from the blood collection sectors from the  National Institute of Infectious Diseases Evandro Chagas  and  Nova Iguaçu General Hospital. We also thank  Michelle das Neves  and  Priscilla Guimarães  from the CD4 count and viral load clinical services of the  Laboratory of AIDS and Molecular Immunology, IOC-Fiocruz. Finally, we thank  Thaize Quiroga , MSc from the  Flow Cytometry Sorting Core Facility ,  IOC-FIOCRUZ , and the  Flow Cytometry Core of the Central Laboratory  from the  National Institute of Infectious Diseases Evandro Chagas-FIOCRUZ.
During the cold seasons across the study period , the total number of counts was 526, with a median number of counts per year of 20 (between 8 and 54). Flock size was calculated as the mean of the maximum value for each month. The quarterly mean of daily precipitation and of daily temperature were calculated from homogenized (discontinuity-corrected) climate data provided by the Swiss National Basic Climatologic Network for the Grand Saint Bernard station (2,470 m above sea level) (see  Begert, Schlegel, & Kirchhofer, 2005  for homogenization method). The quarterly mean of snowpack was calculated for each year of the study based on ISBA-Crocus models  (Vionnet et al. , 2012)  (elevational range 1,200-1,500 m) for winter and spring only. Lastly, the "snow anomaly" was calculated as the residual of a linear regression between snow cover and precipitation in cold seasons, so that the negative and positive values represent proxies of rainy and snowy seasons, respectively.
| Capture-recapture data analysis We performed a multistate goodness-of-fit test using U-care  (Choquet, Lebreton, Gimenez, Reboulet, & Pradel, 2009)  on a simple dataset including data only on adults and specific seasons (winter and summer). These tests revealed very high over-dispersion and trap-dependency (Appendix S1), thus we chose to develop a more complex model in line with the structure of the dataset and the biology of the studied species. The dataset was made up of 3,115 winter captures or resightings and 3,323 spring resightings at the Le Tour site, and 2,722 resightings in summer sites, for a total of 9,160 captures/resightings. This dataset allowed us to fit a model with a relatively high number of parameters. A multi-event capture-recapture (CR) modeling framework  (Pradel, 2005)  accounting for dispersal between monitored and unmonitored sites was used to test predictions (see below).
g. ,  Cayuela et al. , 2014). The combination of age (three age classes), site (three monitored and one unmonitored site) and resighting heterogeneity classes (high and low resighting probability) resulted in 24 "living" states and one "dead" state (Figure  S4  in Appendix S1). The CR data were analyzed on a seasonal basis, that is, in quarters corresponding to the four seasons defined above. From its initial state at first capture, an individual could transit from state to state between quarter t and t + 1 according to four steps (Figure  S4  in Appendix S1). The first step was the probability of surviving or not from t to t + 1. Survival could depend on an individual's state, notably its age and the summer site it occupied, and could also vary between quarters or between years. The second step was the probability of changing age class, which was not estimated and was forced to occur during the first (juvenile to immature) and second (immature to adult) autumn-winter transition.
Vegetation and Soil The land cover properties are specified according to the 1-km resolution ECOCLIMAP-II database  (Faroux et al. , 2013). More than 500 land cover units are derived from the Corine Land Cover map for the year 2000 at 100-m resolution over Europe and from the Global Land Cover 2000 database elsewhere. Theses land cover units are aggregated to the model resolution into 12 subgrid land tiles (Table  1  and Figure  2a ) in order to account for land cover heterogeneities in both ISBA versions.
7-5. 0 μm) vegetation and vegetation-free albedos are then retrieved from the MODIS 2001 to 2010 period at 1-km resolution by using a Kalman Filter-based method. The total snow-free land surface albedo in a grid cell is computed via an arithmetic average between vegetation and vegetation-free albedos weighted by the fraction of vegetation. This fraction of vegetation in the grid cell is specified for each subgrid land tile (Table  1 ). The soil textural properties, such as clay and sand, are given by the Harmonized World Soil Database (HWSD) at a 1-km resolution (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012; http://webarchive. iiasa. ac. at/ Research/LUC/External-World-soil-database/HTML/). We also use the percentage of organic carbon f oc (%) and the bulk density of the mineral matter ρ m (kg/m 3 ) over two soil horizons (0-0.
3 and 0. 3-1 m) from HWSD to calculate the SOC (kg/m 2 ) content (Figure  2d ) and the soil carbon density profile ρ soc (kg/m 3 ; section 2. 2. 5) in each horizon as follows: SOC ¼ ρ s Δz s f oc =100 (17a) ρ s ¼ 100 f oc =α om ρ om þ 100-f oc =α om ρ m -1 (17b) where ρ s (kg/m 3 ) is the total soil bulk density  (Adams, 1973) , α om is the dimensionless van-Bemmelen value of 0. 58 generally used to convert organic matter content to organic carbon content, ρ om (kg/m 3 ) is the average bulk density of the organic matter set to 224 kg/m 3 , and Δz s (m) is the thickness of each HWSD soil horizon. Note that the SOC content has been set to 0 kg/m 2 over desert and glacier areas defined in the ECOCLIMAP-2 database.
Subgrid Topography The topographic index used in both ISBA versions to determine the surface runoff via the TOPMODEL approach is specified from the HYDRO1K database (https://lta. cr. usgs. gov/HYDRO1K) at 1-km resolution. This distribution is calculated via a three-parameter gamma function derived from the mean, standard deviation, and skewness of the HYDRO1K actual distribution  (Decharme et al. , 2006 (Decharme et al. , , 2013)). CTRIP parameterizations require topographic information at a global scale and a high resolution. We use the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010; https://topotools. cr. usgs. gov/ gmted_viewer/) at a 7. 5-arc-second (~250 m) spatial resolution  (Danielson & Gesch, 2011)  to derive (1) the river bed elevation needed to calculate the river elevation  (Vergnes & Decharme, 2012) ; (  2 ) the river bed slope needed by Manning's formula to compute the river flow velocity  (Decharme et al.
The river width (Figure  2f ) is calculated using two steps. First, we use the following empirical relationship between the mean river width Ω riv (m) and the mean annual discharges Q yr (m 3 /s) in each CTRIP grid cell: Ω riv ¼ max Ω min ; α W ×Q β W yr ( 21 ) where α W and β W are equal to 5. 41 and 0. 59, respectively, and Ω min (m) is the minimum river width in a grid cell set to 30 m. This relationship was developed by  Vergnes et al. (2014)  using a large set of 1,873 gauging stations across France with access to river width, river depth, and mean annual discharge measurements. The mean annual discharge in each CTRIP grid cell is estimated at the global scale using the GGHYDRO global runoff database (http://people. trentu. ca/~gcogley/glaciology/glglgghy. htm) from  Cogley (2003). Second, we use a combination of two sources of satellite estimates: the Global Width Database for Large Rivers (GWD-LR; http://hydro.
iis. u-tokyo. ac. jp/~yamadai/GWD-LR/) at 0. 5°resolution  (Yamazaki et al. , 2014)  and the Global Lakes and Wetland Database (GLWD; http://wp. geog. mcgill. ca/hydrolab/glwd/), which gives the lake fraction distribution at 0. 5°resolution  (Lehner & Döll, 2004). The final river width in CTRIP is then given by W riv ¼ max Ω riv ; W space À Á (22a) W space ¼ W GWD ∀W GWD >0 f lake A ctrip L r ∀W GWD ¼ 0 8 < : (22a) where f lake is the GLWD lake fraction in each grid cell and W GWD (m) the GWD-LR river width estimate. Finally, the river bankfull depth, h riv (m), which controls the beginning of the river flood, is determined as follows, using an empirical relationship between observations of river width and maximum river height also proposed by  Vergnes et al. (2014) : h riv ¼ α h ×W β h riv (23) where α h and β h are equal to 1.
, 2011 (Schneider et al. , , 2014)). Second, we use the Tier-2 Water Resources Re-analysis (WRR2) forcing from the E2O project, available at 0. 25°r esolution and re-gridded to 1°resolution using conservative remapping. This data set is directly based on the 3-hourly ERA-Interim reanalysis (https://www. ecmwf. int/en/forecasts/datasets/reanalysis-datasets/ era-interim) over the 1979-2014 period. In the E2O project, only the 2-m air temperature and humidity and the surface pressure have been corrected for differences between topography at 0. 25°resolution and ERA-Interim topography at ~70-km resolution. Monthly precipitations have been hybridized with observations using the Multi-Source Weighted-Ensemble Precipitation (MSWEP; http://www. gloh2o. org) data set  (Beck et al. , 2017). The MSWEP-based E2O precipitation rate is globally larger than GPCC-based PGF data (Figure  3 ).
Evaluation of the Global Water Budget 6. 1. Major Land Surface Reservoirs 6. 1. 1. Floodplains To evaluate floodplain inundation simulations in sfxcm6, we use satellite-derived inundation estimates over the 1993-2004 period from  Prigent et al. (2007). This global data set is derived from a complementary suite of satellite observations that quantifies the monthly variations of the distribution of surface water extent at ~25-km resolution  (Papa et al. , 2010; Prigent et al. , 2007). In this study, this data set is remapped at 0. 5°resolution and corrected from lakes, bogs, fens, and mire areas of the GLWD database, as well as from monthly irrigated areas using the Monthly Irrigated and Rainfed Crop Areas around the year 2000 product  (Portmann et al. , 2010). Details on the correction method can be found in  Decharme et al.
Here we compare sfxcm5 and sfxcm6 simulations with available global data sets such as the Global Land Evaporation Amsterdam Model (GLEAM; https://www. gleam. eu/) root zone soil moisture and the ERA-Interim/Land (ERA-I/Land) total soil moisture reanalysis (https://www. ecmwf. int/ en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim-land). GLEAM data are available over the 1980-2016 period at a monthly frequency and at 0. 5°spatial resolution. GLEAM estimates global root zone soil moisture and evapotranspiration from a simple land surface model based on satellite forcing and soil moisture assimilation  (Martens et al. , 2017). The depth of the root zone in GLEAM is a function of land cover type: 0. 1 m for bare soil, 1 m for low vegetation, and 2. 5 m for tall vegetation.
River Discharges River discharges allow the water budget to be assessed over large areas and are an excellent validation test for land surface models, considering that the precipitation forcing is as perfect as possible. To perform such a comparison, we gathered a large set of daily discharge measurements at 20,164 stations from 1900 to present. These in situ measurements were provided primarily by the Global Runoff Data Centre (GRDC; http://www. bafg. de/GRDC/EN/Home/homepage_node. html), along with the USGS (http://waterdata. usgs. gov/nwis/ sw) over the United States, the HyBAm over the Amazon basin (http://www. ore-hybam. org), the French Hydro database (http://www. eaufrance. fr) over France, and the streamflow time series of the Parana river at Rosario from  (Antico et al. , 2018). From this entire network, we selected 698 stations from 1979 to 2010 that meet four conditions: (1) only one station per TRIP grid cell; if several stations are located in the same grid cell, the one with a CTRIP drainage Journal of Advances in Modeling Earth Systems area closest to the observed drainage area is kept; (2) observed upstream drainage areas must correspond to at least two atmospheric forcing grid cells at 1°(larger 24,000 km 2) ; (3) stations where the error between observed and CTRIP drainage areas is larger than 15% are removed; and (4) daily measurements must cover at least 1,825 days (5 years) over the 1979-2010 period.
Only one SNV in SUPT20H: c. 73A>T (p. Lys25*), presented a complete segregation with RA in an extended pedigree with early-onset cases. In summary, we identified in this study a new variant associated with RA in SUPT20H gene. This gene belongs to several biological pathways like macro-autophagy and monocyte/macrophage differentiation, which contribute to RA pathogenesis. In addition, these results showed that analyzing rare variants using a family-based approach is a strategy that allows to identify RA risk loci, even with a small dataset.
In all the 43 genes, the variant with the highest score is the rare family-specific SNV that was selected previously. Thereafter, we will refer to this variant as the candidate variant. To validate the leading effect of the candidate variants, we performed again the burden test by including first, only the candidate SNVs, and then, all variants except these candidates (results in Table  3 ). All the 43 candidate SNVs, except SUSD5: c. 526C>T (pVAAST p-value = 0. 21) and SMYD5: c. 662G>A (pVAAST p-value = 0. 067), were significantly associated with RA (pVAAST p-value � 0. 037). Concerning the test without the candidate SNV, three genes in the first group (20%), SUSD5, SMYD5 and MNS1, were still significantly associated with RA (pVAAST p-value � 0. 036). So, several rare variants in our dataset contributed to the association of these genes with RA.
Validation in extended pedigrees We selected the top 10 RA candidate variants, with the strongest association with RA, for resequencing to confirm their family specificity and their co-segregation with RA (Table  4 ). We analyzed further 8 out of 10 variants for which we were able to produce sequences and, we validated the familial specificity for all of them except PCCA. Four variants (in SUPT20H, SLC9A6, TIMM44 and NEK1 genes) were shared by all affected members in the families who carried them. However, only one candidate variant showed a complete segregation with RA. This variant is a heterozygous non-sense SNV, c. 73A>T (p. Lys25 � ), located in SUPT20H gene which introduces a premature codon stop at the beginning of the gene (in the fourth exon) (Fig  3 ). This variant, not reported to date, has been deposited in ClinVar database prior to publication.
Indeed, different studies reviewed by Dai and Hu demonstrated the role of autophagy in RA  [28]. In addition to SUPT20H gene, we identified 42 genes significantly associated with RA (pvalue < 0. 05). We evaluated their association with the test implemented in pVAAST. This software combines a rare variant association test (RVAT) with a linkage score, thus it takes into account for the familial structure of our data. It allowed us to provide all the variants called in our samples without any a priori filtering based on biological knowledge, avoiding loss of information. Indeed, the RVAT incorporates a score evaluating for each variant the amino acid substitution impact and phylogenetic conservation information, allowing it to be robust in the presence of neutral and common variants  [12]. Indeed, depending on the chosen test, variants with opposite effects can lead to a decreasing power for the association test  [29].
So, as expected, the score of the association test with all identified variants (familial data and control data) and with the candidate variant only were identical (see Table  1 ). We can highlight two limitations in this association analysis. First, we added as controls the variants of 98 individuals extracted from the 1000 genomes project database. And, although we chose individuals from the same population as our pedigrees (European) and sequenced on the same platform (Illumina), two aspects could have introduced bias: the exome capturing kit and the variant calling workflow. The first aspect can introduce variants specific to one of the dataset because not targeted by the kit used for the other one. To partially overcome this issue, we filtered out variants from the 1000 genomes dataset present outside the boundaries of our targeted exomes.
But, due to the wide range of capturing kits used to produce the data presents in the 1000 genomes database, we did not filter the variants according to the regions targeted by the 1000 genomes project. For the second aspect, it has previously been shown that despite a few variant specific to a variant calling software, the concordance between the called variants is high (92% between HaplotypeCaller  [31] , Samtools  [32]  and FreeBayes  [33] )  [34]. In addition, the PCA, performed on all the individuals included in the association test (Fig  2 ), shows that the observed genetic variability is mostly related to the intra-family specificity, not to the choice of the 98 CEU controls, comforting our choice of controls. Second, individuals of the discovery set have different clinical characteristics (such as sex, age at diagnostic, and ACPA status).