dataesr/scientific-paragraphs-categorization huggingface dataset:
Phylogenetic analyses In order to perform a phylogenetic analysis of the genus Yersinia, a representative subset was built from the 3878 Yersinia genomes available in our database. First, this set of isolates was partitioned into 30 clusters based on a UPGMA clustering from the pairwise dissimilarities between allelic profiles. Second, as these 30 clusters were of varying size n (n=1 to 1849 isolates), each was partitioned into 7 log 10 n UPGMA sub-clusters, and the isolate with the smallest number of missing alleles was selected within each subcluster, leading to a final set of 236 isolates representative of the genus Yersinia. For each of the 500 genes of the cgMLST scheme, the corresponding allele sequences were translated, and a multiple sequence alignment was performed using mafft v7.
407  [38]. The 500 multiple sequence alignments were concatenated into a unique supermatrix of 139 050 amino acid characters that was used to infer a maximumlikelihood phylogenetic tree using iq-tree v1. 6. 3  [39]  with the evolutionary model JTT+F+R5 selected by minimizing the BIC criterion  [40]. Following a similar approach, datasets within Y. enterocolitica and Y. pseudotuberculosis species led to two supermatrices of characters for 246 (+28 outgroup) and 294 (+20 outgroup) representative isolates, respectively, that were phylogenetically analysed with the evolutionary models JTT+F+R3 and JTT+F+R2, respectively.
BIGSdb A database was created for the genus Yersinia in the Institut Pasteur MLST and whole-genome MLST resource (https:// bigsdb. pasteur. fr), which uses the BIGSdb software tool, designed to store and analyse sequence data for bacterial isolates  [29, 30]. All de novo assembled genomes were uploaded into the isolates database, and the reference alleles of the 500 cgMLST loci were defined in the linked database of reference sequences and profiles ('seqdef '). Using BIGSdb functionalities, a scan of the genome sequence was performed for each isolate using parameters (min 80 % identity, min 80 % alignment, blastn word size of 20 nt) to check for the presence of each core gene and to determine its allele number. Allelic profiles identifiers (cgST) were defined for genomes with 50 or fewer missing alleles.
pestis and Y. wautersii), were individualized as distinct groups. Therefore, genomes of our dataset were classified into 28 groups in total. Note that this classification is not a taxonomic proposal, but is used instead as an operational classification. The correspondence between phenotypic and genome-based phylogenetic classifications is summarized in Table  2  and provided for individual isolates in Table  S1. We hereafter summarize these observations. within each group were all >96. 6 %. In contrast, the maximum ANI value observed between the members of these species and their closest relatives was 94. 7 % (Table  S1 ). These observations are concordant with the species status of these nine Yersinia members. Isolates identified as Y. mollaretii based on phenotypic characterization fell into two main lineages, one of them being further subdivided into two sublineages (Fig.
9%). Second, two additional lineages, which we call NEW 3 and NEW 4 (Fig. 1 ), were made up by isolates phenotypically characterized as Y. enterocolitica biotype 1A. These lineages comprised one and ten isolates, respectively, and were closely related to Y. enterocolitica. ANI values between the two lineages were 94. 7 % on average (>99. 3 % within NEW 4), and were 93. 9 and 94. 2 %, respectively, with Y. enterocolitica (Table  S2 ). Altogether, these results show that lineages NEW 1 to NEW 4 may represent four additional Yersinia species. Of note, NEW 1 was recently described as the new species Yersinia hibernica  [10] : the ANI between its type strain (accession numbers of its chromosome and plasmid are CP032487 and CP032488) and strain IP37048 (NEW 1) is 99. 9 %.
To evaluate this method for strain identification, 1843 Yersinia isolates received prospectively at the YNRL between 2016 and mid-2017 were analysed by the automatic cgMLST taxonomic assignment procedure. Details of the isolates used for this evaluation are presented in Tables  3  and 4. For 1814 (98. 4 %) isolates, genotypic assignment was consistent with phenotypic characterization (Table  3 ). Taxonomic assignments showed that the validation dataset comprised 669 non-pathogenic Yersinia isolates belonging to nine different species, 1116 pathogenic Y. enterocolitica isolates and 27 Y. pseudotuberculosis isolates. Only 29 (1. 6 %) isolates received a genotypic assignment that was discordant with phenotypic identification. Based on phenotypic characterization, these isolates belong mainly  to two species: Y.
The phenotypic Vitek GN card or MS (MALDI-TOF) represent alternative identification methods, but are not fully reliable for Yersinia species  [52, 53]. In addition, phenotypic identification at the infra-specific level and, thus, prediction of a strain's pathogenic potential, is not possible. For instance, neither species of the Y. pseudotuberculosis complex nor biotype 1A nor 1B strains of Y. enterocolitica can be distinguished using MALDI-TOF MS  [11, 17, 20]. Moreover, reference database updates would be needed for these methods to allow identification of novel taxa. We implemented a genus-wide cgMLST genotyping system based on 500 highly conserved protein-encoding genes that were assembled with no nucleotide ambiguities, had no paralogue and presented limited size variation. These stringent criteria were used to maximize genotyping reliability and, therefore, led to a subset of genes much smaller than the actual core genome of Yersinia.
Acknowledgements We thank  Sylvie Brémont  for performing phenotypic characterization of isolates, and  Andreea Alexandru ,  Maud Vanpeene  and  Vincent Enouf  (sequencing core facility,  Institut Pasteur)  for set-up of the DNA preparation conditions and for genomic sequencing. This work used the computational and storage services (TARS cluster) provided by the  IT department at the Institut Pasteur, Paris. We acknowledge the continuous support of  Keith Jolley  ( Oxford University ) for the development of the BIGSdb web application.
Data bibliography 1. Savin C, Martin L, Bouchier C, Filali S, Chenau J et al. GenBank BioProjects PRJEB3982, PRJEB3984, PRJEB3985, PRJEB3986, PRJEB3987, PRJEB3988, PRJEB3989, PRJEB3990, PRJEB3991 and PRJEB3992 (2014). 2. Reuter S, Connor TR, Barquist L, Walker D, Feltwell T et al. ENA studies PRJEB2116 and PRJEB2117 (2014). 3. Nguyen SV, Muthappa DM, Hurley D, Donoghue O, McCabe E et al. GenBank accession numbers CP032487 and CP032488 (2019).
For these simulations, we assumed P = 0. 1, E = 0. 05 and M in (0. 1, 0. 5, 1). In each generation, we recorded the number of trimorphic demes, the number of demes that had lost only the S-, M-and Lmorphs and the number of demes that had fixed the S-, M-or Lmorph. The global frequencies of the different morphs were also recorded each generation. The model was implemented in C, using the GNU Scientific Library for all functions related to random sampling, and is freely available at https://github. com/ popgenomics/triStyle.
Acknowledgements We thank  Spencer Barrett ,  Stephen Wright  and two anonymous reviewers for helpful comments on the manuscript. We particularly thank the organisers of the special issue for their invitation to contribute this article. JRP acknowledges  Spencer Barrett  for his scientific inspiration over the years, not least in prompting thoughts about tristyly in metapopulations when he was a postdoctoral placement in the  Barrett laboratory. Numerical results presented in this paper were carried out using the cluster from the  HPC M eso-centre of Lille University. JRP is supported by funding from the  Swiss National Science Foundation , and the  University of Lausanne.
Additionally, 8 mice were run on a 4 week version of the overlapping condition, with 3 weeks of each 25 trials (5 days, Mo-Fri) and a final trial on Wednesday of week 4. The second week's Monday (trial 26) as well as the final trial (trial 76) were run with the same configuration as the previous sample trial to function as 3-day and 5-day test trials. All animals were included in the analysis. Data acquisition. We developed an in-house program for training and scoring. The Object Scorer reads previously prepared training sheets with the object and locations for each trial of each animal, presents this information at the beginning of each trial to the experimenter (see  Fig 1) , and automatically extracts exploration times from the manually scored videos. Therefore, the operator cannot keep track of which animal is in which condition and which is the stable versus moved object for each trial, and he or she can be considered blind.
Source code for the Object Scorer software is available at https://github. com/MemDynLab/Score. Statistical analysis. The discrimination index (DI) used to assess memory performance was calculated as the difference in time exploring the novel object location and stable location divided by the total exploration time. This results in a score ranging from -1 (preference for the stable location) to +1 (preference for the moving object location). A score of 0 indicates no preference for either object location. Total exploration time and DI were assessed with repeated-measure ANOVAs with factors condition and trial (6 trials) in rats. Additionally, in rats, the DI of the test trial was assessed with a repeated-measure ANOVA for condition. Due to the different training schedule in mice (5 trials per day for 4 days but only 1 test trial on day 5), the sample trials were separately assessed by including the factors condition, trial, and day across the 20 sample trials in mice in an repeated-measure ANOVA.
The miniscope was attached to the baseplate, and proper focus of the field of view was established by using the nVista acquisition software (Inscopix, Inc). The animal was then placed in an empty cage for 10 min before the start of training. This cage also provided as a resting cage during the intertrial intervals. Recordings were taken from each 5-min sample trial and during each intertrial interval for another 5 min. To reduce risk of photobleaching, the blue LED was switched off for at least 5 min in between recording sessions. After the animal had completed the training day, the protocol was repeated for the next animal. In vivo electrophysiology in the medial prefrontal cortex and hippocampus in rats. Lister Hooded rats of 2 months of age (Charles River) were group housed with a 12-hr:12-hr light:dark cycle, and water and food were provided ad libitum.
Once the target area was reached, the recording sessions started for all the conditions of object location task previously mentioned (stable, overlapping, and random). Using the Open Ephys recording system combined with Bonsai software, electrophysiological recording took place during all training and test trials in the arena, as well as during 45 min before training (presleep), 45 min intertrial intervals, and 3 hrs after training in the sleep box. The head stage plugging prior to recordings did not cause any stress on animals due to the well handling and "Weeto" cereals animal fondness.
Model fitting. The aim here was to find for each individual subject, rat or mouse, the values of the model M's parameter set θ (here α and β) that best fit the data by maximizing the log-likelihood LL of the trial-by-trial sequence of proportions of time spent exploring object o 1 relative to object o 2 , as observed in the experimental dataset D: y opt ¼ arg max y flogðPðDjM; yÞÞg ð4Þ LL opt ¼ max y flogðDjM; yÞg ð5Þ We fitted the model separately on the data for each condition (stable, random, overlapping) in order to observe potential differences in the optimized parameters between conditions. In order to search for the model's LL opt and θ opt , we employed a gradient-descent method with Matlab's fminsearch function initialized at 1,000 initial parameter-sets (α,β) randomly sampled in [0;1]×]-1;1[.
Acknowledgments We would like to thank  Jop Niessen  for creating the dwell time maps and all students who have performed the object space task for this and other projects and who have convinced us of the reliability of the task:  Gu ¨lberk Bayraktar ,  Nikkie Cornelisse ,  Hussein Ghareh ,  Federico Giuliani ,  Sidney Hulzebos ,  Joanne Igoli ,  Stefanos Loizou ,  Luc Nijssen ,  Luca Reinink ,  Olaf Stoutjesdijk ,  Minou Verhaeg ,  Thijs Groenveld ,  Thom Joosten ,  Liz van den Brand ,  Violeta Caragea ,  Robin Lempens ,  Rian Kraan ,  Koen van den Berg ,  Wessel Fanselow ,  Daan van den Tillaart, Nikki van.
We took advantage of recent developments of large biodiversity databases and museum collections and we compiled a database of over 19 million records of flower visitor occurrences (Supplementary Table  1 ), spanning the period 1960-2016. This database includes 2023 European species from the 4 main insect orders of pollinators: Hymenoptera, Diptera, Lepidoptera and Coleoptera (Extended Data Fig. 1 ). Numerous species exhibit distinct modes in their phenology, either because they are multivoltine (i. e. multiple generations per year) or because the phenology differs between sexes or social casts. Since different modes from a species can shift in a different direction, we studied each mode separately, leading to 2248 phenology modes (see Methods). For each phenology mode, we estimated changes in mean flight date and flight period length over the years by modeling the mean and variance of collection dates (see Methods).
Assembling data on flower visitor occurrences in time and space European flower visitors mainly belong to four insect orders -Coleoptera, Diptera, Hymenoptera and Lepidoptera  10. We first looked for occurrence data (i. e. sighting at a given date and location) of species that belong to these insect orders and that are defined as floricolous in scientific or grey literature. We restricted our search to European species listed in Fauna Europaea  33. Data are from 15 distinctive sources, summarized in the Supplementary Table  1 , with a high proportion from the Global Biodiversity Information Facility (GBIF). After the removal of duplicates (same species, date and locality), the database initially included about 30 million of occurrences between 34° and 72° of latitude North and between -15° and 32° of longitude.
Several modes were initially detected for 494 species. For each of these species, we checked in scientific and grey literature whether a multimodal phenology was expected. In 208 cases, there was no strong biological support of existing multimodal phenology and we thus considered these species had one single mode. After this step, 288 remaining species showed a multimodal phenology (nbmax > 1). We applied the second step only for these species. The second step of the method was to attribute each record to a specific mode. To do so we used clustering Gaussian mixture-models implemented in the mclust R package  34  , considering a number of gaussians in 1 to nbmax. This clustering model allow us to initialize the attribution of each record to a given mode. Using the classification given by these clustering models, we run linear mixed-effects models, similar to the one described in equation (  1 ) but with the addition of a mode effect ( ): ( We kept the number of modes that minimize the Bayesian Information Criterion (BIC) of this linear mixed-effects model.
We then manually changed the mode of poorly predicted points. If the change improved the likelihood of this mixed-effects model we retained it and continued this process iteratively. We stopped the process when changing the mode of poorly predicted points did not further improve the likelihood of the model. The R script of the full method is available here: https://github. com/f-duchenne/Flower-visitors-phenology. Although the mode effect ( ) is independent from spatial variables and altitude in equation (  2 ), our method still allows to take into account spatial and altitudinal variation in the number of modes (Extended Data Fig. 7 ). We confronted the relevance of detected modes regarding what we know on the biology of species. We found that our method distribute records among modes in a highly consistent way.
Database after selection process Following the separation of distinct phenological modes for each species and the removal of larval records, we selected phenologies (or phenological modes) with at least 400 records during the period 1960-2016 and with at least 40 records from the period 1960-1980, to be able to study phenological shifts between early and more recent period. We removed species (n=30) with phenology peaking during winter by excluding species with a mean flight date before 60 or after 306 Julian days. Studying the phenology of such species raises methodological questions that we will not address here. We also removed records with imprecise localization (above 1km²) except for small countries (Luxembourg, Belgium, Switzerland, Netherlands, Denmark, Lichtenstein, Monaco, Andorra and Kosovo).
Thus, our dataset includes some records with imprecise localization (above 1km²) but they represent less than 0. 1% of the final dataset. This selection process lead to 19 845 792 occurrence records with 2248 phenologies for 2023 species (Supplementary Table  1 ). The repartition of records among insect orders and throughout the study period is presented in Extended Data Fig. 1. Supplementary Table  1  indicates the amount of data coming from the various data sources.
Estimating species phenological shifts To estimate changes in both the mean flight date (MFD) and the flight period length (FPL), we modeled jointly the mean and the variance of collection dates using the dispmod R package  35  which performs two nested linear models, one for the mean and one for the variance. Due to computational limits it was not possible to use one model including the whole dataset, modelling both, MFD shifts and FPL changes, and modelling spatial effects properly for each species. Thus, we studied each species and phenology mode separately. For each phenology modes, the model for the mean collection date was: (3) is the Julian day of the observation k, is the grand mean (intercept), is the time effect on the mean collection date as well as on its variation across latitude ( ) and longitude (.
Moreover, because these polytomies were not too old relative to the entire phylogeny, it should not affect strongly our results. Because they are not present in our phylogeny three families of Diptera (Heleomyzidae, Limoniidae and Pediciidae) and two Lepidoptera species (Sphrageidus similis, Lymantriidae, and Agria desoptilete, Lycaenidae) were excluded from phylogenetic analysis. We estimated phylogenetic signal in phenological shifts using Pagel's λ 40 implemented in the phylosignal R package  41  , because it is much more robust to polytomies than Blomberg's K  42.
Finally, in order to have more robust values, we repeated this overlap calculation after shifting segmentation of the latitude and the longitude by 1. 25°, 2. 5° and by 3. 75°. Then we averaged values obtained by 1. 25°×1. 25° grid cells for both measures, overlap within and among orders. In order to study the seasonal dynamic of overlap changes, we calculated a proxy of the phenological overlaps day by day in 1980 and in 2016 for each grid cell (Extended Data Fig. 5 ). We do not use exactly the same calculation of overlap as previously for computational reasons. To simplify the calculation method, we aggregated predicted phenologies at the order level to get a density distribution by order, henceforth called order phenologies, as presented in Fig. 3a. Then we calculated the pairwise overlap among order phenologies day by day for both years, 1980 and 2016, and for each grid cell.
We also evaluated the day-by-day density value for each order phenologies for both years, 1980 and 2016, and for each grid cell. This density value is a proxy of the phenological overlap within order, because we assume that every species has the same constant abundance. Then we calculate the daily changes of these both indexes between 1980 and 2016 (Extended Data Fig. 5 ). We did so for one grid pattern only (i. e. without sliding windows).
D. , C. F. and E. T. conceived the project. F. D. assembled dataset and performed the statistical analysis. F. D. , C. F. , E. T. , D. M. and M. E. interpreted the results and M. D. , M. Pe. , J. S. P. , M. Po. & P. V. provided data, biological expertise on studied species. F. D. wrote the paper with contributions from all authors.
Acknowledgements We especially thank  Nicolas Loeuille  and  Tom Van Dooren  for their comments on this work,  Emmanuelle Porcher ,  Adrien Perrard ,  Elsa Teulière ,  Emile Faure ,  Elise Kerdoncuff ,  Benoit Perez  and  Théophile Olivier  for fruitful discussions. We also thank all naturalists who provided data to complete our database, especially  Lionel Casset ,  Fred Chevaillot ,  Laurent Guilbaud  ( INRA Avignon ),  Jean-Laurent Hentz  ( Nature du Gard ) and  Guillaume Lemoine , as well as  Christophe Daugeron ,  Rodolphe Rougerie ,  Claire Villemant ,  Emmanuel Delfosse ,  Jérôme Barbut  and  Olivier Montreuil  who helped us to access insect collections of the  French National Natural History Museum. This project was funded by the  Ministère de la Transition Ecologique et Solidaire as part of the project  " What is the sensitivity of pollinators to global warming in France " (convention n° SJ 3-17) lead by  Christophe Daugeron  and  Colin Fontaine , and by the  Institut de la Transition Ecologique, Sorbonne Université , as part of the project Yapludsaison.
F. D. , C. F. and E. T. conceived the project. F. D. assembled dataset and performed the statistical analysis. F. D. , C. F. , E. T. , D. M. and M. E. interpreted the results and M. D. , M. Pe. , J. S. P. , M. Po. & P. V. provided data, biological expertise on studied species. F. D. wrote the paper with contributions from all authors.
For global annotation of the libraries (Supplementary files 5 and 7), we used the release 5. 49 of fasta reference files available in Flybase, including transposon sequences (dmel-all-transposon_r5. 49. fasta) and the release 20 of miRNA sequences from miRBase (http://www. mirbase. org/). Sequence length distributions, small RNA mapping and small RNA overlap signatures were generated from bowtie alignments using Python and R (http://www. r-project. org/) scripts, which were wrapped and run in Galaxy instance publicly from ARTbio platform available at http://mississippi. fr. Tools and workflows used in this study may be downloaded from this Galaxy instance. For library comparisons, read counts were normalized to one million miRNA (Supplementary files 5 and 7). A second normalization, performed using the total number of small RNAs matching the D.
melanogaster genome (release 5. 49), gave similar results (Supplementary files 5 and 7). For small RNA mapping (Figures  2  and 4 , Figure  2 -figure supplements 1 and 2), we took into account only 23-29 nt RNA reads that uniquely aligned to reference sequences. Logos were calculated using Weblogo  (Crooks et al. , 2004)   2009 ) and detailed in  Antoniewski (2014). Thus, for each sequencing dataset, we collected all of the 23-29 nt RNA reads matching P(lacW) or the 42AB locus whose 5' ends overlapped with another 23-29 nt RNA read on the opposite strand. Then, for each possible overlap of 1 to 29 nt, the number of read pairs was counted. To plot the overlap signatures, a z-score was calculated by computing, for each overlap of 1 to i nucleotides, the number O(i) of read pairs and converting the value using the formula z(i) = (O(i)-mean(O))/standard deviation (O).
The percentage of reads containing a 'A' at the tenth position was calculated within the paired 23-29 nt RNA matching the reference sequence as described in de Vanssay et al. (  2012 ) (BX2 transgene in Figure  2C  and 42AB in Figure  2 -figure supplement 1B). GRH49 (BX2*) was previously analyzed in  Hermant et al. (2015). Small RNA sequences and project have been deposited at the GEO under accession number GSE116122.