dataesr/scientific-paragraphs-categorization huggingface dataset:
Target Words Selection The evaluation dataset that we use in this paper includes 100 preselected words, so we are able to apply clustering for each of them. However in any practical application of semantic change detection, performing clustering for every word in the corpus would not be feasible in terms of computing time. Thus, we investigate a preliminary step to select a set of words that may have undergone semantic change. We use several scalable metrics, that we compare in Section 5. 1. A first set of metrics relies on the computation of a variation measure: for each word, similarly to  [21] , we compute the average token embedding on the full corpus. Then, we compute the cosine distance between each token embedding and this average embedding. Finally, we take the mean of all these cosine distances as the variation coefficient of the word.
Preliminary Step: Target Word Selection Table  1  shows the Pearson and Spearman correlations between the various preliminary drift measures defined in section 4. 2, and the human-annotated drifts. All metrics are computed using the full dataset from the 1960s to the 1990s, except averaging which is only computed between the first decade (1960s) and the last one (1990s) without using the textual data in between. The measure with the highest correlation with the human annotations is averaging, which measures semantic drift between the first and the last time step. This intuitively makes sense as it resembles the way the evaluation dataset was annotated by only considering the first and the last decade. Variation by decade also shows good results; it is a measure of the evolution of the level of variation of a word usage through time.
Impact of Fine-tuning Figure  2  shows the comparison of fine-tuning influence for two best clustering methods (affinity propagation, and k-means with k=5). Interestingly, a light fine-tuning (just for one epoch) decreases the performance of both methods (in terms of Spearman correlation) in comparison to no fine-tuning at all (zero epochs). After that, the length of fine-tuning until up to 5 epochs is linearly correlated with the performance increase. Fine-tuning the model for five epochs appears optimal. After that, the performance for both methods starts decreasing, most likely because of over-fitting due to the reduced size of the fine-tuning dataset compared to the training data. The impact of fine-tuning on the k-means clustering is stronger than on the affinity propagation. The difference between model's performance on 5 epochs is negligible.
FUTURE WORK One major issue in quantifying semantic change over time is the lack of gold standard data. In this work we used the 100 humanannotated words from  [13] ; but to extend this study to a wider range of words for other corpora and other time slices we can construct datasets that simulates synthetic change as done in  [30]. This involves constructing "pseudowords" that are concatenations of two real words, and then injecting it in place of the real words with controlled probability distributions over a series of time slices. We also plan to investigate how the clusters found by the methods in this work can be used to interpret the different usages of a word in a specific time slice. The initial experiments on this subject have already been conducted with the two-stage clustering, which removes the smallest clusters, containing one or two instances.
ACKNOWLEDGMENTS We are grateful to  Andrey Kutuzov  for valuable discussions during this paper preparation. We also thank  Pr. Alexandre Allauzen  from  ESPCI -Université Paris Dauphine  and  Pr. Asanobu Kitamoto  from  National Institute of Informatics (Tokyo)  for their advises. This work has been partly supported by the  European Union 's  Horizon 2020 research and innovation programme  under grants  770299  (NewsEye) and 825153 (EMBEDDIA).
More details of the sampling procedure and the software used for bioinformatic processing of the metagenome data are summarized in Table  S1. The code used in this study for data processing is available via https://www. proto cols. io/view/ allel e-varia bilit y-bkhqkt5w. In short, all raw sequence data were quality trimmed and co-assemblies were created for each time-series with the exception of the Baltic Sea metagenomes, where the publicly available BARM assembly  (Alneberg et al. , 2018)  was used. CDS (coding DNA sequence) regions were identified by gene calling and subsequently clustered by 100% amino acid sequence similarity using the cd-hit software  (Li & Godzik, 2006). Functional annotation of the CDS regions was performed using the KEGG database  (Kanehisa et al. , 2007). We defined all gene variants coding for the same gene orthologue but differing in their amino acid sequence as alleles of this gene.
However, we believe that such inaccuracy will be similar among the considered gene categories and should therefore not bias our downstream analyses in a certain direction. The quality-trimmed data were mapped on the respective co-assemblies using bowtie2 set to very-sensitive-local  (Langmead & Salzberg, 2012)  and summarized using the featurecounts software  (Liao et al. , 2014). For downstream analyses, we focused on KEGG orthologues that were present in fewer than two-thirds of all prokaryotic genomes in the KEGG database (available at: https://www. kegg. jp/ kegg/downl oad/; release 78. 1; downloaded: May 2016) containing data from 3,865 prokaryotic genomes), which we refer to as auxiliary genes. Our decision to exclude gene orthologues that are present in most of the prokaryotic genomes was motivated by the fact that a high proportion of these genes encode proteins involved in essential functions, such as replication, transcription or translation (Table  S3 ).
| Gene copy number and seasonal metacommunity size For each gene orthologue identified via the KEGG ontology, the average per-cell gene copy number among carrier genomes was determined based on the average gene copy number in prokaryote organisms of all entries coding for the respective gene orthologue (KEGG; release 78. 1). In the remainder of this paper we refer to the average copy number of a gene across all carrier genomes in the KEGG database simply as "gene copy number. " While the presence of different communities at each of the examined sites will result in site-specific gene copy numbers, this information cannot be determined without access to the full sequence information for each genome present at this site. However, metagenome data do not allow us to reconstruct full genome data from all community members and we therefore use the full genome data of all prokaryotic genomes stored in the KEGG database to obtain a value for the gene copy number of individual genes that was then used to approximate the true values in the communities from all sample sites.
We thereby assume that the ranking from single-copy genes to high-copy-number genes will be similar across different ecosystems. The gene-specific seasonal metacommunity size was assessed by dividing the total number of alleles for the respective gene orthologue detected throughout each metagenome time-series divided by the gene copy number of this gene orthologue. To match data from the downstream allele variability analyses, we estimated the metacommunity size only after subsampling allele count-data for the respective genes to 500 counts per individual metagenome in each time-series. The resulting parameter estimates the seasonal metacommunity richness of carrier organisms for each gene orthologue and represents a parameter for temporal gamma diversity. The term metacommunity has been defined as local communities that are connected by dispersal of multiple potentially interacting species  (Wilson, 1992) , while for this study we use the term seasonal metacommunity to describe several communities from the same location but across different seasons.
| Subcellular location The subcellular location of all entries in the KEGG database (release 78. 1) was determined based on the presence of signal peptides in the amino acid sequences via the psortb (version 3. 0) software  (Yu et al. , 2010). We considered all sequences affiliating according to the KEGG taxonomy with the phyla Firmicutes or Actinobacteria as originating from gram-positive bacteria  (Gontang et al. , 2007)  and the remaining bacterial sequences as originating from gram-negative bacteria. We applied the settings for gram-negative bacteria, gram-positive bacteria or archaea for the KEGG entries affiliating with either of these groups, respectively. To test Hypothesis 3, gene orthologues were defined as cytoplasmatic if >80% of all entries were identified as coding for cytoplasmatic proteins.
An R script with the code to estimate gene allele variability and perform the described statistical tests is available via https://github. com/sarab eier/allele. varia bility. Gene orthologues that exceeded 500 counts in the individual samples for all time-series, and that were therefore included into the analyses, were ranked in each time-series by their allele variability. The variance of the gene variability rank was estimated in order to identify genes with similar or very distinct variability across all five time-series.
| Overall environmental heterogeneity An index for environmental heterogeneity for each of the time-series was estimated by summing the coefficient of variance for temperature, salinity and chlorophyll a concentration (Table  S2 ) as published elsewhere  (Maurice et al. , 2013). Environmental data for MEN were retrieved from the North Temperate Lakes Long Term Ecological Research site, which also documents detailed experimental methods (https://lter. limno logy. wisc. edu). Salinity data for Lake MEN were not available, but we assumed the coefficient of salinity variance to be zero in this freshwater system without connection to marine waters. Environmental data for BATS and HOT were available via the homepage of the Bermuda Atlantic Time-Series Study (http:// bats. bios. edu/) and via the Hawaii Ocean Time-Series data provided html), respectively.
| RE SULTS Mean values for auxiliary gene variability ranged from 0. 40 in the HOT time-series to 0. 52 in the SOLA time-series. Mean allele variability of phylogenetic marker genes used as a proxy for the variability in taxonomic community composition exhibited similar values, ranging from 0. 40 to 0. 53 (Table  1 ). The ranking of mean allele variability in the five time-series was nearly identical for auxiliary genes and phylogenetic marker genes (Table  1 ). An index of environmental heterogeneity based on the temporal variability of temperature, salinity and chlorophyll a concentration in the five time-series (Table  1 ) indicated a higher positive correlation to the mean allele variability of auxiliary compared to core genes (Pearson's r =. 80 and r =. 63, respectively). However, due to the low number of data points (n = 5), none of these correlations was significant, apart from high correlation coefficients (p =.
10 and p =. 26, respectively). In all five time-series, we identified a gene orthologue encoding an outer membrane iron receptor protein (KEGG ID: K02014) as the temporally most variable auxiliary gene (Figure  1 ; Table  S5 ). The allele variability of other genes was less consistent across the five data sets. For instance, the alleles of restrictios enzymes encoded by the genes hsdR and hasdM (K01153 and K03427) were highly variable in SOLA, LMO and MEN, but belonged to genes with low variability in the offshore sites BATS and HOT (Table  S5 ). Gene variability assessed via the multivariate dispersion metric  (Anderson et al. , 2006)  based on the Bray-Curtis dissimilarity index will, analogously to a simple pair-wise comparison, be zero if the composition of subunits (in our case the alleles of a gene) in all compared samples is identical.
Consequently, new alleles and their carrier species can emerge in samples that were not present at the start of the time-series. Hence in our study, seasonal metacommunity size was defined as the richness of carrier species of the investigated genes measured across all sample dates within each time-series. Possibly, the increasing importance of stochastic effects with increasing seasonal metacommunity size could have contributed to the positive relationship between seasonal metacommunity size and temporal variability observed across our five metagenome time-series. However, another mechanism may in our case further strengthen this relationship, because our analyses are based on data from natural environments where the functional attributes should reflect environmental heterogeneity. A larger seasonal metacommunity size for a specific gene could essentially be caused by higher diversification rates for this gene as a consequence of naturally occurring high temporal heterogeneity of parameters that impact functioning of this particular gene  (McArthur et al.
For instance, parameters such as vitamin availability that often arise from biotic activity and indicate biological interactions may be essential for community assembly  (Herren & McMahon, 2017)  and accordingly also for temporal allele variability. The indicated correlation between the heterogeneity index and average allele variability therefore suggests that in our study temperature, salinity and chlorophyll a concentrations may have influenced the variability of other parameters not considered that have an impact on allele variability. Interestingly, within all time-series, a gene encoding multicopy outer membrane iron receptor proteins (K02014, including FhuA) was detected as the gene with highest temporal allele variability (Figure  1 ; Table  S5 ). Iron is a cofactor in many central reactions of organisms, such as respiration and photosynthesis, and is therefore essential for the viability of all cells  (Holm et al.
, 1996). Proteins summarized as the KEGG orthologue K02014 interact with siderophore-bound iron, and the variability of siderophore types excreted by different siderophore producers may partly explain the high allele variability of K02014. However, besides their function as iron transporters, these proteins also serve as targets for phage infection. This was highlighted by the "Ferrojan Horse Hypothesis," according to which phages incorporate iron atoms into their tail to be recognized and enter host cells via K02014 transporters  (Bonnain et al. , 2016). A recent metagenome-enabled study demonstrated that iron-binding motives were present in 87% of unigenes encoding marine viral tail proteins. Furthermore, the corresponding viral contigs were distributed ubiquitously and with high abundance in marine sites worldwide sampled during the Tara Ocean Expedition  (Caputi et al.
, 2019). These findings indicate that "ferrojan horse" phages may play an important role in the ecology of marine systems. It has been proposed previously that the interaction of prokaryotes with phages is a major mechanism for maintaining prokaryote diversity, particularly due to the variability of genes that interact with phage infections, such as K02014  (Cordero & Polz, 2014; Rodriguez-Valera et al. , 2009). The large allele variability of K02014 across all five time-series is thus a further indication for a possibly high relevance of predation pressure from "ferrojan horse" phage infections in aquatic environments for the community assembly. Similar to K02014, hsdR and hsdM (K01153 and K03427) that encode restriction enzymes interacting with phages, exhibited very high allele variability in LMO, SOLA and MEN.
Correlations of auxiliary gene variability with seasonal metacommunity size (log-transformed) or gene copy number (logtransformed) in the five metagenome time-series. Blue open data points and the blue trend line indicate cytoplasmatic genes, orange filled data points and the orange trend line indicate noncytoplasmatic genes, and light grey open data points indicate genes for which subcellular location could not be determined following our working definition. The black trend line represents all data points independent of their subcellular location. The highlighted enlarged data point indicates a gene orthologue encoding an outer membrane iron receptor protein (K02014), which in all time-series exhibited maximal temporal variability. All plotted data points present mean values obtained after 1,000 permutations repeating the subsampling processes number range covered, sequencing depth and data processing; Tables.
(CSP 1977) , which generated the Lake Mendota time-series. K. D. M. received funding from the  United States National Science Foundation Microbial Observatories Program  ( MCB-0702395 ), the  Long-Term Ecological Research Program  ( NTL-LTER DEB-1440297 ) and an  INSPIRE Award  ( DEB-1344254 ). P. G. was supported by the  Agence Nationale de la Recherche (ANR)  through the project  EUREKA  ( ANR-14-CE02-0004-01 ). Open access funding enabled and organized by Projekt DEAL.
Funding information Deutsche Forschungsgemeinschaft , Grant/ Award Number:  BE5937/2-1 ;  United States National Science Foundation , Grant/Award Number:  MCB-0702395 ;  Joint Genome Institute , Grant/Award Number:  CSP1977 ;  Agence Nationale de la Recherche , Grant/ Award Number:  ANR-14-CE02-0004-01 ; Long- Term Ecological Research Program , Grant/Award Number:  NTL-LTER DEB-1440297 ;  INSPIRE Award, Grant/Award Number :  DEB-1344254 ;  Swedish Research Council (VR)  BioProject  PRJNA273799 ,  JGI IMG/M study ID: Gs0063444 ).
E 1 1 Correlations of auxiliary gene variability with seasonal metacommunity size (log-transformed) or gene copy number (logtransformed) in the five metagenome time-series. Blue open data points and the blue trend line indicate cytoplasmatic genes, orange filled data points and the orange trend line indicate noncytoplasmatic genes, and light grey open data points indicate genes for which subcellular location could not be determined following our working definition. The black trend line represents all data points independent of their subcellular location. The highlighted enlarged data point indicates a gene orthologue encoding an outer membrane iron receptor protein (K02014), which in all time-series exhibited maximal temporal variability. All plotted data points present mean values obtained after 1,000 permutations repeating the subsampling processes number range covered, sequencing depth and data processing; Tables.
5 simulations with ESMs participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5) to 171 ± 12 Gt C yr -1 , compared to the unconstrained model range of 156-247 Gt C yr -1. In a second step, we use a machine learning model to constrain gridded future absolute GPP and gridded fractional GPP change in two independent approaches. For this, observational data are fed into the machine learning algorithm that has been trained on CMIP5 data to learn relationships between present-day physically relevant diagnostics and the target variable. In a leave-one-model-out cross-validation approach, the machine learning model shows superior performance to the CMIP5 ensemble mean. Our approach predicts an increased GPP change in northern high latitudes compared to regions closer to the equator.
Moreover, instead of considering a single globally or regionally averaged value per climate model, our approach uses spatiotemporally gridded data for every ESM in the ensemble. This dramatically increases the available number of training data (even if they are not statistically independent) and uniquely enables the statistical model to identify and exploit regional processes. Furthermore, the large training database allows us to perform extensive out-of-sample testing, which is crucial to ensure the model robustness. In the second phase (predicting with the trained GBRT model), we use observational data products to constrain the RCP 8. 5 CMIP5 projection of GPP. To do this, observational data for all predictors x of the historical climate (for every grid cell/month) are fed into the trained GBRT model in order to get a spatially constrained projection of GPP (see black dots in Figure  2b ).
The raw CMIP5 and observation-driven data sets are preprocessed before entering the machine learning algorithm (see Text S1). In particular, increases greater than 300% in the fractional GPP change that is used as target variable in Step 2b are removed. These unrealistically high values (making up about 5. 9% of the data) occur in places with very small absolute GPP values in the historical simulation that are used in the denominator of the derivation of the fractional GPP change (see Equation  1 ). In addition, we randomly select 25% of the training data and use it as hold-out test data. This data set is excluded from the training and can therefore be used for an independent validation of our statistical model.
Observation-Based GBRT Prediction of Absolute GPP In this section, we use the different statistical models to predict absolute GPP at the end of the 21st century. For this, we feed observation-based data (see Table  1 ) into the regression models GBRT, GBRT-1D, LASSO, and LASSO-1D. The result of Step 1 on the gridded data is illustrated in Figure  5a , which shows the ratio between rMMM and MMM. This ratio is almost constant over the whole globe, which is not surprising due to the global nature of the rescaling (see Equation  2 ). The use of the global emergent constraint predicts a slightly lower GPP increase over the 21st century (39%) than the unweighted CMIP5 ensemble mean (43%). The corresponding ratio 39 % /43 % ≈ 0. 91 is approximately equal to the mean value of rMMM/MMM (0. 92). Consequently, the global estimate of rMMM gives a lower GPP at the end of the century (179 Gt C yr -1 ) as MMM (198 Gt C yr -1 ), which can be interpreted as a correction of the ESMs' overestimation of future GPP to changes in the atmospheric CO 2 concentration.
These issues could be tested in the future by including other (even offline) models, which account for the aforementioned processes. A true validation will only be possible once we experience the changed condition or could be based on paleo data. In summary, our approach does not address the coupled system directly but hinges on the climate predictions of the individual ESMs and hence is not able to implicitly correct respective biases of the CMIP5 ensemble. The presented approach based on machine learning is not limited to projections of the future carbon cycle, which was used as an illustrative example. Indeed, since its only prerequisites are the availability of gridded climate model data and gridded observation-based data, the GBRT algorithm can be applied to any variable of interest if physically relevant diagnostics that influence the target variable are known.
This opens a wide range of possibilities for constraining uncertainties in projected variable scenarios. Additionally, maps of the local relative feature importance as shown in Figure  7c  can be used to reveal connections between different Earth system variables that are currently unknown. With the proposed data-driven approach, a possibly extended application could be to constrain projections of not just one but several variables simultaneously. Moreover, the concept of our method is not limited to GBRT but can be used with any other (machine learning) regression algorithm. 10. 5281/zenodo. 3387139. ESMValTool and ESMValCore are developed on the GitHub repositories available at https://github. com/ESMValGroup. All scripts are implemented in Python and use the following packages for their core functionalities: sklearn  (Pedregosa et al.
, 2011) , xgboost  (Chen & Guestrin, 2016) , and lime  (Ribeiro et al. , 2016). CMIP5 model output (see Table  S1 ) is available through the official nodes of the Earth System Grid Foundation (ESGF) and can be directly used within the ESMValTool (e. g. , https://esgfdata. dkrz. de/projects/esgf-dkrz/). Download instructions and preprocessing scripts for observation-driven products (CRU, ERA-Interim, LAI3g, FLUXNET-MTE, and Scripps-CO2-KUM) are included in the ESMValTool.