dataesr/scientific-paragraphs-categorization huggingface dataset:
Control beads (100 and 400 nm) were supplied by Malvern Instruments. EV lysates were prepared in RIPA buffer containing a protease inhibitor cocktail (Roche Applied Science) and quantified using the Bradford method. Next, the lysates were subjected to SDS-PAGE and transferred to nitrocellulose membranes (Whatman). Membranes were then incubated with the primary antibodies overnight at 4°C. The following primary antibodies were used: calnexin (clone C4731, Sigma-Aldrich; 1:2000 dilution), TSG101 (ab125011, Abcam; 1:1000 dilution) and CD63 (ab59479, Abcam; 1:1000 dilution). After incubation with the specific HRP-conjugated antibody (Vector; 1:10,000 dilution), the chemiluminescent signal was detected using Clarity and/or Clarity Max (Bio-Rad, Italy). The autoradiographic films thus obtained were quantified using ImageJ2 software.
In order to exclude that protein-miRNA complexes could co-precipitate with EVs during the isolation procedure, we used a combination of density gradient separation and treatment of isolated vesicles with proteinase K 20 μg/ml to confirm that miR-21-5p is truly loaded into EVs. MiR-21-5p levels were not significantly different between the two conditions (miR-21-5p relative expression; control = 0. 383 ± 0. 050; proteinase K = 0. 259 ± 0. 063; n = 3; p = 0. 185), confirming the intra-vesicular presence of the miRNAs. Data regarding the EV-related aspects of the experiments conducted on HUVECs have been deposited to the EV-TRACK knowledge base (http://evtrack. org) [44]  under accession number EV190105.
Size selection allowed keeping 130-160 bp fragments. After ethanol precipitation, the library pool was quantified with Agilent High Sensitivity DNA Kit, diluted to 1. 8 pM and sequenced using NextSeq® 500/550 High Output Kit v2 (75 cycles) (Illumina; FC-404-2005) on the Illumina NextSeq500 platform. Raw base-call data generated by the Illumina NextSeq 500 system were demultiplexed using Illumina BaseSpace Sequence Hub (https://basespace. illumina. com/home/ index) and converted to FASTQ format. After a quality check with FastQC (https://www. bioinformatics. babra ham. ac. uk/projects/fastqc/), the adapter sequences were trimmed using Cutadapt (http://cutadapt. readthedocs. io/ en/stable/index. html), which also removed sequences < 10 nucleotides. Reads were mapped using the STAR algorithm (https://www.
ncbi. nlm. nih. gov/pubmed/ 23104886). The reference genome consisted of human miRNA sequences from the miRbase 21 database  [45]. Raw counts from mapped reads were obtained using the htseq-count script from the HTSeq tools (http://wwwhuber. embl. de/HTSeq/doc/overview. html). Raw counts were further normalised using DESeq2 bioconductor package (http://bioconductor. org/packages/ release/bioc/html/DESeq2. html). NGS raw data are available through European Nucleotide Archive (http://www. ebi. ac. uk/) under accession number PRJEB33703. MiRNAs expressed in SEN cells and their EVs were identified to obtain a normalised expression above the 40th percentile in at least one sample of the group.
Statistical analysis of NGS data NGS data were analysed using Genespring GX software v14. 8 (Agilent Technologies). The miRNAs showing a differential expression were identified by comparing SEN and CON cells using a fold change ≥ 1. 5 filter and FDR 5% at moderated t-test with Benjamini-Hochberg correction. A fold change ≥ 1. 3 with FDR 20% or 30% was used in lEV and sEV analysis, respectively. Cluster analysis was performed using GeneSpring GX software with Manhattan correlation as a similarity measure.
Methylome analysis of loci coding for SA miRNAs The sequences and genomic coordinates (GRCh38) of predicted hairpin precursor miRNAs were retrieved from miRbase 21 (Supplementary Table  7 ). The CpG nucleotide clusters (CpG islands) associated with miRNA genes within a distance of 2 kb were identified by visual inspection of the miRNA gene neighbourhood using the University of California Santa Cruz (UCSC) Genome Browser. Putative proximal promoter regions within 2 kb of the TSS were obtained from the database http://bicresources. jcbose. ac. in/zhumur/pirna quest/.
After a wash in PBS, cells were stained with nuclear  HOECHST 33, 342 (cat. no. H-3570; Molecular Probes, USA) at 1:10,000 dilution in PBS for 5 min, coverslipped with Vectashield mounting medium (H-1200, Vector Laboratories) and viewed in a fluorescence microscope. Omission of the primary antibody resulted in lack of labelling, confirming the specificity of the antibody. The presence of senescenceassociated heterochromatin foci (SAHF) was assessed by measuring the coefficient of variation (CV) of the HOECHST staining intensity in each nucleus. At least 200 nuclei have been assessed for each condition and each replicate using the using the CellProfiler image analysis software, version 3. 1. 9 [48]. The Ki67 proliferative index was expressed as a percentage based on the number of Ki67-positive cells.
These data showed that in elderly people circulating EVs are characterised by increased miR-21-5p levels compared with either younger adults and healthy centenarians, suggesting that the EV miR-21-5p content can be a marker of cellular senescence and inflammaging. This finding lends support to our group's earlier report of an inverse U-shaped age-related trend of miR-21-5p purified from total plasma  [57]. According to IPA, SIRT1 and miR-21-5p are spokes in a network that includes miRNA loading on factor 1 with a score > 0. 60 or < -0. 60. TNF and TP53 can be considered as hubs. The association of the miRNAs loaded on factor 1 with the IPA database "diseases and disorders" revealed a significant enrichment of the sub-category "inflammation of organ", confirming the hypothesis that an inverse U-shaped trend reflects high circulating levels of the biomarkers of inflammaging in elderly individuals (Figure  7(d) ).
Data from four subjects in each age group are displayed as mean ± SD. *, p < 0. 05; **, p < 0. 01; ***, p < 0. 001 from one-way ANOVA followed by Tukey's multiple comparison test for pairwise comparisons. (d) Ingenuity Pathway Analysis (IPA) on miRNAs loading on factor 1 with a score > 0. 60 or < -0. 60. Association of this dataset with the "diseases and disorders" IPA database revealed a significant enrichment of the subcategory "inflammation of organ".
| Metabarcoding data Metabarcoding data for Tara Oceans samples were exploited for the present study. Total nucleic acids (DNA + RNA) were extracted from all the samples, and the hypervariable V9 region of the nuclear 18S ribosomal DNA was amplified through PCR (see  Alberti et al. , 2017; de Vargas et al. , 2015). A quality filtering based on reads quality checks and a minimum number of occurrences of three copies in at least two different samples was implemented to reduce PCR and sequencing errors. Within this dataset, 237,565 V9 diatom-assigned ribotypes were detected  (de Vargas et al. , 2015). We focus on the 20-180 μm size fraction dataset, which contains 183 net samples encompassing 125 stations sampled at the subsurface (5 m) and 58 of them at deep chlorophyll maximum depth (DCM).
| Bioinformatics pipeline The initial taxonomic assignation of reads was confirmed and refined using a custom version of the Protist Ribosomal Reference database (PR 2 ;  Guillou et al. , 2013)   showing similarity >90% over a query coverage with the reference > 109 bp. Reference and environmental sequences were aligned with mafft v. 7  (Katoh, Rozewicki & Yamada, 2017)  using the experimental service for large numbers of highly similar and short sequences (https://mafft. cbrc. jp/align ment/serve r/large. html?aug31). Poorly aligned sequences and hard-to-align blocks within the alignment were removed. Ribotypes taxonomically annotated to diatoms were furthermore filtered through a phylogenetic approach. Phylogenetic analyses were performed using the approximately-maximum likelihood method  (Yang, 1994)  implemented in the fasttree2 software  (Price, Dehal & Arkin, 2010).
| Morphology-based data Morphology-based counting was implemented for subsamples of the same size fraction used for molecular analyses (20-180 μm), obtained from a few to more than 100 L of seawater, at surface and DCM depth from the stations in the Cape Agulhas region, the South Atlantic transect, South Pacific Ocean and the Southern Ocean  (Malviya et al. , 2016). Additional samples were analysed from the Atlantic, Indian and Pacific Oceans, the Mediterranean Sea and the Arctic. Up to 3 mL of each net sample was placed in an Utermöhl chamber. The whole sedimentation chamber bottom was observed and cells (> 5 µm) were identified up to the species level whenever possible using a light inverted microscope (Axiophot200, Carl Zeiss, Oberkochen, Germany) at 400 × magnification  (Utermöhl, 1958).
A false discovery rate p-value adjustment of the p-value was implemented. The optimal correlation between morphology-based richness and any barcode-based richness at different thresholds was selected as the most significant correlation (adjusted p-value <. 05) with the highest correlation coefficient (rho). This optimal correlation drove the choice of the molecular data type other than of the filtering threshold to be applied to the whole dataset. This latter filtered information will be used for all downstream analysis in this paper. The filtering process is schematized in Supporting Information Figure  S1. The cumulative threshold globally applied corresponded to a different absolute threshold for each station, computed as the maximum relative abundance of the discarded OTUs of each station.
| Environmental data Nine environmental variables were considered to investigate the processes behind diatom richness dynamics. Variables included descriptors of hydrodynamic mixing, nutrient availability, temperature and chlorophyll a concentration. Local confluence and mixing were estimated through the finite-size Lyapunov exponent, a measure of the front intensification rate computed as the backward-in-time relative separation of water parcels from altimetry-derived surface currents  (Lehahn, d'Ovidio & Kohen, 2018 ). An altimetry-based three-days advection scheme was also applied to 25-km resolution infrared sea surface temperature (SST) images to estimate the SST gradients at km scales at sampling sites  (d'Ovidio, De Monte, Alvain, Dandonneau & Levy, 2010). Temperature, salinity, silicate, phosphate and nitrate availability were extracted at 5-m depth from the World Ocean Atlas 2013 (WOA13) database, while iron availability was derived using the pelagic interactions scheme for carbon and ecosystem studies (PISCES)-v2 model  (Aumont, Ethé, Tagliabue, Bopp & Gehlen, 2015).
| Earth system models WOA and PISCES-v2 data were used to obtain present-day global environmental conditions. Future environmental conditions were derived by computing the mean of six Earth system models over the 2006-2015 (present-day) and 2090-2099 (end of the century) decades under the greenhouse gas emission scenario RCP8. 5. The models included in the analysis are IPSL-CM5A-LR/MR, GFDL-ESM2G/M, MPI-ESM-LR/MR, CESM1-BGC, HadGEM2-ES and NorESM1. We considered the average value of the models as it tends to smooth errors and incongruences between models  (Bopp et al. , 2013). The delta between the end of the century and the present-day mean models was computed (Supporting Information Figure  S2 ) and added to the WOA and PISCES v2 global scale data to obtain future environmental conditions.
| Machine learning modelling integration In order to model surface diatom richness, as measured by the optimally filtered molecular information, we integrated four machine learning approaches: boosted regression tree (BRT), random forest (RF), fully connected neural network (NN) and generalized additive models (GAM). The parameterization of each model was optimized (Supporting Information Text S1). The final models were then trained on the whole dataset. Model performance was investigated by computing the mean residuals and the Pearson correlations between the observed and predicted richnesses for the whole dataset (Supporting Information Table  S2 ). Models with a cross-validated root-mean-square deviation below 50 and mean residuals below 30 were selected as significant. Two sets of models were built and optimized: one including all the nine variables previously described to investigate the roles of the variables on the model and a second one excluding the hydrodynamic parameters (SST gradient and Lyapunov exponent) to allow the projection at global scale, as we lack global scale data for these two variables.
In both cases, three out of four models were considered as significant, excluding the GAM model from both exercises (Supporting Information Table  S2 ). The projection of the seven-variable models was performed at the global scale using both the present time and the end of the century conditions. The variable importance  (Fisher, Rudin & Dominici, 2018)  was computed for the nine-variable models using the DALEX R package  (Biecek, 2018)  as the difference between the loss function calculated for validation data with every variable being shuffled and the loss function calculated for the validation dataset. Furthermore, with the iBreakDown R package  (Gosiewska & Biecek, 2019)  we computed the variables contribution to single predictions and we averaged it across models. We used the ninevariable models to evaluate Tara Oceans stations predictions and the seven-variable models to evaluate global scale predictions.
Using the seven-variable models, the local variable contributions were computed on global scale data seven more times (one per variable), using each time the future conditions of all the parameters except one, which kept present-day condition. Finally, the difference between the future-present-day variation of richness using all the future conditions and the future-present-day variation of richness using as future conditions all the future parameters except one was computed. All data-mining and statistical analyses were performed in R (version 3. 4. 1; R Core Team, 2017).
| Integration of morphological and metabarcoding counts Our first step was to strengthen the consistency among diatom richness estimates from molecular and morphological data. There have been previous efforts to compare metabarcoding-and morphology-based information for diatoms (e. g. ,  Malviya et al. , 2016; Piredda et al. , 2018; Zimmermann, et al. , 2015). However, while previous studies limited their analysis to a descriptive comparison check, here we follow a different approach, proposing a reproducible procedure to reasonably reconcile the two measures. The Tara Oceans metabarcode dataset, following a stringent quality filtering, contains 174,267 ribotypes annotated as Bacillariophyta. This study focuses on diatoms from 183 sampling stations addressing the 20-180 µm size fraction, because for those net samples morphology-based data adequately cover richness, thus allowing a reliable comparison with metabarcoding data.
, 2018). From the quantitative viewpoint, a remarkable match of relative abundances between morphological and metabarcoding information has been observed for diatoms  (Malviya et al. , 2016; Piredda et al. , 2016 Piredda et al. , , 2018)) , although this is not usually the case for protists  (Abad et al. , 2016; Massana et al. , 2015). The rationale of the approach herein proposed was to obtain a sufficient covariance between the results of the two methods in order to have a similar spatial pattern of richness while still retaining the high resolution of the metabarcoding. In addition, this will allow comparison with patterns based on microscopic counts (e. g. , Righetti, Vogt, Gruber, Psomas  & Zimmermann, 2019). With these goals, the datasets obtained through different clustering algorithms were progressively filtered (see Materials and methods, Supporting Information Figure  S1 ), and the resulting richness correlated to the morphology data (Figure  1a , Supporting Information Figure  S3a  total abundance of OTUs.
| The identity of discarded OTUs Under the hypothesis that filtering mainly discarded genetic variants of more abundant haplotypes, we investigated the nature of the excluded OTUs by assigning both retained and discarded OTUs to named species. Only 24% of the OTUs filtered out belonged to species that were not detected in the retained dataset, while the vast majority of them were presumably variants of more abundant haplotypes at the same station (Figure  2 ). Further support for this hypothesis comes from the strong relationship between the abundance of diatoms in a sample and the number of discarded OTUs (Figure  3 ). This relationship suggests a higher intraspecific variability in diatom-rich stations, which can be explained by (a) a lower detectability of variants in samples with low abundance of diatom reads  (Elbrecht, Vamos, Steinke & Leese, 2018) , (b) high intraspecific diversity in conditions favourable to diatom growth as a consequence of cell proliferation, sexual reproduction and resting stage germination  (Godhe & Rynearson, 2017; Lebret, Kritzberg, Figueroa & Rengefors, 2012) , or both.
As this correction did not substantially impact the results (Supporting Information Figure  S5b ), downstream analyses have been performed on uncorrected filtered richness, so as to keep a better comparability with morphology-based datasets. We investigated the spatial patterns of diatom richness at the global scale (Figure  4 ). As designed, the filtering procedure only F I G U R E 2 Overview of the identity and distribution of the filtering discards. The boxplots show the discarded operational taxonomic units (OTUs) percentage over the total number of diatom OTUs discarded. The plot aggregates this information calculated in every station and then divided in three classes. The three classes are: (a) the discarded OTUs assigned to species still represented in a station (green), (b) discarded OTUs assigned to species that are not present in the retained set in a given station (yellow) and (c) non-annotated discarded OTUs (grey) [Colour figure can be viewed at wileyonlinelibrary.
com] slightly lowered the metabarcoding resolution but allowed the richness pattern to be reconciled with that obtained from morphologybased data. Indeed, both sets of data showed higher richness in the Equatorial Pacific and the Agulhas regions (Figure  4b, c ), while the peaks of richness that were molecularly identified at the poles by the unfiltered metabarcoding data (Figure  4a ) are now bevelled to the medium-low values observed in the morphology-based data. To extend the distribution of diatom richness calculated from the Tara Oceans stations to the global scale, we applied an integrated machine learning modelling approach, using as predictor variables physical parameters (temperature and salinity), chemical parameters (nitrate, iron, silicate, phosphate) and a proxy of the trophic status of the system (chlorophyll a; see Materials and methods).
| Towards a unified picture The comparison between metabarcode-and morphology-based diatom richness estimates in the global Tara Oceans survey displayed over a 20:1 ratio between the former and the latter. Such imbalance highlights the need to clearly specify which kind of richness is herein being studied. Since our goal was to characterize spatial patterns of diatom richness globally, we decided to integrate both types of information assuming that morphological and molecular data could not be uncorrelated. To that aim, we selected the threshold below which the information provided by the two estimates was consistent: this allowed us to include most of the rare biosphere and to reproduce the expected unimodal relationship between diatom abundance and observed richness. This approach is not purely methodological since our filtering procedure allowed a large number of the OTUs to be kept, those above the empirically-found threshold of abundance, which are cumulatively structured in space very much like the canonical, morphology-based species.
This, in turn, suggests that this reconciled OTU dataset follows the same ecological patterns of the canonical species on which our knowledge on diatom ecology is largely based, while being one order of magnitude richer. From an evolutionary perspective, species richness may result either from neutral or adaptive processes. While the former could pass undetected by microscopy counts, the latter might lead to different morphotypes and thus increase also microscopy-based richness. We posit that microscopic counts always tend to underestimate or miss the rare biosphere. Furthermore, even without equating OTUs to morphotypes, the OTUs' distribution is a better representation of the diversity of a sample, in the sense of the number of distinct taxonomic units, than microscopic counts. Global estimated patterns of phytoplankton richness, although not numerous, have been obtained by model simulations (e.
g. ,  Barton et al. , 2010) , or from direct observations (e. g. ,  Righetti et al. , 2019; or Ibarbalz et al. , 2019  based on the same dataset). Our spatial patterns are based on global-scale data using different machine learning methods that consider not only the intrinsic structure of spatial variance, as in geostatistical methods, but also other variables that might affect the richness, but without assuming any driving mechanism as in biogeochemical models. Our maps mimic the global richness patterns of phytoplankton distributions obtained by biogeochemical models (e. g. ,  Barton et al. , 2010) , albeit with a few discrepancies. We interpret some of them as due to inadequate coverage of Tara Oceans sampling across specific regions, particularly when those regions have characteristics not present elsewhere (e.
| CON CLUDING REMARK S Machine learning tools are useful to integrate extensive, multivariate datasets to set the starting point of a mechanistic interpretation of patterns that should complement the evidence provided by them. Molecular data strongly suggest that diatoms are less exposed to a temperature-dependent latitudinal decrease in richness and, supported also by microscopic counts, 'the variation of richness with latitude deviates from the pattern of most plankton. Overall, our analysis confirms and reinforces that even for a single phytoplankton group, the apparent paradox of many species coexisting in the same water parcel holds true  (Hutchinson, 1961). There is no unique driver of such patterns, with non-equilibrium (first scenario), hump-shaped richness-productivity relationship (second scenario), and latitudinal dependence (third scenario) always concurring, albeit with different geographical patterns, to modulate richness.
(c) Reconciled richness pattern derived by the filtering of the Swarm metabarcoding in sizeclass 20-180 μm optimized over the morphological observations. All maps employ the Mollweide's projection [Colour figure can be viewed at wileyonlinelibrary. com] Global prediction of diatom richness. (a) World scale map (Mollweide's projection) of predicted present-day annual mean diatom richness computed at each grid point as the average of the predicted richness of the random forest (RF) model. Present-day environmental conditions are provided by the World Ocean Atlas database with the addition of pelagic interactions scheme for carbon and ecosystem studies-v2 biogeochemical model information for iron. Panel (b) shows the latitudinal distribution of the same measure of diatom richness depicted in panel (a).