dataesr/scientific-paragraphs-categorization huggingface dataset:
Counts from expressed genes were rst normalised with the calcNormFactors() function  [32]. Then, voom() function from limma R package (v. 3. 50. 0) was used to t a generalised linear regression model to correct the data with the group as a xed effect  [33]. The p-values were adjusted for multiple testing using the Benjamini and Hochberg procedure  [34]. Only DE miRNAs with an adjusted P-value < 0. 05 were used for the downstream pathway analysis.
In silico functional analysis Using only DE miRNAs, a predictive functional analysis was performed to evaluate the in uence of the DE miRNAs on biological processes and pathway. In silico analyses were performed using human orthologs (hsa-), instead of bovine orthologs (bta-), since information about miRNAs activity is more detailed and extensive in human than bovine. In this particular step, we used miRWalk 2. 0 (http://mirwalk. umm. uniheidelberg. de), OmicsNet 2. 0 (https://www. omicsnet. ca) and Cytoscape v. 3. 9. 1 with ClueGO plugin v. 2. 5. 9 (https://cytoscape. org) softwares. First, the most deposited mature form (-3p or -5p) of DE miRNAs were checked in miRBase database and modi ed accordingly in the DE miRNA list with the aim of retrieving the most relevant results with the functional analysis.
Subsequently, the partially modi ed DE miRNAs list was used in miRWalk to retrieve all the miRNA-gene interactions. In particular, only validated interactions deposited in miRTarBase database and with biding position data (related to 3UTR, 5UTR and CDS) were selected. In parallel, to identify functional biological categories related to DE miRNAs, Panther Biological Processes analysis was performed using OmicsNet. Then, the categorization of biological processes was manually performed. Biological processes were clustered in two main macro-categories: cell life processes (including cell cycle, regulation of cell cycle, cell proliferation and apoptotic process) and gene expression machinery processes (including regulation of transcription by RNA polymerase II, Transcription DNA-templated, mRNA processing, mRNA splicing via splicesome, and regulation of translation and protein phosphorylation).
Annotations and differential analyses In total, 2127 Bos taurus annotated miRNAs (bta-miRNAs) were detected. Using a cut-off value for the SCC parameter of 200,000 cells/ml, a total of 1997 miRNAs were differentially expressed (DE) between group L and group H. Among them, 1267 were upregulated and 730 downregulated including 68 miRNAs with a false discovery rate (FDR) < 0. 05. Further investigations were performed using these 68 DE miRNAs. They were mostly downregulated. Intriguingly, only 1 miRNA, bta-miR-361-3p, was upregulated with 2. 7-fold change (FC) value. Among the downregulated miRNAs, a total of 15 had a FC < -2. The complete list of DE miRNAs with related FC and p-value data is reported in Additional le 1.
Functionalanalysis The functional analysis was conducted using human orthologues (hsa-) instead of bovine (bta-), since human databases are more complete and more informative regarding miRNA function. The targets for 17 known miRNAs were obtained from OmicsNet with Panther database. In addition, the regulation of immune system processes was investigated using ClueGO plugin in Cytoscape. The main immune processes signi cantly in uenced by the DE miRNAs were mainly related to mediated immunity including immunoglobulin productions and lymphocyte regulation. Therefore, the miRNA-gene networks were displayed using Cytoscape and three separated reactomes were generated considering cell life (Fig. 2 ), gene expression (Fig. 3 ) and immunity processes (Fig. 4 ). Two large nodes of regulation of miRNA-gene interactions involved miR-503-5p and miR-455-3p both cell life and gene expression processes.
The abovementioned studies investigated miRNAs role during different mastitis models working with mammary gland biopsies  [35, 37, 38] , cultured Mac-T cells  [19]  or blood  [36]. Therefore, an adequate comparison is not really possible. Moreover, in the dataset of this study only two samples were positive to S. uberis and this may not be enough to highlight differences between the groups. Differently, Tzelos and colleagues investigated the expression of miR-223 during mastitis, but a statistically signi cant different expression was not observed between mastitic and healthy cows  [40]. In addition, in Tzelos's study miR-223 was investigated both on whole and skim milk samples and not on milk EVs affecting also in this case the possibility of an adequate comparison. However, according to human medicine literature, miR-223 is either expressed almost exclusively or highly enriched in several subsets of white blood cells including T-and B-lymphocytes, neutrophils, mast cells and monocytes  [41].
3b ). Therefore, we deem pyroclastic flow emplacement and caldera collapse as unlikely tsunami source mechanisms. The existence of a 2500-m-wide and 500-m-deep crater, however, provides clear evidence for a highly explosive eruption phase. The seismic and bathymetric datasets suggest a link between this explosion and the flank deformation. At the top of the remnant cone, the detachment surface intersects prominent ridges that are the seafloor expression of feeder dykes (Fig. 4a, d ). Their present exposure is probably due to their higher mechanical stability relative to the surrounding volcanoclastic deposits, which have been eroded since the 1650 eruption. While the explosion represents an obvious tsunami source mechanism, the previously unconsidered tsunamigenic potential of the flank deformation requires careful evaluation.
Chain of events during the 1650 eruption The combination of high-resolution marine geophysical datasets and historical eyewitness accounts offers a unique opportunity to reconstruct the sequence of volcanic processes that triggered the 1650 explosion and tsunami. The 1650 cone consists of highly vesicular pumice, which was deposited as fallout from the eruption column, where many of the large pumice clasts floated at the sea surface before becoming water-saturated and sinking  18  (Fig. 7a, b ). The 1650 cone had a volume of approximately 4 to 6 km³ (ref. 16), which was deposited in a very short timeframe of only 2 weeks, based on the onset of visually detectable eruption products in eyewitness accounts  18, 20, 21. This rapid deposition has likely precluded significant consolidation of the pumice deposits covering the pre-1650 topography.
Toward improved flank stability monitoring and volcanogenic tsunami early-warning systems The 3D seismic analysis of Kolumbo highlights the complexity of tsunami genesis during volcanic eruptions, with flank deformation being the nucleus of a chain of events in a critical system. Internal deformation within volcanic flanks can progress over long periods, as seen by 3D seismic data of Ritter Island's failed flank 8 and pre-collapse InSAR data for Anak Krakatau  7. The most prominent historical sector collapse triggered the 1980 eruption of Mount St. Helens, which was preceded by intrusion-controlled flank deformation. Although not all deformation of a volcanic flank will lead to a catastrophic collapse, deformation was a key precursor to these historical volcanic collapse events. Current monitoring strategies targeting the detection of flank deformation are based on satellite data and repeated photogrammetric imaging  7  , methods that cannot be used to survey submerged volcanic edifices.
3D seismic survey The high-resolution P-Cable 3D seismic dataset was collected during research cruise POS538 aboard R/V Poseidon in October 2019 (ref. 46). The receiver array consisted of 16 streamers (each with 8 channels, spaced 1. 5625 m apart) attached to a cross cable of ~200 m length. The source was a 45in 3 -45in 3 generator-injector airgun (harmonic mode), fired every 4 s, with a nominal ship speed of 3. 5 knots through water. Processing included geometry definition, trace editing (removing dead and very noisy traces), static time corrections, anomalous amplitude attenuation, frequency-wavenumber (FK) filtering and trace interpolation. We then applied adaptive subtraction multiple suppression on common receiver gathers using a 1D pattern-matching algorithm. The resulting data were sorted to common midpoints on a grid with cell sizes of 3.
Numerical tsunami simulations We used the numerical simulations code VolcFlow for all presented simulations (details in refs. 26-29,47). VolcFlow simulates the slide as a fluid and couples the interaction between the slide and the water based on a depth-averaged approximation on a topography-linked coordinate system. VolcFlow applies a plastic behavior with turbulence as rheology, which is primarily controlled by yield strength and density. The code uses general depth-averaged equations of mass and momentum conservation  43, 48. The slide was defined by reconstructing a pre-slide topography (Fig. 5 ), which was achieved by using the geometry of the remaining cone and constraints about the pre-slide geometry (i. e. , that the cone breached the sea surface) and mapping the detachment surface in the 3D seismic datasets as slide plane, which had to be converted into the time domain using a seismic velocity of 1.
Autonomous underwater vehicle (AUV) mapping During cruise POS510 (ref. 44), detailed bathymetric data were acquired through 7 planned dives by AUV Abyss, which is equipped with the dual frequency (200/400 kHz) RESON Seabat 7125 multibeam. During the acquisition, the nominal frequency of 200 kHz was applied, and all the logged data were further processed using the open-source software package MB-System. First, the navigation was adjusted using the interactive graphical interface program MBnavadjust and then merged with the bathymetric data. The following step included identifying and cleaning erroneous beams; therefore, dedicated filters were applied and manual cleaning through 2D and 3D editors. Finally, the processed bathymetric data were gridded, yielding a spatial resolution of 2 m, using the function mbgrid and, more specifically, the Gaussian weighted mean algorithm.
Acknowledgements We thank the masters and crews of  RV Poseidon  for their support during research cruise POS538. We thank  Mark Hannington  and  Sven Petersen  for providing access to the AUV bathymetry dataset. Thanks also to  Dirk Klaeschen  for giving valuable advice at several key stages of seismic processing. We also thank  Schlumberger  and  IHS  for granting educational licenses to their software.
The swallowtail genus Papilio (Lepidoptera: Papilionidae) is species rich, distributed worldwide, and has broad morphological habits and ecological niches. Because of its elevated species richness, it has been historically difficult to reconstruct a densely sampled phylogeny for this clade. Here we provide a taxonomic working list for the genus, resulting in 235 Papilio species, and assemble a molecular dataset of seven gene fragments representing ca. 80% of the currently described diversity. Phylogenetic analyses reconstructed a robust tree with highly supported relationships within subgenera, although a few nodes in the early history of the Old World Papilio remain unresolved. Contrasting with previous results, we found that Papilio alexanor is sister to all Old World Papilio and that the subgenus Eleppone is no longer monotypic.
The molecular data mainly came from previous studies (e. g. Zakharov et al. , 2004a Zakharov et al. , , 2004b;; Condamine et al. , 2012 Condamine et al. , , 2013a Condamine et al. , , 2013b;; Lewis et al. , 2015; Wu et al. , 2015; Owens et al. , 2017 Owens et al. , , 2020;; Allio et al. , 2021; Joshi and Kunte, 2022). We used five mitochondrial gene fragments (COI, COII, ND1, ND5 and rRNA 16S) and two nuclear gene fragments (EF-1a and Wg), chosen based on their availability among Papilio species. We aligned the DNA sequences for the rRNA 16S using MAFFT 7. 110  (Katoh and Standley, 2013)  with the E-INS-i algorithm, while we aligned the coding genes using MACSE 2. 00  (Ranwez et al. , 2011)  with the alignSequences subprogram and default options. All the resulting alignments were checked for codon stops and eventually refined by eye with Mesquite 3.
Inferring phylogenetic relationships We performed both maximum likelihood (ML) and Bayesian Inference (BI) to reconstruct phylogenetic relationships. Although phylogenomic approaches have used amino acids on Papilionidae at the genus level  (Allio et al. , 2020) , all our analyses relied on nucleotides for the tree inference of Papilio because the molecular dataset is limited to seven gene fragments and converting it to amino acids will result in a smaller dataset with limited phylogenetic information given the species diversity of the focal clade. ML inference was implemented with IQ-TREE 2. 1. 2 (Minh et al. , 2020)  using ModelFinder to select the best-fit partition scheme and the best-fitting substitution model for each partition (-m MFP + MERGE option,  Chernomor et al. , 2016; Kalyaanamoorthy et al.
, 2018). BS values and UFBS values were considered strong when higher than 70% and 95%, respectively. We also assessed branch support using Shimodaira-Hasegawa approximate likelihood ratio test (SH-aLRT;  Guindon et al. , 2010) , with SH-aLRT values above 80% considered as strong support for a clade. BI analyses were performed using MrBayes 3. 2. 7a (Ronquist et al. , 2012). Because ModelFinder includes more substitution models than BEAST or MrBayes, the most likely partitioning scheme and substitution models can be altered, the molecular dataset was analyzed with Parti-tionFinder 2. 1. 1  (Lanfear et al. , 2017)  to estimate the best partition scheme with an initial subset of 19 possible partitions (all three codon positions for coding sequences and non-coding genes treated separately). Partitions and corresponding optimal substitution models were searched using the greedy algorithm, the mrbayes set of models, and the Bayesian Information Criterion (BIC) to compare the fit of different models.
7. 1  (Rambaut et al. , 2018). The runs had to have values of ADSF approaching zero, PSRF close to 1. 0 and ESS above 200 to assume convergence. All trees that predated the time needed to reach a log-likelihood plateau were discarded as burn-in, and the remaining samples were used to generate a 50% majority rule consensus tree (option halfcompat). Branch support was estimated and PP PART ≥ 0. 95 was considered to indicate strong support for a given clade  (Douady et al. , 2003). All MrBayes analyses were performed on the computer cluster CIPRES Science Gateway  (Miller et al. , 2015) , using BEAGLE  (Ayres et al. , 2012)  with default parameters. BI was also conducted using PhyloBayes MPI 1. 8  (Lartillot et al. , 2013)  under the CAT-GTR-Γ 4 mixture model  (Lartillot and Philippe, 2004) , which has proven to perform well on large molecular datasets (e.
g. Allio et al. , 2020). The analyses were conducted on the nucleotide dataset. For each analysis, two independent MCMC starting from a random tree were run at least 10,000 cycles, with trees and associated model parameters being sampled every 10 cycles. The initial 2,000 trees sampled in each MCMC run were discarded as burn-in after checking for convergence in both likelihood and model parameters (tracecomp subprogram), and clade posterior probability (bpcomp subprogram). We checked the mean (meandiff) discrepancy observed across all bipartitions, considering that meandiff < 0. 01 indicates convergence. The 50% majority-rule Bayesian consensus tree and the associated posterior probabilities (PP CAT ) were then computed from the remaining trees using bpcomp.
Estimation of divergence times We first tested the hypothesis of a molecular clock with PATHd8  (Britton et al. , 2007). Since a strict molecular clock was not supported for 72. 7% of the nodes in this dataset at P < 0. 05, a Bayesian relaxedclock approach considering rate variations across lineages was employed to estimate divergence times  (Drummond et al. , 2006). Bayesian MCMC analyses implemented in BEAST 1. 10. 4 (Suchard et al. , 2018)  were performed to approximate the posterior distribution of rates and divergence times and infer their credibility intervals. Molecular dating analyses can be sensitive to several parameters or priors  (Alfaro and Holder, 2006; Drummond and Bouckaert, 2015). For instance the prior governing the branching process  (Condamine et al. , 2015)  or the number of molecular clocks  (Angelis et al.
, 2018; Condamine et al. , 2018a)  may impact estimates of divergence times. Accordingly, we conducted analyses with the following non-default settings and priors: (1) the molecular dataset was partitioned according to the best-fit scheme recovered by PartitionFinder, (2) substitution models were set based on ModelFinder estimation in IQ-TREE (see above), (3) a birth-death tree prior was set for the branching process, and (4) one clock model was set for the mitochondrial partitions, and one clock for the nuclear partitions. Finally, we let the tree topology be estimated by BEAST (classic operator mix) but enforced the monophyly of subfamilies, tribes, and genus Papilio following the results of phylogenomic analyses  (Allio et al. , 2020). Bayesian relaxed clocks were set up with an uncorrelated lognormal distribution clock model, with the mean set to a uniform prior between 0 and 0.
Hypothesis tests on various systematic positions Marginal likelihood (SS) Bayes factor Unconstrained MCMC each and with a sampling frequency of 5,000 generations. We discarded the first 10% of generations as burn-in and checked for convergence using Tracer (ESS > 200). We combined the four runs using LogCombiner. Using TreeAnnotator, we computed the maximum clade credibility tree with median ages and the 95% credibility intervals (CI) at each node. The xml files are available in FigShare (Data 2: https://doi. org/10. 6084/m9. figshare. 13318124. v1).
Inference of historical biogeography We estimated the ancestral ranges of origin and geographic range evolution for Papilio using the ML approach of dispersal-extinctioncladogenesis (DEC,  Ree and Smith, 2008)  as implemented in the DEC eXtended version (DECX, Beeravolu and Condamine, 2016; available at: https://github. com/champost/DECX). To infer the biogeographic history of a clade, DEC requires a time-calibrated tree, the current distribution of each species for a set of geographic areas, and a time-stratified geographic model that is represented by connectivity and dispersal scalar matrices for specified time intervals spanning the entire evolutionary history of the group. DECX allows classical vicariance as a cladogenetic event by using temporally flexible constraints on the connectivity between any two given areas following the movement of landmasses and dispersal opportunity over time.
Investigating heterogeneity of diversification rates To provide an assessment of diversification rates through time, we used the ML approach of  Morlon et al. (2011) , tested and automated in  Mazet et al. (2023)  and implemented in the R-package RPANDA 2. 0  (Morlon et al. , 2016). This method aims at capturing heterogeneity of diversification by allowing preselected subclades to follow different birth-death models from the deeper pruned tree (the backbone). For both subclades and backbones (different backbones are tested because of the different combination of shifts), speciation and/or extinction rates can change exponentially through time: λ(t) = λ 0 × e αt with λ 0 denoting speciation at present, α the trend of rate variation of speciation through time t. Extinction rate can exceed speciation, meaning that diversification rates can be negative  (Morlon et al.
These internal branches were the source of discrepancy between the Bayesian analyses (under both the mixture model and the partitioned analyses) and ML analyses. The main discrepancy was that Clades 3a and 3b were recovered as sisters in ML analyses and weakly supported (BP PART = 23, and UFBS PART = 58; Fig. S1-2 ) but not recovered in Bayesian analyses (Fig. 3 , Fig. S3 ). Despite a substantially increased taxon sampling compared to previous studies (37% of total species diversity in  Condamine et al. , 2012  versus 78% in this study), these results suggest that the systematic backbone of Papilio is still not definitely resolved, which impedes a higher-level systematic revision of the genus. We also acknowledge the limited size of the Sanger-based dataset and Clade 3a was always composed of the subgenus Euchenor (Joshi and Kunte, 2022) including two New Guinean species (Papilio euchenor and P.
Acknowledgements This project has received funding from the  European Research Council (ERC)  under the European Union's  Horizon 2020 research and innovation programme  (project  GAIA , agreement no. 851188 ), and an ' Investissement d'Avenir ' grant from the  Agence Nationale de la Recherche  (project  CASMA, CEBA, ref. ANR-10-LABX-25-01 ) to F. L. C. ; a  Ramanujan Fellowship  from the  Dept. of Science and Technology, Govt. of India , to K. K. , and support of the  Dept. of Atomic Energy, Govt. of India , under project nos. 12-RandD-TFR-5. 04-0800and 12-RandD-TFR-5. 04-0900to   TIFR/NCBS ; a  Natural Sciences and Engineering Research Council of Canada (NSERC)   Discovery Grant  to  F. A. H. S. Sequencing  of most  Indian Papilio  was done at  NCBS  under research and voucher specimen collection permits issued to K.
Observations and model ensembles The observations used here for the attribution are the E-OBS v23e dataset of daily minimum temperatures  (Cornes et al. , 2018). The above indices are calculated for each grid point, spaced every 0. 25 • in this dataset, and then averaged over the study area. For the attribution of the frost event, we use five model ensembles. Each simulation of each ensemble was biasadjusted using the cumulative distribution function transform (CDFt) method  (Vrac et al. , 2016) , using the daily minimum and the daily average temperatures from E-OBS over the 1981-2020 period. Bias correction is an important step here, since GDD calculations use a threshold. This method was assessed for use in climate services in  Bartók et al. (2019)  and showed good performance. We used statistics of pooled ensembles, using data until 2021 for the GEV fit of the distributions.
0 simulations. To keep the ensemble balanced we retained a maximum of three realizations per model. Not all CMIP6 models could be processed at the time of the study. Models are detailed in the Supplement and constitute an ensemble of 45 simulations. The fourth ensemble used is a set of 10 SST-forced High-ResMIP simulations  (Haarsma et al. , 2016). For the historical time period , the SST and sea ice forcings used are based on the observed dataset, and for the future time period (2015-2050) the SST and sea ice are derived from CMIP5 RCP8. 5 simulations and a scenario as close to RCP8. 5 as possible within CMIP6. The analysis of this ensemble was carried out using the observed GMST as for the observations. The fifth ensemble is the same set of models run in coupled mode, and the model GSATs were used.
We conclude that at a local scale, variability is dominating trend signals (Table  2 ). Results here differ from the rapid attribution analysis  (Vautard et al. , 2021)  in the completion and adjustment of the E-OBS dataset by the producers. This led to slightly different values for the observed indices in 2021. For instance, the estimation of the TNnGDD250 index-based return period was estimated here to be 8 years instead of 12 years in the rapid attribution. However, the results are qualitatively similar to those found in the preliminary analysis.
• C d -1 , with an average value on 5 April of 259 • C d -1 over the study domain. This value is high for this calendar day (rank = 14th since 1921 in the E-OBS extended dataset), but the record value was obtained in 2020, with a mean GDD of 320 • C d -1. Given the range of values taken in the domain, we considered 3 thresholds for GDD: 250 • C d -1 as a central value and 150 and 350 • C d -1.