dataesr/scientific-paragraphs-categorization huggingface dataset:
Staining intensity was reported by calculating the optical density of staining using MATLAB R2020a. Image analysis for Ki67 and TGM-1 antibodies were calculated using 2 fields of view from 2 biological samples with a total 2 donors per antibody. For Filaggrin, 2 fields of view per sample, 3 biological samples per donor, 5 donors with N = 30 in total were used for calculations. All the image analysis for filaggrin samples used the same parameters. Multiplex cytokine assay and ELISA. In this study, supernatant (media secreted from the explants) from the Day 1 and Day 6 timepoint was collected to evaluate 20 cytokines and chemokines. Media from the Day 6 timepoint was accumulated from the most recent media change on Day 4 of culture. Multiplex cytokine array (Inflammation 20-Plex Human ProcartaPlex™ Paneld, Thermofisher, Waltham, MA, USA) was performed using the Luminex™ MAGPIX™ Instrument System (Thermofisher, Waltham, MA, USA) according to manufacturer instructions.
qPCR reactions were run using validated Taqman gene expression assays. The plates were run using the Life Technologies QuantStudio 12K Flex instrument. Each gene was assayed in duplicate. Primers for K5, K14, NOTCH1, PCNA, and PPARD were used according to Applied Biosystems protocol. Results were normal-ized using reference gene PPIA, obtained from Applied Biosystems. Final results were interpreted as normalized RQ values to the control treatment group. Statistical analysis. Data was assessed using Student's t-test, and two-way ANOVA with multiple comparison Holm-Sidak test. The differences between groups were compared using Graphpad Prism 9. 0. 1 (GraphPad Software Inc. , La Jolla, CA, USA).
Motor evaluation Spontaneous locomotion activity was evaluated using the open-field test. Freely moving mice were placed in the center of a white polyvinyl chloride open-field apparatus (50 cm long × 50 cm wide × 45 cm high) for 15 min. The video tracking software (Viewpoint, Lyon, France) allowed us to determine virtual areas: a central square (30 cm long) and a peripheral zone (corridor of 10 cm wide). Total distance travelled as an index of locomotion and time spent in the central area as an index of anxiety-related behavior were recorded using the video tracking system. Motor coordination was determined using a standard mouse rotarod (TSE, Homburg, Germany) using the standard operating procedures described by the EUMORPHIA group program  (24). Briefly, a training phase identifying mice able to stay on the rod at 4 rpm for 60 s was first realized.
Cognitive evaluation Working memory was evaluated using the Y-maze according to the standard operating procedure described by the EUMORPHIA group program  (24). The maze was made of three equal arms (40 cm long × 10 cm wide × 16 cm high) in black polyvinyl chloride radiating at 120°from each other, providing the shape of Y. Mice were placed at the end of one arm of the maze, facing the wall, and allowed to explore freely for a 10-min session. Latency to leave the first arm, total number, and sequence of entries into each arm were scored for each mouse. An arm entry was counted when the mouse had all four paws inside the arm. If a mouse completed fewer than eight arm entries during the entire session, its data were excluded from further analysis. Spontaneous alternation was defined as entries into all three arms on three consecutive choices.
Statistical analysis All statistical analyses were carried out using GraphPad Prism 6 software (GraphPad Inc. , San Diego, CA, USA). Equality of variance and normal distribution were evaluated before each analysis to determine the type of analysis to be performed: parametric or non-parametric. For body weight and EAE score longitudinal monitoring, data were analyzed using a two-way analysis of variance (ANOVA) with groups (CTL + vehicle, EAE + Vehicle, CTL + GSE, and EAE + GSE) and day post-induction (DPI, D0-D30) as main factors and with DPI defined as a repetitive measure. Posthoc comparisons were made with Tukey's comparison test for multiple comparisons between groups for each DPI. For motor, sensitivity, and cognitive behavioral tests, data were analyzed using two-way ANOVA with groups (CTL + vehicle, EAE + Vehicle, CTL + GSE, and EAE + GSE) and day post-induction (DPI, D0-D30) as the main factors, with DPI being defined as a repetitive measure.
Finally, anxiety-related behaviors were evaluated by 1) the time to enter the first new arm using the Y-maze test, 2) the time spent in the central zone using the open-field test, and 3) the time spent in the open arms using the EPM. In all three tests, no significant difference between groups was observed (Tukey's test, p > 0. 05) (Figures  2H-J ). Taken together, these data suggest that GSE as well as our model of EAE mice had no effect on motor coordination, locomotion activity, working memory, and anxiety-related behaviors. in EAE + GSE (Tukey's test, p > 0. 05) (Figures  3A, E ). At D30, during the chronic phase of EAE, there was no statistical difference in MDA between groups (Tukey's test, p > 0. 05) (Figures  3C, G ). In the brain, an increase in the content of carbonyl protein was observed in EAE mice (46% at D20 Tukey's test, p < 0.
Finally, this study did not evaluate the production of reactive species nor the aberrant redox signaling leading to oxidized macromolecules due to the choice of our therapeutic protocol. Treatment with GSE was initiated at D10, when the EAE development was already well ongoing (first symptoms are already observed). We then considered the characterization of the first step of oxidative stress development in EAE mice beyond the scope of this article. However, the evaluation in the near future of the effect of GSE (using earlier treatment initiation) on the production of reactive species and aberrant redox signaling would be of huge interest. As a hallmark of MS pathology, the demyelination process was assessed by Western blotting quantification of MBP and CNPase protein expression. Our data fully support CNS demyelination during the acute and chronic stages of EAE, in support of literature  (41, 42).
1% formic acid) and solvent B (HPLC-MS grade acetonitrile/0. 08% formic acid). For untargeted metabolites, the gradient run time was 17. 5 min: 5 min at 2% B followed by an increase to 10% B for 7 min, then an increase to 80% B for 0. 5 min, and to 100% for 1 min, then down to 2% B for 4 min. LC-MS and data-dependant LC-MS/MS data were acquired in positive ion mode in the m/z = 100-1300 range. Analysis quality control was ensured by injecting a mixture of all samples every 10 runs to monitor separation and check for cross-contamination. For SCFA, the gradient run time was 7 min: 3 min at 20% B followed by a linear increase to 50% B for 0. 5 min, increase to 100% B for 1 min and down to 20% B for 2. 5 min. LC-MS data were acquired in negative ion mode in the m/z = 100-250 range. The retention times and m/z ratios of the SCFA 3-NPH derivatives and corresponding stable deuterated internal standards are provided (See Additional file 1: Supplementary Table  4 ).
Differences were regarded significant at P < 0. 05. Genomic data were visualized and analyzed using the ampvis2 R-package in RStudio and the web-based tool MicrobiomeAnalyst. The Shannon diversity index was analyzed to compare further the treatment effects. Beta diversity was evaluated by principal coordinate analysis (PCoA) based on the Bray-Curtis distance  [25]. To obtain species-level classification, sequences of dominating OTUs were analyzed in the EzBioCloud database  [26]. PERMANOVA analysis, using a permutation method under a reduced model, was used to study the significant differences (P < 0. 05) between the different groups. Correlation network analysis was based on the SparCC algorithm, with the permutation set at 100, P value at 0. 05 and a correlation threshold of 0. 3 at the genus level.
Acknowledgements The authors warmly thank  Mr. Antoine St-Cyr  and  Sollio Agriculture  for the help with the successful completion of this study. Liya Zhang  acknowledges the support of the  China Scholarship Council  through a Ph. D. fund scholarship. We thank the  Plateau technique de spectrométrie de masse bioorganique  ( Muséum National d'Histoire Naturelle, Paris, France ) for acquiring the mass spectra.
Faecalibacterium prausnitzii is abundant in the healthy human intestinal microbiota, and the absence or scarcity of this bacterium has been linked with inflammatory diseases and metabolic disorders. F. prausnitzii thus shows promise as a next-generation probiotic for use in restoring the balance of the gut microbial flora and, due to its strong anti-inflammatory properties, for the treatment of certain pathological conditions. However, very little information is available about gene function and regulation in this species. Here, we utilized a systems biology approach-weighted gene co-expression network analysis (WGCNA)-to analyze gene expression in three publicly available RNAseq datasets from F. prausnitzii strain A2-165, all obtained in different laboratory conditions. The co-expression network was then subdivided into 24 co-expression gene modules.
To this end, expression data can be used to accurately group genes into functional modules based on co-expression patterns. In particular, the algorithm developed for weighted gene co-expression network analysis (WGCNA) enables the construction of gene networks through consideration of the coexpression patterns between two genes as well as the overlap between neighboring genes  [27]. Highly correlated genes are clustered into larger modules based on similarities in their expression profiles, and members of a given module are often involved in similar functional processes  [28]. To date, WGCNA has been successfully used to construct gene co-expression networks in several bacteria, including Mycobacterium tuberculosis  [29] , Escherichia coli  [30] , Lactococcus lactis  [31] , Vibrio cholerae  [32] , and Streptococcus oralis  [33].
In the current study, gene expression data were obtained for F. prausnitzii A2-165 from three public RNAseq datasets. Co-expression networks were analyzed using WGCNA and genes were clustered into modules with similar expression patterns. For each module, we investigated enrichment in biological functions. The structure of the network was also analyzed to identify ovh/). The shinyapp is now freely accessible to the public via the URL https://shiny. migale. inrae. fr/app/ faeprau. hub and bottleneck genes; of the 25 we identified, the majority were located in the region of six mobile genetic elements that could promote intracellular or intercellular DNA mobility. This study represents the first time that a gene expression network has been constructed for F. prausnitzii. In addition, an R-Shiny application is available online to make the data and the analyses accessible and usable to a wide audience (https://shiny.
Data sources In the Lebas dataset  [22] , the authors explored the effects on the transcriptomic profile of F. prausnitzii A2-165 of treatment with supernatant from Lactococcus lactis subsp. lactis CNCM I-1631, Lactococcus lactis subsp. cremoris CNCM I-3558, Lactobacillus paracasei CNCM I-3689, and Streptococcus thermophilus CNCM I-3862. The raw RNA-seq fastq files are available at ArrayExpress under project E-MTAB-9387. In the D'hoe dataset, the authors sequenced RNA from F. prausnitzii A2-165 monocultures at three different time points  [21]. The RNAseq results were deposited in the Short Read Archive under the study identifier SRP136465. In the Kang dataset  [23] , the authors compared the transcriptomic profiles of F. prausnitzii A2-165 cultured using different carbohydrates (galactose, fructose, glucose, N-acetylglucosamine) as the sole carbon source.
Data processing The workflow for this study is presented in S1 Fig. A read count table was obtained using the raw sequencing reads. After trimming reads with fastp v. 0. 20. 0 (default parameters,  [34] ), fastq-formatted reads were aligned to the genome of F. prausnitzii strain A2-165 (Genome Assembly ASM273414v1) using BWA v. 0. 7. 17  [35] , allowing a single mismatch in the read. Then, sam-formatted alignments were sorted and converted to bam output files using SAMtools v. 1. 10 [36]. The number of reads per transcript from each sample was counted using HTSeqCount v. 0. 12. 4  [37]  and GFF-formatted gene annotations downloaded from NCBI. We checked the distribution of raw counts and performed principal component analysis in each dataset (S2 Fig). Gene expression values were normalized using the DEseq2 package v.
1. 34. 0, in R  [38]. The count table, containing 2,950 genes (S1 Table ), was filtered to eliminate nonexpressed genes. The resulting final dataset contained 2,902 genes and was further processed with WGCNA  [39]. The quality of the expression matrix was evaluated by hierarchical clustering based on the distance between different samples, measured using Spearman's correlation. No outliers were detected (S3 Fig).
Construction of weighted gene co-expression network WGCNA is widely used in systems biology to construct a scale-free network from gene expression data  [41]. In this study, we used the WGCNA package in R (v. 1. 71-3). In brief, Pearson's correlation matrices were calculated for all pairs of genes and were subsequently transformed into matrices of connection strengths using a power function; in this case, the power was β = 17. We selected pairs of genes with distance correlation coefficients greater than 0. 5 from the three datasets. The connection strengths were then used to calculate topology overlap (TO)  [42] , which measures the connectivity of a pair of genes. Hierarchical average linkage clustering  [43]  based on TO was used to identify gene co-expression modules, which groups genes with similar patterns of expression.
Functional analysis and networks After extracting the list of genes present in each module, we used eggNOG-mapper, a tool for functional annotation based on precomputed orthology assignments  [44]. Functional enrichment analysis was performed in each module with STRING (version 9. 0; https://string-db. org/ ) in March 2022. Enrichments with strength > 0. 5 were selected. The resulting interaction networks contained direct (physical) and indirect (functional) interactions, derived from numerous sources including experimental repositories and computational prediction methods. Nodes with a high degree of connectivity are called "hubs"; they interact with several other genes and may thus play a central role in the network. Instead, "betweenness" measures the total number of shortest nonredundant paths passing through a given node or edge.
Comprehensive transcriptomic profiles from RNAseq data Raw data from previous gene expression profiling analyse  [21] [22] [23]  were preprocessed using a similar workflow (S1 Fig). Details on these experiments and the samples used are presented in Table  1. To investigate differences in the transcriptomic profiles among the different conditions analyzed, we conducted a PLS-DA. This revealed a clear separation along the gene expression datasets from the different publications along the first two principal components, which explained 40% and 9% of the observed covariance (Fig  1A ). While component 1 was sufficient to distinguish the Lebas dataset, the second component was necessary for separation of the D'hoe and Kang datasets. In total, we identified 50 genes that could be used to discriminate between the Lebas and Kang samples (component 1) and 50 genes that could be used to distinguish the D'hoe and Kang datasets (component 2) (Fig  1B ).
Exploration of transcriptomic data with R-Shiny We developed an R-Shiny application to facilitate deeper exploration of the transcriptomic data from each of the three studies. For example, using this approach on the Lebas dataset, we observed that treatment of F. prausnitzii with supernatant from L. lactis subsp. lactis generated more transcriptional downregulation than upregulation relative to the control (see our R-Shiny application), confirming the findings of the original study  [22]. In the Kang dataset, a comparison between bacteria raised on glucose or fructose as the sole carbon source highlighted that the most upregulated gene was CG447_08360, which encodes a putative PTS glucose transporter subunit IIBC (Fig 2A  and 2B ). This gene was also highly upregulated in the presence of galactose.
As the transcription PTS genes is usually regulated in response to substrate availability  [48] , this suggests that CG447_08360 encodes a galactose/glucose transporter. Finally, the D'hoe dataset provided the opportunity to track F. prausnitzii gene expression across three growth time-points. In the R-Shiny application, it was possible to compare, for example, expression at time 15h with that at time 3h and to highlight the most upregulated and downregulated genes between these time-points (Fig 2C  -2H ). Notably, serine/ threonine and ferrous iron transporters were highly upregulated at the earlier time compared to later, which suggests that available serine, threonine, and ferrous iron are imported into the cell during early growth (Fig 2D  and 2E ). It was also interesting to observe the upregulation of CG447_14545 at time 3h (Fig 2F ); this gene encodes a multi-antimicrobial extrusion (MATE) efflux transporter, which plays a role in antibiotic resistance as well as potentially in bacterialhost interactions and intercellular signaling  [49, 50].
Instead, time-point 15h was marked by strong upregulation of CG447_09920 and CG447_06950, which encode energy-coupling factor (ECF) transport S component and DNA polymerase IV, respectively (Fig 2G  and 2H ). ECF transporters are responsible for vitamin uptake in prokaryotes  [51] , while DNA polymerase IV is part of the SOS response to DNA damage in bacteria and contributes significantly to cell fitness in late stationary phase cultures in the absence of any exogenous DNA damage  [52]. Our R-Shiny application thus provided new insight into each dataset and reveals genes that may play a role in the fitness of F. prausnitzii under the different growth conditions employed.
Construction of weighted gene co-expression network To choose the best soft-thresholding power, we analyzed network topology and determined that 17 was the lowest appropriate power value, with a scale-free topology fit index of 0. 95 and a relatively high average connectivity (Fig  3A ). WGCNA revealed 24 distinct modules representing highly co-expressed networks of genes (Fig 3B  -3D ). The list of genes in each module of the WGCNA analysis is available in S3 Table. Two modules, turquoise and blue, were more notably larger than the rest, with 568 and 431 genes, respectively (Table  2 ). The grey module contained the non-clustering genes, of which there were only 36, suggesting that WGCNA performed well under the parameters applied here. Individual modules can represent independent units that are responsible for certain biological functions.
By linking the gene expression of the modules with the experimental conditions employed, it may be possible to identify modules that have important functions under certain conditions. To this end, we used WGCNA to correlate each module with the conditions used in the three publications  (   24 samples were then visualized (Fig  4B ). Several large modules represent expression changes that clearly delineate between Lebas on one side and D'hoe and Kang on the other. Indeed, the blue, cyan, pink, brown, and purple modules (Pearson correlation coefficient � 0. 75) correspond to genes that were upregulated in the D'hoe and Kang datasets, while the green, salmon, red, and yellow modules (Pearson correlation coefficient � -0. 75) represent genes that were upregulated in the Lebas dataset. This suggests that the expression patterns of the genes in these modules reflect the ways in which the growth of F.
Functional and pathway enrichment analysis The top biological processes associated with the large brown and blue modules were "structural constituent of ribosome" and "regulation of translation", respectively (Table  3 ). Several modules were enriched in metabolic functions such as arginine, histidine, cobalamin, or fatty acid biosynthesis (green, red, lightgreen, or royal blue modules, respectively; Table  3  and Fig   5). We also identified clusters of genes related to bacteriophage function, chaperones, stress response, or SOS response, suggesting the involvement of functions related to adaptive responses (darkred, tan and black modules). Overall, we found that genes involved in similar pathways or with the same biological function tended to belong to the same expression cluster with WGCNA, which supports the relevance and utility of this integrative bioinformatics approach in this context.
In their original study, Lebas et al. showed that F. prausnitzii A2-165 responds to cell-free supernatants from lactic acid bacteria by downregulating mobilome genes and upregulating cell-wall related genes. Here, we observed that the purple module (downregulated in Lebas relative to D'hoe and Kang) contained clusters of genes related to transposon and mobilization proteins, while conversely, the yellow module (upregulated in Lebas relative to D'hoe and Kang) was linked with processes associated to the cell wall and the regulation of cell shape. In addition, we found that some genes involved in mobilization or cell-wall synthesis form potential networks of interaction (Table  3 ). In this way, the results from our WGCNA strengthen the observations by Lebas et al.
Data Availability Statement: Not applicable. Acknowledgments: Imen Ghzaiel (I. G. ) received financial support from Nutrition Méditerranéenne et Santé (NMS; N. L. ), Institut Européen des Antioxydants (IEA; S. M. ) and was awarded the NMS prize (Romanée-Conti) in 2021. I. G. also received financial support from ABASIM (Association Bourguignonne pour les Applications des Sciences de l'Information en Médecine; Dijon, France). Mohamed Ksila (M. K. ) received financial support from PHC Utique (G. L. /T. G. and O. M. -K. ; 2021-2022; code CMCU22G089/code Campus France 47608V). A. G. A. , A. Z. , F. B. and G. L. are members of the International Natural Product Sciences Taskforce (INPST: https://inpst. net/ accessed on 20 April 2022).
In Silico Metabolite Prediction Sibiriline putative metabolites were predicted using Sibiriline SMILES string (OC(C=C1)=CC=C1C2NC3=NC=CC=C3C2) through the free online software GLORYx  [31]  and Biotransformer 3. 0  [32]. Briefly, GLORYx integrates machine learning-based site of metabolism prediction with reaction rule sets to predict and classify the putative structures of metabolites that could be formed by phase I and/or phase II metabolism. Biotransformer 3. 0 is an open-access software tool that supports the rapid, accurate, and comprehensive prediction of the metabolism of small molecules in both mammals and environmental microorganisms. All metabolite structures were generated using ChemDraw software (PerkinElmer, Inc. , Waltham, MA, USA).
MN generation. Spectral data allowed us to generate MN using a semi-quantitative approach. Data acquisition, processing (i. e. , MS data conversion, preprocessing, MS1 annotation, and generation of molecular networks), visualization, and network analysis have already been described  [19]. Briefly, raw data were converted to an open MS format (. mzXML) with ProteoWizard's MSConvert module  [34]. The mzXML files were then preprocessed (deconvolution, de-isotoping, alignment, gap-filling) with MZmine 2 software  [35]. The single mgf output file was then loaded on the Global Natural Products Social networking (GNPS) web-based platform in order to generate the multi-matrix molecular network  [18]. For the use of high-resolution data, the basic parameters were modified to m/z 0. 02 for the mass tolerance of precursor and fragment ions used for MS/MS spectral library searching, and m/z 0.