dataesr/scientific-paragraphs-categorization huggingface dataset:
There were no overlaps between overrepresented Pfam domains, suggesting the colonization stage may be correlated with phylogenetic proximity and consequently with the species' Pfam contents. The same enrichment analysis using MCL counts (Table  S3 ) revealed 316 and 171 enriched clusters in early and late colonizers, respectively, and the majority of clusters showed no functional annotations (early -77%, late -89%). The clusters uniquely enriched for the early colonizers included three CAZymes clusters (CBM14, cutinase CE5, and SUN bglucosidase GH132) and other clusters with functional annotations such as STAND proteins (PF17111), methyltransferase (PF13489), ankyrins (PF12796), fungal lectin (PF07938), heterokaryon incompatibility (HET) (PF06985), CFEM (PF05730), arthropod defensin (PF01097), BTB/POZ (PF00651), F-box (PF00646), cysteine-rich secretory protein (PF00188), transcription factors (PF00172, PF00096), proteases (PF00082, PF00026), protein kinase (PF00069), p450 (PF00067), and heat-shock proteins (PF04119, PF00012).
Interestingly, clusters containing hsp70 (PF00012), STAND proteins (PF17111), and ankyrins (PF12796) were enriched in both early and late colonizers (Table  S2 ). In late colonizers, we found enriched CAZymes, an LPMO (AA9) and CMB1. The nodes of the phylogenetic tree in Fig. 1 (a) show the number of gene families expanded and contracted (family-wide P ≤ 0. 01) using the program CAFE  (De Bie et al. , 2006). Among the pyrophilous fungi, only P. domesticum and G. carbonaria showed higher ratios of expansions : contractions (4. 19 and 1. 22, respectively). For the rapidly evolving families (Viterbi P ≤ 0. 01) in the pyrophilous genomes (Table  S4 ), only P. omphalodes showed a ratio of expansions : contractions < 1, totaling 114 expanded and 56 contracted families. Among the expanded gene families, we found CAZymes and other enzymes potentially involved in the degradation of PyOM, such as LPMO (AA9 -Pezech1), gluco/chitooligosaccharide oxidase (AA7_dist -Pezech1), endoglucanase (GH6-CMB1 -Pezech1), GMC oxidoreductase (AA3_2 -Pyrdom1, Pyrom1), cutinase (CE5 -Tri-pra1), cytochrome p450 (Tripra1, Pyrdom1), endo-b-1,4glucanase (GH7 -Wilmi1), and xylanase (GH10-CBM1 -Wilmi1).
We performed the same analysis using Pfam domains and found 78 overrepresented and 98 underrepresented entries (Fig. 4c ). In this case, we can see a clear separation between pyrophilous Pezizales and Agaricales fungi. Pfam domains prevalent in Pezizales genomes show functions involved in defense against microbial attack and stress, such as chitin binding peritrophin-A (PF01607 or CBM14), arthropod defensin (PF01097), tocopherol cyclase (PF14249), ricin-type beta-trefoil lectin (PF00652), heat shock protein 9/12 (PF04119), fungal fucosespecific lectin (PF07938) and fungalysin metallopeptidase (PF02128). The phylo-PCA (Fig. 4d ) does not show a separation of pyrophilous fungi but still has a strong signal separating plant pathogens and ectomycorrhizal fungi. Saprotrophy, the main pyrophilous underlying lifestyle, is more functionally diverse and scattered over the plot.
Charcoal modulates developmental genes in P. domesticum Pyronema spp. are known to form abundant fruiting bodies on burned soil a few weeks after large forest fires  (Bruns et al. , 2020) , and P. omphalodes (syn. P. confluens) is an established model organism for the analysis of cell biology and fruiting body development in filamentous fungi  (Traeger et al. , 2013). Therefore, we hypothesized that exposure to the products of pyrolysis could trigger sexual development in these fungi. To check this, we have performed a gene co-expression analysis of the two available datasets (see the Materials and Methods section). A study by  Fischer et al. (2021)  focused on the metabolic pathways used by P. domesticum to utilize charcoal as its sole carbon source. Since the two strains are very closely related, we found direct orthologs for 73.
However,  Fischer et al. (2021)  showed that Pyronema species, which are the most dominant early post-fire colonist, also have the ability to metabolize more recalcitrant forms of highly aromatic carbon. Our future investigations will address whether other pyrophilous fungi share this ability. The gene content and enrichment of gene families and functional domains suggest an inferior competitive ability of early pyrophilous colonizers. Specifically, we found that they have lower numbers of secondary metabolite clusters and other parasitic signatures such as chitinases and proteases. The only potentially defense-related expanded families found mainly in the Pyronematacae family are the CBM14 (peritrophin-A -PF01607) and arthropod defensin (PF01097). The CBM14 domain is mainly present in Metazoa, especially in insects and fungi.
4 ). Since the heatmap color scale in Fig. 4  represents a z-score, it appears that some species have low counts, but in reality it means that they are on the lower side of the distribution. The development of mature fruiting bodies follows specific genetically encoded programs that determine the species-specific morphologies  (Krizs an et al. , 2019). To examine whether charcoal could induce sexual development in P. domesticum, the first responder found at the Rim Fire site in 2013 and in experimental pyrocosms  (Bruns et al. , 2020) , we compared transcriptome datasets from two studies where P. domesticum was grown on charcoal  (Fischer et al. , 2021)  and sexual reproduction was induced by white/blue light  (Traeger et al. , 2013). We found some modules containing genes induced in both the charcoal and sexual reproduction conditions through co-expression analysis, and they were enriched in GO terms like amino acid biosynthesis, response to oxidative stress, and defense response.
, 2015) and by eliminating inhibitory compounds present in the substrate  (De Groot et al. , 1998). Interestingly, cutinase (CE5) is an expanded family in pyrophilous fungi and is differentially expressed on charcoal. Since they are not plant pathogens, and this strain of P. domesticum can consume charcoal as a carbon source  (Fischer et al. , 2021) , this family might be involved in the degradation of hydrophobic products of combustion. In summary, we performed a comparative genomics analysis of pyrophilous fungi and we found the enrichment of gene families involved in response to stress and degradation of PyOM. Also, we found that the protein length and GC3 content of pyrophilous fungi lie between those of nonpyrophilous and thermophile fungi. The transcriptome analysis found overlapping modules of genes co-expressed in charcoal and light-induced sexual development conditions in P.
(Miyauchi et al. , 2020)pyrophilous genomes. of complete BUSCO models using the Ascomycetes dataset. N50, count of smallest number of contigs whose length sum makes up half of genome size. L50, weighted median statistic such that 50% of the entire assembly is contained in scaffolds equal to or larger than this value. Wilcoxina mikolae was sequenced in a different study(Miyauchi et al. , 2020).
Acknowledgements We are grateful to  Prof. Irina Druzhinina  for kindly offering advice on fungal ecology and the structure of the manuscript. The work (project  10. 46936/10. 25585/60001080 )by the  US Department of Energy (DOE) Joint Genome Institute , a  DOE Office of Science User Facility , is supported by the  Office of Science of the US DOE  under Contract no. DE-AC02-05CH11231. The work was also funded by the  Department of Energy  grants  DE-SC0016365  to TDB and IVG and  DE-SC0020351  to TDB, MT and IVG.
Data availability Genome assemblies and annotations are available at MycoCosm (https://mycocosm. jgi. doe. gov) (Grigoriev et al. , 2014)  and have been deposited at DDBJ/ENA/GenBank under the following accessions: Pyronema domesticum, CBS 144463 (JAJJXL0000 00000); Pyronema omphalodes, CBS 144459 (JAJJXM00000 0000); Geopyxis carbonaria, CBS 144460 (JAJJXG000000000); Tricharina praecox, CBS 144465 (JAJJXF000000000); Morchella snyderi, CBS 144464 (JAJJXB000000000); Peziza echinospora, CBS 144458 (JAJLQT000000000).
Similarly,  Crabot et al. (2020)  found that the temporal variability of community structure was related to the frequency and duration of drying events. However, few studies have examined temporal β diversity of stream invertebrate communities over large spatial and temporal scales. Here, we complement recent findings from  Mouton et al. (2020) , who reported decreasing spatial taxonomic β diversity (taxonomic homogenisation) but increasing spatial functional β diversity (functional differentiation) among river macroinvertebrate assemblages. We capitalise on the same high-resolution time-series datasets, comprising macroinvertebrate communities collected annually from 1991 to 2016, at 64 mainstem river sites across New Zealand's two mainland islands. The two studies are complementary in the sense that  Mouton et al.
1 m 2 area in the sampler down to a depth of ca 10 cm and from as many substrate types as possible. Individuals were later identified in the laboratory, to the lowest practicable taxonomic level (species = 26%; genus = 47%, family = 21%; Quinn and Hickey 1990). The same taxonomic resolution was maintained throughout the entire period to allow analyses of long-term changes  (Scarsbrook et al. 2000). We described each macroinvertebrate taxon (n = 113) using sixteen functional traits related to morphology, lifehistory, dispersal strategies and resource acquisition methods extracted from the New Zealand freshwater macroinvertebrate trait database (NIWA, <https://niwa. co. nz/sites/niwa. co. nz/files/nz_trait_database_v19_2_18. xlsx>). This trait database has been explicitly developed for New Zealand's standardised freshwater macroinvertebrate sampling protocols  (Dolédec et al.
e. the coefficient of variation × 100; Fick and Hijmans 2017), mean air temperature and air temperature seasonality (i. e. the standard deviation of the mean). We applied a principal component analysis (PCA; the dudi. pca function in the Ade4 package ver. 1. 7-15 Dray and Siberchicot 2020) individually to our sets of trends in air temperature and precipitation variables. Based on the correlations among each trend in climate variable and the axes of the PCAs, we created synthetic indices of climate change using the first two axes of each climate PCA. In the PCA of trends in air temperature variables, the first axis (35. 9% of variation explained) described a gradient of increasing mean air temperature (TMean) and the second (33. 0%) described changes in air temperature seasonality (TSeas). For precipitation variables, the first axis (28.
2010). For reach-scale descriptors, we used: 1) the slope (degrees) of the stream segment at each sampling site (SegSlope;  Leathwick et al. 2010) , 2) the predicted wetted river width (m) at the 7-day mean annual low flow (WidthMALF; Booker 2015), 3) the estimated proportion of riparian shading at each river segment (SegRipShade; measured from satellite imagery by  Leathwick et al. 2010) , 4) temporal changes in water quality, 5) temporal changes in flow and 6) temporal changes in substrate size. Temporal changes in water-quality, flow and substrate size were all estimated as slopes of variables (described below) and years. For temporal changes in waterquality, we compiled data from 1991 to 2016 of median annual values (mg m -3 ) of nitrate (NO 3 -N), ammoniacalnitrogen (NH 4 -N), dissolved reactive phosphorus (DRP) and dissolved oxygen (DO) concentrations and water clarity (metres; CLAR).
These data originate from samples collected monthly, at the same sites as those sampled for macroinvertebrate communities  (Davies-Colley et al. 2011). We log 10transformed water-quality variables and calculated temporal changes for each variable at each site. We ordinated temporal changes in water-quality using PCA and retained the first axis (33. 3%) as a synthetic indicator of temporal changes in water quality (WQ). For changes in flow, we compiled mean daily flow data from hydrological gauges located at each site, which we divided by upstream catchment area to obtain a measure of run-off per day  (Vaughan and Gotelli 2019). We then calculated mean annual and seasonal (for winter, spring and summer) flow and flow coefficient of variation, as was done for precipitation. We ordinated temporal changes using PCA and kept the first two axes (Flow PC1, 54.
Temporal β diversity indices We used the temporal β diversity index from Legendre (2019) (TBI; the TBI function in the adespatial package ver. 0. 3-8) to measure temporal changes in the taxonomic and functional composition of each of the 64 communities. This index has been specifically developed 1) to calculate differences in assemblage composition (dissimilarity) between two time periods from pairwise distances, and 2) to decompose these temporal differences into indices of local colonisations and extirpations. Macroinvertebrate abundances were log(x + 1)-transformed and pairwise distances in community composition were calculated using the percentage difference index of dissimilarity (%diff; Odum 1950, also known as the Bray-Curtis index). To estimate temporal changes in functional composition (i.
This method allows to examine whether changes in composition of each assemblage (relative to the baseline year) consistently increase (positive slope) or decrease (negative slope) over time  (Antão et al. 2020). To test our first expectations (E 1a and E 1b ), we tested for relationships between trends in temporal β diversity, including colonisations and extirpations with latitude using linearmixed effects models (the lme function in the nlme package ver. 3. 1-152  (Pinheiro et al. 2017 )) setting island (North Island versus South Island) as random effect.
Drivers of temporal β diversity We performed hierarchical generalised additive mixed effects models (GAMMs; the gamm function in the mgcv package ver. 1. 8-33  (Wood and Wood 2015) ) to relate temporal changes in each of the six indices of trends in temporal β diversity to our set of environmental descriptors. Predictor variables were Box-Cox transformed prior to analyses  (Box and Cox 1964;  the BoxCoxTrans function in the caret package ver. 6. 0-84  (Kuhn et al. 2020 )), and standardised to zero mean, one unit variance (using the decostand function in the vegan package ver. 2. 5-6  (Oksanen et al. 2014) ). GAMMs were fitted with catchments nested within islands as random effects to account for spatial structure in data  (Dormann et al. 2007 ). We used regression splines to account for potential non-linear relationships, restricted maximum likelihood (REML) to optimize the parameter estimates and assumed Gaussian-type distribution errors.
Relationships with taxonomy and functional traits To test our third expectations (H 3a and H 3b ), we started by relating changes in community composition to the taxonomic classification of our taxa (H 3a ). For this purpose, we compiled a taxonomic classification of the studied taxa (using phylum, subphylum, class, subclass, family, subfamily, order, genus and species names), from the National Center for Biotechnology Information (using the taxonomy function in the myTAI package  Drost et al. 2018). From this classification, we created a taxonomic tree by calculating taxonomic distances between taxa (using Gower's distance and the taxa-2dist function in vegan). We used hierarchal clustering (the hclust function in R) to produce a tree from the taxonomic dissimilarity matrix. We checked for potential relationships between taxonomic distance and Euclidean dissimilarity in population size and/or species range shifts using a Mantel test  (Mantel 1967)  with 9999 permutations (the mantel.
Sample Preparation and NMR-and MS-Based Analysis NMR profiling was accomplished with a 500 MHz Bruker AVANCE III 1 H NMR spectrometer (Advance III, Bruker, Germany) on an automatic sample changer with cooling capability. The NMR spectra were acquired via two complementary experiments on each plasma sample: one-dimensional (1D) 1 H nuclear Overhauser effect spectroscopy NOESY1D presat (NOESY1dgppr sequence)  [25]  and Carr-Purcell-Meiboom-Gill (CPMG presat)  [26]. The spectral width was set at 10 kHz. A Fourier transformation was carried out on the NMR data with the line broadening (LB) set at 0. 3 Hz. The spectra were phased and baseline corrected, then all spectra were aligned on a glucose doublet at 5. 23 ppm. Finally, each 1D NMR spectrum was sliced into buckets of 0. 001 ppm, containing NMR signals.
These steps were accomplished using NMRPipe software  [27]. Before the statistical analysis, the NMR spectra bins were normalized to the median of intensities  [28]. The assignment of the 1 H NMR spectra was performed using two-dimensional (2D) experiment, Chenomx software and the Human Metabolome Database (HMDB) (version 4. 0)  [29]. Raw MS data were collected, converted to mzML format files using ProteoWizard software  [30]  and processed in each ionization mode (ESI+ or ESI-) as a single batch using the XCMS pipeline  [31] , which was implemented on MetaboAnalyst 5. 0  [32]. This generated MS feature tables. The features were filtered to remove those with more than 50% non-detected responses in samples before further analysis. The missing values were replaced with the average of the corresponding non-missing values in the k (here, k = 5) closest features in terms of the Euclidean distance of the responses across all samples.
Data Analysis The processed 1D 1 H NOESY NMR data were used for multivariate analysis. An orthogonal partial least squares discriminant analysis (OPLS-DA) was carried out using an in-house MATLAB OPLS script, which was based on the Trygg and Wold method  [33]. The statistical analysis and visualization were performed on MATLAB ® R2016b for macOS (Mathworks, Natick, MA, USA). The quality parameters of the models, the explained variance (R2Y) and the predictability of the model (Q2Y) were calculated. Q2Y was computed via 7-fold cross-validation and confirmed by exploring the impacts of permutations in the dataset rows. Variables missing more than 30% of values in the clinical data were removed from analysis and the missing values were imputed using the predictive mean matching method with the R package mice  [34].
Box and whisker plots were generated using the R packages ggplot2 (https://ggplot2. tidyverse. org/, accessed on 1 October 2022) and ggprism  [35]. A metabolome-wide association study (MWAS) was accomplished using a linear regression framework with the R package xmsPANDA (https://github. com/kuppal2 /xmsPANDA, accessed on 1 October 2022). For each feature, the log 2 transformed response was used to test the relationships across the different stages of MAFLD, including HC, ST and NASH, as continuous variables (HC = 0, ST = 1 and NASH = 2). The model included adjustments for age (continuous) and gender (categorical), which are known to influence plasma metabolomes  [36, 37]. To account for the multiple comparisons, Benjamini and Hochberg-based false discovery rate (FDR) thresholds of 0. 2 and 0.
Identities were assigned using the evidence scoring for matching to adducts at a ±5 ppm mass tolerance for metabolomic data and a ±10 ppm for lipidomics data, with the scoring threshold set to 0. 05. Metabolic activity patterns associated with MAFLD were characterized using Mummichog. A multiblock analysis was carried out using the R package mixOmics  [39]. The metabolomic and lipidomic features and clinical characteristics selected by MWAS were included in the multiblock model to identify highly correlated multi-omics signatures that discriminated the MAFLD stages: HC, ST and NASH. Additional details are provided in the Supplementary Materials.
Study Design and the Characteristics of HC, ST and NASH Patients A simplified diagram of the study design is shown in Figure  1A. The NMR data and the LC-MS-based metabolomic and lipidomic data were collected from the baseline plasma of 85 (87. 6%), 80 (82. 5%) and 82 (84. 5%) patients, respectively, due to the limited number of plasma samples. Among the 82 patients (lipidomics), there were 66 females (80%) and 16 males (20%) (Table  S1 ). The diagnoses of HC (Figure  1C ), ST (Figure  1D ) and NASH (Figure  1E ) were established histologically using liver biopsy specimens (Figure  1C-E ). Figure  1B  shows the positive correlation between plasma and liver TGs, demonstrating that elevated triglyceride levels, in addition to liver TG content, could be a key feature of MAFLD in obese patients.
64, an R2Ycum of 0. 48 (R2 is the indicator of how the model fits the data; the closer to 1, the better the fit) and a Q2cum of 0. 39 (Q2 is the capacity of the model to correctly class a new dataset; the closer to 1, the better the model predictability). The CV-AUC was 0. 896, which meant the probability that the OPLS-DA model scored a randomly chosen NASH patient who was classed as a NASH patient was 89. 6% higher than a randomly chosen NASH patient who was classed as an HC patient. All of which indicated that our OPLS-DA model for the classification of HC and NASH patients was fairly stable. The model variables that contributed most significantly to the classification of HC and NASH patients are represented in Figure  2C. An OPLS-DA covariance plot was generated, which restored the form of the NMR spectra and colored them according to their correlation scores with NASH).
Multiblock Integrative Analysis of MAFLD To identify highly correlated multi-omics markers that could discriminate between HC, ST and NASH patients, we performed a multiblock analysis  [39]  by integrating clinical, metabolomic and lipidomic data (Figure  S3A ). The circos plot in Figure  5A  shows the correlations between the variables of different blocks, represented on the side quadrants. It can be seen that the clinical variables were almost entirely negatively correlated with the other two DataFrames (metabolomic and lipidomic data), while the metabolomic and lipidomic variables were more mixed and were positively correlated with each other. The correlation cutoff was set at 0. 7 (cutoff = 0. 7). As shown in Figure  S3A , the first components from the metabolomic and lipidomic datasets and the lipidomic and clinical datasets were highly correlated with each other (as indicated by the large numbers in the bottom left of the figure).
Enriched Metabolic Pathways Associated with MAFLD To uncover the potential biological mechanisms underlying the MAFLD-associated MS features from our MWAS, pathway enrichment analyses were performed using Mummichog, which directly predicts metabolic activity from HRMS data without a priori MS feature identification. Mummichog has been shown to improve biological interpretation while balancing false positive and false negative errors in the enrichment of highdimensional mass spectral data  [45]. Our Mummichog analysis identified the significant enrichment of 16 pathways (Tables  1  and 2 ) associated with MAFLD: 7 were identified in both the metabolomic and lipidomic analyses, while 9 of them were only detected in one type of analysis (1 in the metabolomic analysis and 8 in the lipidomic analysis).
Acknowledgments: The authors would like to thank all participants of the study. We acknowledge the  LC-MS facility (CAS Key Laboratory of Separation Science for Analytical Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian, China). We would also like to thank  Yang Ouyang ,  Qiuhui Xuan ,  Zaifang Li ,  Qi Li ,  Xiaolin Wang ,  Nicolas Guedeney ,  Julia Deschamp  and  Agnès Victor Bala , for their technical advice. We gratefully acknowledge  Aiyana Washington  for her diligent proofreading of the manuscript.