dataesr/scientific-paragraphs-categorization huggingface dataset:
The sampling spans more than one year with 51 samples collected from March 2011 to May 2012 from which we co-extracted the macrobiomolecules and analysed the derived MG, MT and MP datasets  35  alongside the physicochemical factors measured at the site. We reconstructed the MG structure of the community, alongside its taxonomy, genetic potential, transcript and protein levels. We employed SVD to extract relevant temporal patterns, which were then clustered into 17 fundamental signals. Those were integrated with collected environmental parameters to build an ARIMA model, augmented with seasonal components, which could explain the observed signals. Multiple models (ARIMA, prophet  36  and NNETAR neural networks model  37  ) were trained to forecast the following five years' signals. Validation was conducted using future time points, i.
A KEGG Orthology group (KO term) was assigned to 55% of the total retrieved ORFs, whilst taxonomic affiliations were assigned to 38. 5%. The number of ORF copies as well as their detected gene expression and protein abundances were determined over the extended dataset (see methods). We found on average 2. 2×10 6 (s. d. 4. 8×10 5 ) ORFs, 9. 1×10 5 (s. d. 1. 7×10 5 ) transcripts and 2. 4×10 5 (s. d. 2. 5×10 4 ) protein groups per sample. However, the vast majority of the genes were not found to be expressed over the entire dataset or were only detected in a few samples, with a maximum of 16. 8×10  6  ORFs detected in one sample. This suggests that a significant fraction of the gene pool in LAO is not specifically required for community function but rather their cumulative functional effort may be compartmentalised, fitting the previous results from Roume et al.
The temporal signals underlying the microbial community Considering that the information necessary to forecast the community dynamics and linked gene expression may be most represented in any biological (e. g. taxonomical or functional representation) or environmental data layer, we decided to include multiple layers in our analysis. Regarding the microbial community we explored multiple taxonomic and functional levels at once and summarised their temporal characteristics. Thus, the three quantification matrices (MG, MT and MP) were used to compute "summary" matrices according to the ORFs' descriptors. Hence we computed one matrix per omic layer for the six formed taxonomic descriptors (Phylum, Class, Order, Family, Genus and Species) and two functional ones (KO terms and pathways). The resulting 27 matrices (3 original and 24 summary) were used to compute the system's eigengenes (EGs) i.
e. the independent patterns underlying the data and their potential time dependency  27. In previous works, the first EG in a time-series has been demonstrated to represent the "steady state" gene expression, encapsulating the largest explained variance (EV), and was removed to perform the timeseries analysis  27. Indeed, the first EG showed the largest variance explained (around 15% in all the datasets), therefore we excluded it from the subsequent analysis. We screened the subsequent EGs for time-dependency (see Methods) selecting a set of 210 EGs and assessed how much of the data variation they explained beyond the first EG (Supplementary Figure  4 ). Considering that the EGs of each matrix are linearly independent (i. e. they do not have redundant information) for each matrix, we expected some level of redundancy by using different types and levels of summarising the information.
Forecasting of future time points From the analysis of the signals identified in the training datasets it is already possible to identify five signal groups: i) alternative basal states, e. g. two alternative stable states of abundance/expression, (S5, S14); ii) perturbation, i. e. standing wave with varying amplitude and frequency, (S4, S8, S15); iii) cyclical, i. e. standing wave with constant amplitude and frequency, (S6, S11, S12); iv) "crashes", i. e. quick shift in the state and reversion to basal state, (S1, S10, S16); v) mixed, i. e. the other factors(Supplementary Figure  6, panels c-f ). Alternative stable states, perturbations and crashes (groups i, ii and iv) are hard to model without observing multiple times the shift and the perturbation events, respectively. Additionally, these scenarios may include permanent shifts into a new community equilibrium or transitory signals in the community that will be eventually resolved (e.
g. a viral infection). To forecast such events and based on results, systematic information on microbial interactions would be required which is beyond the scope of this study. The used modelling also enables us to forecast cyclical events (group iii). The 17 signals were used to train three models (with various parameters) from the package fable  37  and the best performing model on the training set was selected for each of them (see methods). In detail, ARIMA, prophet and neural network models (with up to four Fourier terms for ARIMA and prophet) were trained for the S1-17 using the environmental variables as external regressors. The 51 weeks spanning 2011-2012 data were used as a training set and the model with the smallest Root Mean Square Error (RMSE) was selected for forecasting. A total of 21 new samples were collected in the month of June of the following 5 years to validate the model for the MG and MT data.
Data and code availability The generated MG and MT reads (FASTQ) files, as well as the previously produced data, are available as NCBI BioProject PRJNA230567. The MP data from the PRIDE repository with accession number PXD013655  38. The meta-omics pipeline IMP v3. 0 46 is maintained and developed at the GitLab page: https://git-r3lab. uni. lu/IMP/imp3. The code used in the analysis is available at https://github. com/fdelogu/microforecast, whilst the data required to start the analysis is available on Zeonodo with the doi 10. 5281/zenodo. 7225349.
Collection of environmental variables The environmental variables were collected on site by the researcher(s) whilst they were performing the sampling, which include: dry matter, phosphate, nitrate, ammonium, oxygen, conductivity, pH, temperature and oxygen (Supplementary Table  4 ) following the previously established protocol  38. The other variables were retrieved from the automated data collection routine of the Schifflange BWWTP, which measures online these values and aggregate them as 2h average starting at 1:00 am. Those recordings include the same variables for different parts of the plant (inflow, both vats, outflow) with the addition of other measurements such as the in/out flow volume. For simplicity, we used exclusively the variable pertaining to the inflow, both vats and outflow in this study (Supplementary Table  5 ).
Co-assembly of metagenomics and metatranscriptomics reads All the samples from the training and the test datasets followed the same bioinformatic pipeline. Sample-wise preprocessing of the MG and MT data was performed using IMP v3. 0  46  (https://git-r3lab. uni. lu/IMP/imp3) with custom parameters, i. e. i) Illumina Truseq2 adapters were trimmed, ii) the step involving the filtering of reads of human origin step was omitted for the preprocessing. The reads were corrected using BayesHammer 48 per sample, per omic. The resulting MG and MT reads were assembled with metaSPAdes v3. 13. 1  49  and rnaSPAdes v3. 13. 1  50  respectively. The MG and MT reads of each sample were re-assembled together using the contigs and "highly filtered" transcripts from the first assemblies as trusted contigs.
Binning and clustering The chromosomal prokaryotic subsets of each sample were binned using IMP v3. 0 46 with MaxBin 56 , MetaBAT 57 , binny  58  plus a refinement step with dastool. The resulting bins dereplicated along the entire time series with dRep  59  to create representative metagenome assembled genomes (rMAGs). Similarly, the eukaryotic subsets were binned with MetaBat 60 and dereplicated using dRep  59  without genome quality assessment resulting in rMAGs. All the plasmidial, viral and the unbinned contigs from the eukaryotic and chromosomal prokaryotic subsets were clustered using CD-HIT  61  on each of those subsets. We refer to the subset of the clustered unbinned contigs as representative contigs (rContigs). The collection of the rMAGs and the rContigs constitute the representative database (rDB) of the system.
MG and MT quantification and filtering The filtered MG and MT reads were aligned to the ORF reference set using bwa  64  and sorted using samtools  65. The resulting sorted bam files were processed using bam2hits  66  , and the output split with a maximum number of 100'000 ORFs per subset, whilst respecting the bam2hits read groups. Each subset was quantified with mmseq 66 and mmcollapse  67  , then the quantifications per sample were the-normalized form FPKM, merged and re-normalized to FPKM. Values of gene abundance and expression inferior to 10 -7 were considered equal to 0 and ORFs and transcripts that were not present in at least 20% of the training set were discarded from further analysis.
MP quantification and filtering Raw MP data were retrieved from the PRIDE repository with accession number PXD013655  38  , where the samples were processed as described in Muller et al. 47  , and we reanalyse them. Supplementary Figure  10The  complete set of predicted ORFs was subsetted to obtain smaller sample-specific databases. The MG alignment files generated in the previous step were processed with featurecounts 68 and all the ORFs with a count greater than 0 for the given sample were included in the appropriate sample. Each sample-specific database was concatenated with a cRAP database of contaminants (https://thegpm. org/cRAP; downloaded in July 2019) and the human UniProtKB Reference Proteome (UniProt Consortium, 2021), and decoys were generated by adding the reversed sequences of all protein entries to the databases for the estimation of false discovery rates.
The search was performed using SearchGUI v. 3. 3. 20  69  with the X!Tandem 70 , MS-GF+  71  and Comet 72 as search engines and the following parameters: trypsin was used as the digestion enzyme and a maximum of two missed cleavages was allowed. The tolerance levels for matching to the database were 10 ppm for MS1 and 15 ppm for MS2. Carbamidomethylation of cysteine residues and oxidation of methionines were set as fixed and variable modifications, respectively. Peptides with length between 7 and 60 amino acids, and with a charge state composed between +2 and +4 were considered for identification. The results from SearchGUI were merged using PeptideShaker-1. 16. 45  73  and all identifications were filtered to achieve a protein false discovery rate (FDR) of 1%. The sample-specific peptide-spectrum matches (PSM) obtained for each analysis were then used to calculate dataset-wide protein groups using the Occam subgroup method from the Pout2Prot algorithm  74.
Eigengenes and their analysis The EGs for the training set (samples from 2011-03-21 to 2012-05-03) were computed as singular right eigenvectors obtained with the function svd. The data were normalised according to the basal expression  27  computing the quantification matrices as: For each signal we tested ARIMA with up to four Fourier components, Prophet with up to four seasonal components and a neural network with 10, 20 and 30 nodes in the hidden layer. The models were scored according to their RMSE and the top three were combined (weighted by 1-RMSE) and used as a fourteenth model, for which the RMSE was calculated too. The best of the fourteen models was selected for each signal and used to forecast the test test with the function forecast form the fable package and supplying environmental parameter readings.
variables and the 17 signals as explanatory variables. We then rebuild the test matrices multiplying the forecasted signals over the test set time with the newly calibrated betas and adding the intercept (basal level). The reconstructed samples were compared with the original ones on an individual basis (Supplementary Figure10) and on an average one (Figure4). to P. W. as well as CORE Junior (C15/SR/10404839) grant to EELM. The project received financial support from the Integrated Biobank of Luxembourg with funds from the Luxembourg Ministry of Higher Education and Research. This work was also supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 863664). The work of P. M. was funded by the 'Plan Technologies de la Santé du Gouvernement du Grand-Duché de Luxembourg' through the Luxembourg Centre for Systems Biomedicine (LCSB), University of Luxembourg.
Acknowledgments We thank the  Luxembourg National Research Fund (FNR)  for supporting this work through various funding instruments. Specifically, a  PRIDE doctoral training unit grant  ( PRIDE/15/10907093 ),  CORE grants  ( CORE/17/SM/11689322 ), a  European Union   ERASysAPP grant  ( INTER/SYSAPP/14/05 ), and an  ATTRACT grant  ( A09/03 ) all awarded.
e. , iceberg calving frequency seen in Fig. 6b  and the relationship between annual cumulative AR intensity and melt/runoff seen in Fig. 6c ) is the summation of maximum IVT over the course of a distinct AR landfall event thus incorporating the duration and IVT strength of the storm. Daily cumulative intensity is like the cumulative intensity except we make the summation of maximum IVT for daily periods when an AR landfall is detected. This is used for testing the probability of an extreme state (i. e. , temperature extreme, melt extreme, runoff extreme, 2-day sea-ice removal, high swell seen in Fig. 6a ) Regional climate model and reanalysis datasets. MAR is a regional climate model specifically designed for simulating polar climate  24, 66. MAR atmospheric dynamics are based on the hydrostatic approximation of the primitive equations  67.
Nevertheless, we analyzed the well documented collapse of the Larsen A ice shelf  15  using AVHRR satellite images available on the National Snow and Ice Date Center (NSIDC) website  81  surrounding the date of the Larsen A collapse. Winter calving events were not retrieved with the same accuracy because even though thermal bands could be used during polar night, this imagery technique is more sensitive to cloud cover (even thin cloud are opaque). This lowers the resolution even given favorable cloud cover, thus making the acquisition of data at a daily time scale not straightforward. Even if we present the Larsen C calving event of July 10-12, 2017, we did not consider this event in the final statistics. This prompted us to consider the iceberg retrieval database from Stuart and Long, 2011  82.
Here, we extracted the initial location of icebergs around the Larsen ice shelves, but difficulty was encountered in tracing with confidence the origin of these icebergs to a given ice-shelf because the first detection of an iceberg can occur days after calving and its position can be far from its parent ice-shelf. Nevertheless, the information of proximity with an ice shelf provides some hints on calving in the Larsen area. In the latter database, winter occurrences are also largely less frequent. Similar examination was more limited over the Wilkins ice shelf because cloud cover is generally denser, so collapses and calving events are generally observed a few days after their occurrence. Moreover, contrary to the eastern side, several large events that occurred in the darker autumn and winter periods were not visible in MODIS datasets.
This made the attribution of ARs to calving occurrences more difficult on the western side of the AP. Surface melt surface anomalies over the Larsen ice shelves were also retrieved between October 1 and March 31 every year over 1979-2020  83. A series of microwave radiometer observations (namely Scanning Multichannel Microwave Radiometer (SMMR), the Special Sensor Microwave/Imager (SSM/I) and Special Sensor Microwave Imager/Sounder (SSMIS)) were processed to test our conclusions on AR impact on snow melt (see Supplementary Fig. 7 ). First, a climatology of daily melt surface extent based on mean values for each Julian day mean was computed over the full 41 year-long period. The daily melt anomalies were computed by comparison with this climatology and compared with AR occurrences. Data obtained before 1993 must be considered with caution, in particular with summer 1987-88  31  Finally, to complement our analysis on the relationship between melt pond observations and AR activity obtained using the MODIS image interpretation  34  , we also studied the Landsat8 image database  84.
Acknowledgements This study is part of the Ph. D project of J. D. W. conducted at the  Université Grenoble Alpes. Particular thanks to  Simon Gascoin  for providing help in the MODIS images retrieval and to  Alison Delhasse  for her monitoring of the MAR simulations. We acknowledge support from the  Agence Nationale de la Recherche  projects  ANR-20-CE01-0013  ( ARCA ),  ANR-14-CE01-0001  ( ASUMA ),  ANR-16-CE01-0011  ( EAIIST ),  ANR-15-CE01-0005-01  ( TROIS AS ) and  ANR-15-CE01-0015  ( AC-AHC2 ). I. V. G. thanks  FCT/MCTES  for the financial support to CESAM ( UIDP/50017/2020  and  UIDB/ 50017/2020 ) through national funds and support to project  ATLACE  ( CIRCNA/CAC/ 0273/2019 ) funded by  FEDER , through  COMPETE2020-Programa Operacional Competitividade e Internacionalização (POCI) , and by national funds (OE), through  FCT/ MCTES.
C. A. acknowledges support from the  Fondation Albert 2 de Monaco  under project  Antarctic-Snow  ( 2018-2020 ). Computational resources were provided by the  Consortium des Équipements de Calcul Intensif (CÉCI) , funded by the  F. R. S. -FNRS  under grant no. 2. 5020. 11. J. V. T. acknowledges the support of the  German Ministry for Education and Research (BMBF)  grant number  03F0855.
Data availability The Modèle Atmosphérique Régionale (MAR) data is available at https://doi. org/ 10. 5281/zenodo. 6347190. We acknowledge the use of imagery from the NASA Worldview application (https://worldview. earthdata. nasa. gov), part of the NASA Earth Observing System Data and Information System (EOSDIS). NOAA AVHRR imagery was acquired from the NSIDC  81. The output from the AR detection algorithm can be downloaded from the ARTMIP database at https://www. earthsystemgrid. org/dataset/ ucar. cgd. ccsm4. artmip. tier1. html.
| Studied population and experimental design The focus population (F20089 or CO3 according to  Jullien et al. , 2019)   were successfully multiplied. Out of these, 55 families for each of the two sampling years were randomly chosen in 2012. Seeds from the 110 families were scarified to ease germination and were transferred in Petri dishes with water at room temperature for six hours. We then used two different vernalization treatments (at 5°C during 7 or 14 days) to compare the vernalization requirement between the two years. Five replicates from each vernalization treatment were transferred back to the greenhouse, according to a randomized block design (five blocks and two treatments, adding up to a total of ten replicates per family, 1100 plants in total). Data loggers were placed on each table to monitor temperature and humidity.
| Temporal changes in flowering time Individual flowering times were converted to thermal times following  Bonhomme (2000). The thermal time was calculated as the sum of the mean daily effective temperatures of each day between sowing and the emergence of the first flower, where the mean daily effective temperature is the day's mean temperature minus the base temperature (T b ). We used T b = 5°C, as reported by  Moreau et al. (2007)  for the M. truncatula reference line A17. Plants noted as sick or failing to produce leaves were removed from the datasets (22 individuals removed). Collected measures were tested for normality using quantile-quantile (Q-Q) plots  (Nobre & Singer, 2007). All analyses were conducted using R version 2. 15. 2. We used linear mixed models (lme4 package) to test for a significant change in flowering time between the sampling years.
| Genetic analyses During the multiplication generation in the greenhouse (2011), 200 mg of leaves was sampled from each plant for DNA extraction, using DNeasy Plant Mini Kit (Qiagen). Twenty microsatellite loci were used for genotyping (see the details of amplification reactions and analyses of amplified products in  Jullien et al. , 2019; Siol et al. , 2007). Briefly, samples were prepared by adding 3 μl of diluted PCR products to 16. 5 μl of ultrapure water and 0. 5 μl of the size marker AMM524. Amplified products were analyzed on an ABI prism 3130 Genetic Analyzer, and genotype reading was performed using GeneMapper Software version 5.
| Single-locus analyses assuming independence among loci As a preliminary step, the data were filtered to reduce the percentage of missing data (loci or individuals with >10% missing data were removed), and to discard monomorphic loci. After filtering, the dataset comprised 145 individuals (representing 145 families) and 16 loci (64 from the year 1987 and 81 from the year 2009). We measured the genetic diversity of the population each year using the allelic richness N a -rar  (Hurlbert, 1971)  and the expected heterozygosity H e. In this predominantly selfing population, we expect a strong deviation from Hardy-Weinberg heterozygosity expectations. Thus, for each sampling year, we estimated the inbreeding fixation coefficient F IS and its confidence interval using 5000 bootstraps over loci.
Between year differences for N a -rar , H e and F IS across loci were tested using Wilcoxon signed-rank tests. Analyses were performed in R using the packages adegenet  (Jombart, 2008)  and hierfstat  (Goudet, 2005)  and the program ADZE for rarefaction analyses  (Szpiech et al. , 2008). The percentage of pairs of loci showing significant linkage disequilibrium (LD) was calculated using Genepop  (Rousset, 2008)  with a threshold of 0. 05. Finally, we measured the temporal variance in allele frequencies using the F ST estimator by  Weir and Cockerham (1984). To estimate the effective population size (N e , measured in number of diploid individuals) from the temporal variance of allele frequencies, we used F ST estimates to account for the correlation of alleles identity within individuals due to selfing  (Navascués et al.
| Analyses based on multilocus genotypes We used the program RMES to estimate selfing rates from the distribution of multilocus heterozygosity  (David et al. , 2007). We tested for a difference in selfing rates between years using a likelihood ratio test between models where the selfing rate was constrained to be constant or not. For each sample  (1987 and 2009) , we examined the genetic structure by sorting out the number of multilocus genotypes (thereafter called MLG) and measuring their frequency and redundancy through time using GENETHAPLO (available on GitHub at https://github. com/lauga y/Genet Haplo and described in Appendix S1: Section S1). GENETHAPLO takes into account the uncertainty of the assignment of a genotype to a MLG group due to missing data: In case of ambiguity, an individual is randomly assigned to one of the candidate MLG group with a probability proportional to the MLG group size.
The approach also considers a genotyping error rate: If two individuals differ by less than the error rate, they are considered to belong to the same MLG. After an initial run with an error rate of zero, we checked the distribution of the distances between MLGs. We found an excess of small distances, which could indicate errors in genotype assignation  (Arnaud-Haond & Belkhir, 2007). We corrected this by re-running the program with an error rate of 1/16 (= one mis-read locus). GENETHAPLO also searches for residual heterozygosity (defined as the proportion of heterozygous loci in the multilocus genotype) and evidence for recombination (S1). To identify putative recombination events between MLGs, it uses the genetic distances: a MLG is a recombinant candidate if the sum of its allele differences with two other MLGs ("parental MLGs") equals the number of allele differences between these two parental MLGs.