dataesr/scientific-paragraphs-categorization huggingface dataset:
Cross-comparison of read depths among the EP samples revealed important intraspecies variations in copy number for single-and multicopy genes (see Materials and Methods and see Table  S7  at GitLab [https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental _adaptation_reveals_strain_specific_differences/]). Plotting the gene coverage values for the three L. infantum isolates, the three L. donovani isolates, or the two L. major isolates together with the L. tropica sample, resulted in strong, confined signals at the center of the ternary plots that correspond to genes with equal copy number and thus a 33% distribution across the three axes (Fig. 2 , left panels). Compared to the different reference genomes, we observed important, strain-specific differences in gene copy number that are visualized on these plots by shifts of the signals out of the center.
It is interesting to speculate that this heterogeneity reflects individualized solutions driving fitness gains in vitro. While differences in culture conditions certainly account for some of the observed karyotypic variability, the comparison of two closely related Spanish L. infantum isolates, Linf_LLM45 and Linf_LLM56, reveals a culture-independent component implicated in genomic adaptation. Both isolates were adapted to culture at the same time under the same conditions, yet they showed important differences in karyotype dynamics, with only Linf_LLM56 demonstrating changes in somy levels at EP ϩ 3 (Fig. 3 ; see Table  S4  at GitLab [https://gitlab. pasteur. fr/gbussott/ Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific _differences/]). These strains are genotypically identical (zymodeme MON-1) (see Table S1 at GitLab) and are genetically closely related, with an average nucleotide identity of over 99.
Significantly, the SNV frequency profiles for EP and EP ϩ 3 were largely identical, ruling out the possibility that adaptation occurs through selection of subpopulations that would cause important shifts in SNV frequency distribution (data not shown). Together our results document the highly dynamic nature of karyotype management in Leishmania during environmental adaptation that is likely governed by complex interactions between external cues and intrinsic genetic differences. Dynamic variations in gene copy number during de novo culture adaptation. Plotting of genome-wide sequencing coverage of EP ϩ 3 against EP for all annotated genes resulted in a largely diagonal distribution, suggesting that there are no major CNVs between the two different passages (Fig. 4A ; see Fig. S4  in the supplemental material and see Table  S9  at GitLab [https://gitlab.
MATERIALS AND METHODS Leishmania parasite isolation and culture. Ten Leishmania strains belonging to the L. tropica, L. major, and L. donovani complexes of eight different geographical areas were isolated from infected patients, dogs, or hamsters (see Table  S1  at GitLab [https://gitlab. pasteur. fr/gbussott/Leishmania_genome _dynamics_during_environmental_adaptation_reveals_strain_specific_differences/]). Some strains were cryopreserved in liquid nitrogen prior to culture adaptation until used for this study (see Table  S1  at GitLab). Leishmania isolates were first stabilized in vitro in media that were optimized in the various LeiSHield partner laboratories (stabilization medium; see Table  S2  at GitLab), prior to expansion in classical RPMI culture medium for a defined number of passages (expansion medium).
To determine reproducibility of in vitro genome evolution, duplicate EP ϩ 3 samples (EP ϩ 3. 1 and EP ϩ 3. 2) were generated for the Linf_ZK27, Lmj_1948, Lmj_A445, Ldo_BPK26, and Ltr_16 strains (Fig. S1 ). Culture conditions and time in culture for the 25 samples are detailed in Table  S2  at GitLab. Nucleic acid extraction, sample preparation, and sequencing analysis. Procedures for DNA sample preparation and quality control were standardized using common protocols. Briefly, DNA extraction was performed using DNeasy blood and tissue kits from Qiagen according to the manufacturer's instructions. Nucleic acid concentrations were measured with Qubit, and the DNA quality was evaluated on agarose gel. Between 2 and 5 g of DNA was used for sequencing. The following samples showed small DNA amounts and were thus PCR amplified before sequencing: Ldo_LTB_EP (5 cycles), Ldo_LTB_EP ϩ 3 (5 cycles), Linf_02A_EP (10 cycles), Linf_02A_EP ϩ 3 (5 cycles).
Eventually, PCR and optical duplicates were labeled with Picard MarkDuplicates [version 1. 94(1484)] (https://broadinstitute. github. io/picard/) using the option "VALIDATION_STRINGENCYϭLENIENT. " While the reads were aligned against full assemblies, including unsorted contigs, just the canonical 36 chromosomes were considered for downstream analyses of ploidy estimation and copy number alterations. This filter was necessary because of the high content of repetitive elements and the absence of comparable and high-quality annotations in the contigs. Given that the L. tropica reference genome is still unfinished, the sample Ltr_16 was aligned against the L. major Friedlin genome. Overall, starting from a total of 1,011,803,806 short reads, 952,093,114 were successfully aligned to the respective reference genomes (see Table  S3  at GitLab [https://gitlab.
pasteur. fr/gbussott/Leishmania_genome_dynamics _during_environmental_adaptation_reveals_strain_specific_differences/]). Picard CollectAlignment-SummaryMetrics was used to estimate sequencing and mapping statistics. Comparative genome analysis. Whole-genome sequencing data from the EP Leishmania isolates were processed with Trimmomatic (version 0. 35) (34) to remove low-quality bases (options "LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15") and adapter contaminations (option "ILLUMINACLIP," with values 2:30:12:1:true). Reads that were shorter than 36 bases after filtering were discarded (option "MINLEN:36"). The trimmed reads were assembled with SPAdes (35) (version 3. 7. 0) with the option "careful. " The resulting contigs were used to estimate the average nucleotide identity (ANI) with the dnadiff part of the MUMmer system (version 3.
23)  (36). The analysis included the reference genomes of L. donovani, L. infantum, and L. major that were retrieved from the Sanger database (described above) and reference genomes of L. braziliensis, L. mexicana, and L. panamensis that were retrieved from ENSEMBL Protists release 29  (37). The ANI values were converted to a matrix of distances, which in turn were used for principal-component analysis (PCA) and hierarchical clustering (R hclust function [https://www. r-project. org/]). Chromosome sequencing coverage. For each read alignment file, Samtools view (version 1. 3) and BEDTools genomecov (version 2. 25. 0) (38)  were used to measure the sequencing depth of each nucleotide. Samtools was run with options "-q 50 -F 1028" to discard reads with a low map quality score or potential duplicates, while BEDTools genomecov was run with options "-d -split.
Both the ANOVA P values and the chromosome somy comparisons are reported in Table  S4  at GitLab (https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain _specific_differences/). Gene sequencing coverage. Samtools view (version 1. 3) and BEDTools coverage (version 2. 25. 0) were used to measure the mean sequencing depth of every annotated gene and were run, respectively, with options "-q 50 -F 1028" and "-d -split. " Possible intragenic gap regions were excluded from the calculation of the mean. Then the mean coverage of each gene was normalized by the median coverage of its chromosome. To account for GC content sequencing bias, the coverage values were corrected using a LOESS regression with a 5-fold cross validation to optimize the model span parameter.
Genes supported by reads with a mean mapping quality (MAPQ) score of Ͻ50 were filtered. To enable CNV analysis of gene arrays and genes sharing high sequence identity, we clustered the nucleotide sequences of the annotated genes into groups with cd-hit (version 4. 6)  (39). We used the length difference cutoff option "-s 0. 9. " Then we realigned the clusters with MAFFT (40) and used T-Coffee seq_reformat  (41)  to select a representative gene per cluster (RefGene) showing the highest average sequence similarity to the other cluster members. If two genes had the same average similarity, then the shortest was chosen. We used bwa to build a database containing only the sequences of RefGene, adding Ϯ50 bp of the 5= and 3= ends to ease the read alignment and the quantification of small RefGenes.
We realigned EP samples against this database using bwa mem with the option "-M. " We then quantified the RefGene mean coverage (without considering the Ϯ50-bp extension) with Samtools view and BEDTools coverage using the options "-F 1028" and "-d -split," respectively. Values were normalized by the median coverage of the RefGene's chromosome. Gene groups composed by members located on different chromosomes were negligible and discarded. Genome binning. The reference genomes were divided into contiguous windows of a fixed length, and the sequencing coverage of each window was evaluated and compared across different samples. A window length of 300 bases was used for the shown scatter plots assessing genome-wide CNVs. Both the mean sequencing coverage normalized by the median chromosome coverage and the mean read MAPQ value were computed.
To account for GC content sequencing bias, the coverage values were corrected using a LOESS regression with a 5-fold cross validation to optimize the model span parameter. The windows with a MAPQ score below 50 in either EP or EP ϩ 3. 1 were discarded. Poorly supported windows with a median or mean sequencing depth smaller than one-tenth of the median chromosome coverage both in EP and EP ϩ 3. 1 were also discarded. The windows with an EP ϩ 3/EP coverage ratio outside the axis limits were placed on the edge (value of 3). In the genome browser tracks, the repeat elements and low-complexity regions were predicted with RepeatMasker (version 4. 0. 6) (RepeatModeler software; AFA Smit and R Hubley, RepeatModeler Open-1. 0, 2008 -2015 [http://www. repeatmasker. org]) using the options "-e crossmatch -gff -xsmall -s" in combination with Repbase  (42)  to identify Leishmaniaspecific and ancestral repeats.
" The output was filtered to retain the positions with just one alternate allele with a minimum frequency of 0. 9 and a minimum mean mapping quality of 20 for the reads supporting the reference or the alternative allele. SNVs mapping inside homopolymers (i. e. , simple repeats of the same nucleotide) were filtered using a more stringent parameter, requiring at least 20 reads supporting the variant. The homopolymers were defined as the DNA region spanning Ϯ5 bases from the SNV, with over 40% of identical nucleotides. We discarded SNVs with sequencing coverage above or below 4 median absolute deviations (MADs). The predicted SNVs are reported in Table  S5  at GitLab (https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific _differences/).
Synteny analysis. The synteny analysis was performed with SyntView (46), a software package originally designed to compare microbial genomes. The tool was adapted to browse interactively the genome of four Leishmania reference genomes and explore their syntenic relation: L. infantum JPCM5, L. donovani PBQ7IC8, L. major Friedlin, and L. donovani BPK282A1. This new tool hosting Leishmania syntenic data is publicly available at http://genopole. pasteur. fr/SynTView/flash/Leishmania/ SynWebLinfantum. html. Supplementary table availability. All supplemental tables are publicly available at GitLab at https:// gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals _strain_specific_differences/. Accession number(s). Reads were deposited in the Sequence Read Archive database (SRA) database  (47)  and are publicly available under accession no.
the supplemental material). Linf_ZK27 and Ldo_LTB displayed the most stable karyotypes between EP and EP ϩ 3. As judged by read depth values corresponding to integer or intermediate chromosome copy number values, full or mosaic aneuploidy was observed for four chromosomes in Linf_ZK27 (chromosomes 6, 9, 31, and 35) and six chromosomes in Ldo_LTB (chromosomes13, 15, 20, 23, 31, and 33), which were established at EP and maintained at EP ϩ 3 (Fig. 3; see TableS4at GitLab [https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific_differences/]). All other isolates showed higher intrastrain karyotype instability with both gain and loss of chromosomes observed between EP and EP ϩ 3. Linf_02A represented the most extreme example showing significant changes in read depth for 21 chromosomes (Fig.
ACKNOWLEDGMENTS This study was supported by a seeding grant from the  Institut Pasteur International Department  to the  LeiSHield Consortium , the  EU   FP7  (Kaladrug-R, contract  222895 ), the  Belgian Science Policy Office (TRIT ,  P7/41 ), the  Department of Economy, Science and Innovation in Flanders (ITM-SOFIB) , and the  Flemish Fund for Scientific Research  ( G. 0. B81. 12  to J. C. D. and F. D. ).
34. 1680 4. 07 6. 09 3. 99 Putative amastin-like surface protein LinJ. 36. 0190 3. 1 5. 62 7. 22 Elongation factor 2 a For full data, see TableS7at GitLab (https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific_differences/). bAsterisks indicate genes shown in the right panel of Fig. 2. c ncRNA, noncoding RNA.
Gene_id b Normalized mean read depth of: the supplemental material). Linf_ZK27 and Ldo_LTB displayed the most stable karyotypes between EP and EP ϩ 3. As judged by read depth values corresponding to integer or intermediate chromosome copy number values, full or mosaic aneuploidy was observed for four chromosomes in Linf_ZK27 (chromosomes 6, 9, 31, and 35) and six chromosomes in Ldo_LTB (chromosomes13, 15, 20, 23, 31, and 33), which were established at EP and maintained at EP ϩ 3 (Fig. 3; see TableS4at GitLab [https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific_differences/]). All other isolates showed higher intrastrain karyotype instability with both gain and loss of chromosomes observed between EP and EP ϩ 3. Linf_02A represented the most extreme example showing significant changes in read depth for 21 chromosomes (Fig.
06 0. 89 44 Hypothetical protein a For full data, see TableS7at GitLab (https://gitlab. pasteur. fr/gbussott/Leishmania_genome_dynamics_during_environmental_adaptation_reveals_strain_specific_differences/). PacBio L. donovani LDBPK assembly and annotations were downloaded on 02/05/2017 (ftp://ftp. sanger. ac. uk/pub/project/pathogens/Leishmania/donovani/LdBPKPAC2016beta). bAsterisks indicate genes shown in the right panel of Fig. 2. in.
, 2007; Damm et al. , 2017)  where high spatial variability can be observed  (Durand et al. , 2009b; François et al. , 2014). Third, snow conditions were often simulated using simplified degree day modelling approaches  (Dawson and Scott, 2013; Hendrikx and Hreinsson, 2012)  and neglected the differences between natural snow and groomed or machine made snow properties  (Pons et al. , 2015; Gilaberte-Búrdalo et al. , 2017). The present work aims at producing snow reliability investigations of a wide range of ski resorts in France (Alps and Pyrenees), Spain and Andorra under past and future conditions using state-of-the-art snowpack modelling. We accounted for snow grooming and snowmaking  (Spandre et al. , 2016b)  and used adjusted and downscaled climate projections from the EURO-CORDEX dataset  (Verfaillie et al.
Definition of relevant elevations of ski resorts All data on the geographical location and technical data on ski resorts were extracted from the "BD Stations" database  (François et al. , 2014; Spandre et al. , 2015). Ski Lifts installation and operation in France are supervised by the STRMTG ("Services Techniques de Remontées Mécaniques et Transports Guidés"). The STRMTG is a public service in charge of the safety control of French ski lifts providing authorizations for ski lift operations. The STRMTG manages a database (CAIRN: CAtalogue Informatisé des Remontées Mécaniques Nationales) dedicated to ski lifts which includes technical characteristics of each ski lift such as the ski lift power. The ski lift power is an indicator of the size of a ski lift, defined as the product of the elevation difference between the bottom and the top of a ski lift (in km) and its flow of persons per hour (pers h -1 ), expressed in pers km h -1.
Ski lifts infrastructures in France have a total ski lift power of 977 000 pers km h -1 , 94% of which are included in the present study  (Appendix B). These data are completed with geographical information from the database BDTOPO (25 m of resolution) developed by the French Geographical Institute (IGN, "Institut Géographique National"). Ski resorts operated by a single company are aggregated in the present study  (François et al. , 2014)  e. g. Stations du Mercantour (Isola 2000, Auron, La Colmiane). The following elevations were used to be compared with the snow reliability line: • The ski lifts mean elevation is defined as the average of top and bottom elevations of each ski lift weighted by its ski lift power, being simply referred to as the mean elevation of the ski resort  (François et al.
, 2014; Falk and Vanat, 2016). • The village elevation of a ski resort is defined as the average elevation of all urban areas within 300 m of the bottom of a ski lift, weighted by their surface area  (Breiling and Charamza, 1999). Based on the OpenStreetMap project (http://www. opensnowmap. org/) we estimated the main features for the Spanish and Andorran ski resorts (village elevation, mean ski lift elevation and ski lift power). We extracted the ski slopes from the open source database for all ski resorts in France, Spain and Andorra and estimated these parameters from the area resulting of a positive/negative buffer and intersection to define a spatial representation for every ski resort. Linear models were fitted based on French ski resorts and applied to Spanish and Andorran ski resorts to provide theses features (Figure  1 ).
Climate forcing data The meteorological system SAFRAN  (Durand et al. , 1993)  provides meteorological data (temperature, precipitations, etc. ) for mountain areas of an approximate 1000 km 2 surface referred to as "massif", covering French Alps and Pyrenees, including Spanish and Andorran Pyrenees (Figure  2 ). Within each massif, the meteorological conditions are supposed to be homogeneous and to depend only on the altitude (300 m altitudinal step) with a time resolution of 1 h. SAFRAN forcing data are available for the 1958 -2015 period  (Durand et al. , 2009a; Maris et al. , 2009; Durand et al. , 2012). Computations of snow conditions over the reference period using SAFRAN forcing data are further referred to as "SAFRAN" and can be considered as the reference observational dataset. This study uses the EURO-CORDEX dataset  (Jacob et al.
Spatial variability A significant spatial variability of the snow reliability line can be observed for the reference period  (1986 -2005) , the median elevation of the 70% groomed snow reliability ranging between 1750 m. a. s. l in the Northern Alps, 2000 m. a. s. l in the French Pyrenees, 2250 m. a. s. l in the Southern Alps and up to 2300 m. a. s. l in the Spanish and Andorran Pyrenees (HIST, Figure  3 ). Although a bias can be observed, the spatial variability is consistent between climate models computations over the reference period (HIST) and the reference dataset (SAFRAN). The 90% snow reliability using snowmaking is significantly lower than the 70% groomed snow reliability line (Figure  3 ). The gain thanks to snowmaking on the median elevation ranges between 700 m in the French Pyrenees, 900 m in the Spanish and Andorran Pyrenees, 1000 m in the Northern Alps and up to 1200 m in the Southern Alps.
Reliability of ski resorts Figures  3  and 4  and Table  2  show a bias between the SAFRAN reference dataset and results derived from climate models (HIST) for the reference period. We therefore focus our analysis on the comparison of snow conditions computed by climate models for the reference and future periods. Based on climate models, ski lifts infrastructures were reliable during the reference period  (1986 -2005) , either with natural snow conditions (50% in categories 1, 2 and 3 all together, Table  2 , HIST) or technically (49% in categories 4 and 5 all together). Natural snow conditions in larger ski resorts were more reliable than in the smaller ones with 44 resorts representing 50% of the ski lift power being natural snow reliable and 129 ski resorts also representing 49% of the ski lift power being only technically reliable (Table  2 ).
, 2015). This situation is likely to be more frequent under future climate conditions resulting in increasingly optimistic results compared to the reference period. Additionally, all results computed based on the observational reference dataset and climate models exhibit differences in the reference period (Figures  3  and 4  and Table  2 ). Discrepancies may be due to potential biases of the multivariate distribution of the meteorological variables produced by the adjustment and downscaling method  (Verfaillie et al. , 2017). This could result in potential nonlinear effects due to multiple dependencies especially on temperature, relative humidity, precipitation and wind-speed. Beyond the modelling of snow conditions the main limitations pertain to the snow reliability line approach. Single points representations are considered on flat field i.
Acknowledgements. This study was funded by  Région Rhône-Alpes  (PhD grant of  Pierre Spandre ), and benefited from funding from the  French Ministry for Ecology (MTES)  to the  ADAMONT  project through GICC and  ONERC , from the  Interreg  project  POCTEFA/Clim 'Py, from the  IDEX Univ. Grenoble Alpes Cross Disciplinary Project  " Trajectories ", and from the  Conseil Départemental de l'Isère -Isère Tourisme. The  CGET -Comité de Massif Alpes  funded the creation of the BD Stations database. CNRM/CEN and Irstea are part of  LabEX OSUG@2020  ( ANR10 LABX56 ).
Functional Annotation of Phosphatase Genes Eight families of non-specific phosphatases were previously described in the literature, with an optimal activity in an alkaline (PhoA, PhoD, PhoK, PhoX) or acid (AcpA, PhoN, AphA, NSAPc) pH range. In order to identify their homologs in the genome of Rta, we first selected a set of reference sequences corresponding to bacterial proteins with available 3D structures, one per family (Table  1 ). In addition, we described each family by a specific set of pre-calculated sequence profiles (PSSMs, for Position-Specific Scoring Matrices). These PSSMs were retrieved from the NCBI Conserved Domain Database (version 3. 14)  (Marchler-Bauer et al. , 2015)  according to their similarity to the family reference sequence, using the RPS-Blast program, a standalone version of CD-search  (Marchler-Bauer and Bryant, 2004).
For each family of phosphatases, we removed the profile redundancy as follows: among the PSSMs whose significant hits on the given reference sequence overlap with a mutual coverage exceeding 70% of the length of the longest hit, only one profile (the one with the lowest E-value) was added to the dataset characterizing the family. The final sets of reference PSSMs are listed in Table  1. All translated coding sequences (CDS) of Rta, i. e. , the putative proteome as provided by the automatic structural annotation of the Rta genome  (De Luca et al. , 2011) , were then compared by RPS-Blast with all selected sets of reference PSSMs. In this functional annotation procedure, a given CDS was annotated as a member of a given phosphatase family when it shared at least one significant hit of each PSSM of the family reference dataset.