dataesr/scientific-paragraphs-categorization huggingface dataset:
S31-48 ). Similarly, Anaplasma, Neorickettsia, and Wolbachia (Anaplasmataceae) are among the longest branches in the Rickettsiales (see their Fig. S50 ; see also our Fig. S2 ). All these species are secondarily, and not ancestrally, less compositionally biased, i. e. , they evolved from species with A+T-rich genomes. Moreover, their analyses were based on a rather small dataset that comprised only 18 or 24 mitochondrion-encoded genes (5,583 and 6,643 sites, respectively) and fewer than 41 taxa. These factors may, in combination, have led to the inference of poorly supported trees (e. g. , see their Figs. S31-40), and an artefactual attraction between mitochondria, the Rickettsiales, and the FEMAG I and II groups (i. e. , Fast-Evolving MAG I and II; see their Fig. 4 ). Several previous studies have suggested that mitochondria were either sister to the Rickettsiales  [18] [19] [20]  or phylogenetically embedded in a larger group comprised of both the Rickettsiales and the Holosporaceae  20.
The completeness and contamination of the bins reported by MaxBin v. 2. 2. 2 were assessed with CheckM v. 1. 0. 12  57. Genome bins that were phylogenetically affiliated to the Alphaproteobacteria based on the manual examination of the CheckM reference genome tree (itself based on the concatenation of 43 marker genes) were retained. Reads were then individually mapped onto each alphaproteobacterial genome bin with Bowtie2. All paired and unpaired reads that successfully mapped to the alphaproteobacterial bins were subsequently co-assembled with metaSPAdes. The resulting co-assembly was processed through the Anvi'o metagenomic workflow  58. In brief, reads were mapped to the final metaSPAdes co-assembly with Bowtie2 to obtain contig coverage values. DIAMOND searches  59  of predicted proteins against the NCBI GenBank nr database were done to assign taxonomic affiliations to each contig.
Marker protein selection We built an expanded dataset of mitochondrion-and nucleus-encoded proteins of alphaproteobacterial origin in eukaryotes. For the nucleus-encoded proteins, BLAST  64  similarity searches of all proteins contained in the predicted proteomes of 13 representative eukaryotes were conducted against a database of 176 prokaryotes (136 bacteria and 40 archaea). BLAST hits were clustered into homologous families with a custom Perl script, aligned with MAFFT and the L-INS-I method  65  , and then trimmed with BMGE  66. Phylogenetic trees for each homologous gene family were inferred under the LG model in RAxML v. 8  67. These trees were then sorted based on the criterion that eukaryotes form a clade with alphaproteobacteria. Manual inspection of the trees then followed to remove paralogs and contaminants.
For mitochondrion-encoded genes, mitochondrial clusters of orthologous genes (MitoCOGs)  68  that are widespread among eukaryotes were used. Both mitochondrion-, and nucleus-encoded candidate marker proteins were then compared through BLAST searches against those reported previously by Wang and Wu (2015)  20  and Martijn et al. , (2018)  11. Our dataset encompassed most proteins from these other datasets, with few exceptions. The nonredundant and remaining candidate marker proteins comprising the union of these five datasets, were then further screened phylogenetically. Using a representative eukaryotic (mitochondrial) query for each marker gene, BLAST searches were done against a database that comprises 107 diverse bacteria (representing 27 cultured phyla) and 23 diverse eukaryotes (representing 6 major groups); eukaryotes were selected based on the availability of both mitochondrial and nuclear genomes or transcriptomes (see Table  S7 ).
Homologues were aligned with MAFFT, alignments trimmed with TrimAl 69 and singleprotein trees inferred with IQ-TREE  70. The single-protein trees were inspected visually to remove duplicates, paralogues, and any other visual outlier such as extremely divergent sequences. Singleprotein trees were then re-inferred from the curated alignments and visually inspected. Proteins for which trees showed a sister relationship between eukaryotes and alphaproteobacteria were kept for further analyses. Finally, these candidate marker proteins were annotated and further refined using the EggNOG database and BLASTp searches. The final marker proteins set comprised 108 genes, 64 of which are exclusively nucleus-encoded, 17 are exclusively mitochondrion-encoded, and 27 are both mitochondrionand nucleus-encoded (Fig.
Dataset assembly To increase taxon sampling as much as possible, MAGs reported in Anantharaman et al. , (2016)  33   39  were added to those reconstructed here (see Metagenomic analyses). To improve the quality of our MAG selection, MAGs were analyzed with the CheckM lineage workflow and those with quality values (completeness -5x contamination) lower than 50 were discarded, just as done before by Parks et al. , (2017, 2018)  39, 40. MAGs were then filtered according to their taxonomic affiliation to the Alphaproteobacteria. A phylogenetic tree for all MAGs and all Proteobacteria taxa in the GTDB r89 database  40  was inferred from 120 marker proteins, built-in in the GTDB-Tk software, using IQ-TREE v. 1. 6. 10 70 and the LG4X+F model. To increase phylogenetic accuracy, a second tree was inferred with the LG+PMSF(C60)+G4+F using the LG4X tree as guide.
All MAGs that fell within the Alphaproteobacteria clade in the GTDB-Tk tree were chosen for subsequent analyses. Together, these added up to more than 3,300 alphaproteobacteria. In order to reduce computational burden, Treemmer v. 0. 1b was then used to reduce the number of alphaproteobacterial taxa from the GTDB-TK tree while maximizing phylogenetic diversity  41. The Treemmer analysis was constrained so representatives from major clades, as visually identified, were retained. Finally, a set of reference alphaproteobacteria (formally described species) were added, and long-branching alphaproteobacteria were replaced by short-branching relatives. To retrieve homologues, PSI-BLAST searches with either one, two, or three iterations using representative mitochondrial (eukaryotic) query sequences for each marker protein were done against a database that comprised all carefully selected predicted proteomes.
To remove non-orthologous sequences, homologous protein sets were retrieved for each marker protein, aligned with MAFFT, trimmed with TrimAl and trees inferred with IQ-TREE. The single-protein trees were visually inspected to remove duplicates, paralogues, and any other visual outlier such as extremely divergent sequences. The curated homologous protein sets were finally aligned again with MAFFT v. 7. 3. 10 and the L-INS-I method. To increase phylogenetic signal by removing poorly aligned and non-homologous aligned regions, Divvier v. 1. 0 was used with the -partial and -mincol options  71. Only sites with more than 10% of data were retained. To reduce incongruency among proteins due to, for example, lateral gene transfer, Phylo-MCOA v. 1. 4  72 was employed on single-protein trees with UFBoot2+NNI as branch support which were inferred with IQ-TREE v.
Phylogenetic analyses using site-heterogeneous models For multi-protein phylogenetic analyses on the supermatrix, trees were first inferred in IQ-TREE v. 1. 6. 10 under the LG4X+F model. The resulting site-homogenous tree was then used as a guide tree to infer a new phylogenetic tree under the LG+PMSF(C60)+F+G4 model  75. Consequently, the resulting siteheterogenous tree was used as a guide tree to infer a new phylogenetic tree under the dataset-specific LG+PMSF(MAM60)+F+G4 model. The dataset-specific MAM60 model was estimated using the MAMMaL software  43. This site-heterogeneous mixture model is directly inferred from the dataset analyzed and therefore is more specific than the general C10-60 mixture models. To account for more than 60 (e. g. , C60 or MAM60) amino-acid composition profiles across the data, we used the general UDM128 mixture model as LG+UDM128+G4+F that allows for 128 amino acid composition profiles  76.
The software FunDi was used to estimate functionally divergent sites in the branch that separates the mitochondrial lineage from all other taxa  42. Sites with a probability of being functionally divergent > 0. 5 were removed. Progressive removal of compositionally heterogeneous sites was performed according to the ɀ and the χ 2 metrics/methods as described before  11, 29, 44. Both metrics are designed to estimate compositional heterogeneity per site based on different criteria. Bayesian analyses were conducted with PhyloBayes MPI v1. 8 using the CAT-GTR+G4 model  77, 78. PhyloBayes MCMC chains were run for >20,000 cycles or until convergence between the chains was achieved and the largest discrepancy in posterior probabilities for splits between chains ('max-diff') was <0. 1. Individual chains were summarized into a Bayesian consensus tree using a burn-in of 500 trees and subsampling every 10 trees.
Using these estimates, the values in the class frequency vectors,   (𝑐)  , for any site profile class are modified in the following way to be branch-e-specific class frequencies,   (𝑐𝑒). The modified class frequencies have to satisfy a number of constraints including: This leads to non-linear equations for μ  (ce)  ,    (𝑒)  and   () that are solved numerically for each branch e to generate the modified class frequencies. For each branch and site class c, πj  (ce)  values are used to create a new transition Q (ce) matrix for likelihood calculations for all site patterns over that branch. The same approach is used with frequencies coming all extant taxa to obtain the root frequencies. A software implementation of GFmix is available at https://www. mathstat. dal. ca/~tsusko/software. html. 𝜋 𝑗 (𝑐𝑒) = { 𝜇 ( Partitioning the data matrix for GFmix calculations.
The foregoing framework assumes that for each aligned protein in a given concatenated dataset, the GARP/FIMNKY ratios (be's) for every branch in the tree will be similar. However, for our data matrix this where Nu,v is the total number of taxa for which sequences were available for both proteins (this normalization accounts for the differing amounts of missing data for different proteins). The proteins were then clustered based on  , distances using the UPGMA algorithm in MEGA-X 79 and 10 clusters were chosen as a computationally tractable number of partitions for further analysis. The GFmix model was then applied to these 10 partitions allowing for separate be values and branch lengths for each partition. The overall log-likelihoods for topologies were obtained as the sum of log-likelihoods of that topology over all partitions.
Profile Hidden Markov Model (pHMM) searches To search for bacteriochlorophyll enzymes, a set of 17 custom-made pHMMs for the genes bchB, bchC, bchD, bchE, bchF, bchG, bchH, bchI, bchJ, bchL, bchM, bchN, bchO, bchP, bchX, bchY, bchZ was used against predicted proteomes from the MAGs reconstructed in this study. These pHMMs were created from manually curated sets of bch genes from diverse proteobacteria. The searches were done with the program hmmsearch of the HMMER suite using an E-value cut-off of 1E-25. To search for mitofilindomain containing mic60 genes, the Pfam pHMM for Mitofilin (PF09731) was used with its own GA cut-off value.
Acknowledgements SAM-G is supported by a  EMBO Postdoctoral Fellowship  ( ALTF 21-2020 ). We are thankful to  Bruce Curtis  ( Dalhousie University ) and  Dayana Salas-Leiva  ( Dalhousie University ) for assistance with scripts, and to  Wendy Valencia  ( Harvard University ) and  Camilo Calderon  ( Rutgers University ) for advice on Python and R. This work was supported by the Moore-Simons Project on the  Origin of the Eukaryotic Cell ,  Simons Foundation  grants  735923LPI  (DOI: https://doi. org/10. 46714/735923LPI )awarded to AJR and  GBMF9739  (DOI: https://doi. org/10. 37807/GBMF9739 )awarded to PLG, and  Discovery Grants  from the  Natural Sciences and Engineering Research Council of Canada  awarded to AJR, ES, and  CHS.
Data Availability Sequencing data were deposited in NCBI GenBank under the BioProjects PRJNA315555, PRJNA438773, PRJNAXXXXXX, PRJNAXXXXXX, PRJNAXXXXXX, and PRJNA703749. Assembled metagenomes, novel alphaproteobacterial MAGs, and gene files (unaligned, aligned, and aligned and trimmed) are available at: DOI: 10. 6084/m9. figshare. 14355845. Datasets and phylogenetic trees inferred in this study are available at: DOI: http://dx. doi. org/10. 17632/dnbdzmjjkp. 1. The GFmix model software is available at: https://www. mathstat. dal. ca/~tsusko/software. html.
assumption is not true as different proteins show different degrees of GARP/FIMNKY variation across taxa depending on the location of the corresponding gene (e. g. , nucleus-encoded vs. mitochondrialencoded) and degree of conservation. For this reason, we clustered the proteins in our dataset into groups in the following way. For each protein v and each taxon t we calculated the GARP/FIMINKY ratio, 𝑏 𝑣 (𝑡) = 𝜋 𝐺 (𝑡) 𝜋 𝐹 (𝑡) ⁄ and v in the data matrix as 𝑑 𝑢,𝑣 = ∑ |𝑏 𝑣. Then, we calculated the overall distance between these ratios for every pair of proteins u (𝑡) -𝑏 𝑢 (𝑡) | 𝑡 𝑁 𝑢,𝑣 ⁄.
Cancer cells in similar functional states are found in all glioblastoma, despite the genomic heterogeneity observed between and within these brain tumors. Metabolism being downstream of all signaling pathways regulating cell behaviors, we looked for metabolic weaknesses in link with motility, a key functional state for glioblastoma aggressiveness. A signature--driven data reduction approach highlighted motile cells present in thirty tumors from four independent single--cell transcriptomic datasets. Analyses integrating trajectory modeling disclosed, as characteristic of motile cells, enhanced oxidative stress coupled with mobilization of the cysteine metabolism enzyme 3--Mercaptopyruvate sulfurtransferase (MPST). The soundness of this prediction was verified using migration and invasion assays with patient--derived cells and tissue organoids.
Three major migration/invasion routes have been described: the narrow and tortuous extracellular space of the brain parenchyma  (Scherer 1938) , the perivascular spaces surrounding blood vessels  (Scherer 1938)   (Montana & Sontheimer 2011 )  (Watkins et al 2014)   (Cuddapah et al 2014)  and white matter tracts, notably the corpus callosum, used as a highway for bilateral brain hemisphere invasion  (Scherer 1938)   (Pedersen et al 1995). Both single cell and collective migration have been reported (Vollmann--Zwerenz et al 2020)  (Volovetz et al 2020) , the latter appearing to be facilitated by the formation of an interconnected network of tumor microtubes  (Osswald et al 2015). So far, an integrated view of the metabolic pathways at play in a motile GB cell is lacking. Here, we applied a signature--driven data reduction approach to highlight subpopulations of cells in a motile state present in thirty tumors from four independent single--cell transcriptomic datasets.
Grouping cells according to their potential motile state We used four independent publicly available single--cell transcriptome datasets generated from surgical resections of patients bearing IDH--wild--type GB. Two were obtained with the SMART--seq2 technology  (Neftel et al 2019)   (Darmanis et al 2017)  and two with the 10X Genomics technology  (Neftel et al 2019)   (Pombo Antunes et al 2021). They are hereafter designated as N--S, D--S, N--10X and PA--10X, according to the initials of the first author of the paper reporting its first description and the technology used. The computational analyses were first performed using the N--S dataset, and the robustness of the results evaluated by repeating the analyses with each of the other three datasets. This choice stems from the deeper sequencing depth of SMART--seq2 data compared to 10X Genomics data, and the higher number of malignant cells and distinct patient tumors in the N--S dataset compared to the D--S dataset.
Enrichment of motility--related genes in M HIGH cells We evaluated the coherence of our findings by confronting them with previous experimentally--acquired knowledge. In gene ontology (GO) enrichment analyses, we found an enrichment of motility--related terms among the genes overexpressed at least two--fold in M HIGH cells (log2 fold--change≥2), providing a first support to the relevance of our approach (Fig. 1g , Supplementary Table  3 --sheet2). To further probe this relevance, we applied the same analytical strategy to the other three selected datasets. THBS1 being detected in <1% malignant cells in the PA--10X dataset, it was excluded from the motility signature for that dataset. Each of these analyses resulted in clusters composed of cells coming from different tumors (NMI =0. 058--0.
15) (Supplementary Figure  1C--E ). Differential expression analysis between M HIGH and M LOW cell clusters (Mann Whitney, BH--adjusted p--value < 0. 01) (Supplementary Table  2 --Sheet3--5), followed by GO enrichment analyses of the most overexpressed genes, again highlighted enrichment in motility--related terms (Supplementary Figure  2 , Supplementary Table  3 --sheet3--5). To determine the similarities of M HIGH cell transcriptional profiles from distinct datasets, we compared the lists of their overexpressed genes. The highest overlaps were found between lists of genes derived from datasets obtained with the same sequencing technology. 69. 2% (4852/7010) of the genes overexpressed in M HIGH cells from N--S dataset were also overexpressed in M HIGH cells from D--S dataset, and 86. 0% (6627/7707) of the genes overexpressed in M HIGH cells from N--10X dataset were also overexpressed in M HIGH cells from PA--10X dataset (Supplementary Figure  1F ).
Comparing the four lists of genes overexpressed in M HIGH GB cells, we found that 43 to 53% of the genes overexpressed in M HIGH GB cells from one dataset were also overexpressed in M HIGH cells from the three other datasets (Supplementary Figure  1F , Supplementary Table  2 --Sheet6). These results indicate that our analytical strategy captures cells with similar transcriptional profiles in distinct datasets. To further determine the extent to which our analytical approach captures cells in a motile state, we selected three gene modules previously found to be associated with GB cell motility (Supplementary Table  2 --Sheet7), and determined their enrichment in M HIGH cells and their correlation with the motility score. Most motile cells, including GB cells have been shown to undergo transcriptional reprogramming akin to the so--called epithelial--mesenchymal transformation (EMT) process.
We therefore used an EMT gene module of 14 genes associated with mesenchymal transformation of neural cells and GB cell motility  (Carro et al 2010 )  (Armento et al 2017). We also sought for enrichments in two other gene modules mobilized in invasive GB cells. One consists of 36 genes highlighted by Bhaduri and colleagues with highly correlated expression in an outer radial glia (oRG)--like malignant cell population with increased invasive behavior in GB  (Bhaduri et al 2020). The other consists of 32 genes involved in the regulation of GB cell motility and driven by TEAD transcription factors (Tome--Garcia et al 2018). Each of these three gene modules was found to be enriched in M HIGH cells in all analyzed dataset (Fig. 1h ). We also observed a positive linear relationship between the motility score and the EMT, oRG and TEAD scores (Fig.
Metabolic characteristics of cells with high motile potential We went on to identify metabolic pathways deregulated in GB cells with high motile potential. The genes coding for enzymatic components of the metabolic pathways overexpressed in M HIGH cells were identified in each dataset using the 2019 KEGG list of 1586 metabolic elements  (Kanehisa et al 2017) (Supplementary Table 4--Sheet2). KEGG pathway enrichment analyses highlighted metabolic pathways related to extracellular matrix (ECM) modeling (e. g. glycosaminoglycan synthesis/degradation) and membrane composition rearrangements (e. g. unsaturated fatty acids (FA), FA elongation, GPI-anchor biosynthesis) (Fig. 2a , Supplementary Table  4 --Sheet3--6). In addition to these metabolic pathways expected to be active in motile cells, we also observed enrichment in pathways involved in energy production (Fig.
Overall, the results of this analysis point to an enhanced mobilization of anti--oxidative pathways counteracting the deleterious actions of oxidative compounds derived from the over--mobilization of energy production pathways. Trajectory modeling was used next to delineate the acquisition of metabolic components potentially crucial for a motile cell state. Trajectory reconstruction based on expression of the motility signature genes was achieved with the STREAM python software using the N--S dataset which has the highest sequencing depth. The inferred trajectory corresponded to a branched path, highlighting an intermediary branch (S3--S1 branch) linking cells with low motile potential to cells with high motile potential (Fig. 2b ). Among the metabolic enzymes overexpressed in S3--S1 branch compared to each of the other branches, characterizing therefore the switch between low and high motile states (Supplementary table  5 ), two were overexpressed in M HIGH cells of all datasets considered.
These results suggest that MPST could be necessary for GB cells moving along blood vessels as well as through the brain parenchyma. In the different datasets considered in this study, the edge or core origin of the cells is most often unknown, though they are most likely coming from dense tumor cell regions. Removal of the surrounding parenchyma is usually undertaken with great care by neurosurgeons in order to spare at the best motor and cognitive functions. Although Darmanis and colleagues distinguished malignant cells isolated from the tumor core and the periphery, the cancer cells isolated from the latter are too few to authorize robust computational analyses. It is therefore very likely that our findings relate to dense tumor areas' cells, moving either along the blood vessel network that characterize these highly angiogenic tumors, or across the tumor parenchyma enriched in components favoring cell displacement.
scRNA--seq data pre--processing Four publicly available single--cell transcriptomic datasets from GB patient tumors were downloaded, including the only two publicly--available datasets obtained with the SMART--seq2 technique  (Neftel et   6 ). The last two datasets were selected among possible others because of the availability of the raw UMI counts. We designated them as N--S, D--S, N--10X and PA--10X, respectively. Outlier cells with low--complexity transcriptomes, revealed by plotting the graphical distribution of the number of detected transcripts and genes per cell, were excluded from the analyses as described (Saurty--Seerunghen, 2019). For 10X datasets, detection of over 20% mitochondrial genes per cell was considered as indicative of dying cells leading to their exclusion from the analyses.
Therefore, we retained for further analyses 4916 malignant cells from twenty adult GB patients in the N--S dataset, 1033 malignant cells from four adult GB patients in the D--S dataset, 5797 malignant cells from six adult GB patients in the N--10X dataset and 8666 malignant cells from six adult GB patients in the PA--10X dataset. Log 2 -transformed Counts Per Million (log 2 (CPM+1)) were used for D--S and the 10X datasets, unless otherwise specified. CPM corresponds to the counts of gene--mapped reads normalized by the total number of mapped reads per cell divided by one million, thus allowing comparison of read abundance across libraries of different sizes. For the N--S dataset, gene expression data currently available are in log 2 ((TPM/10) + 1). TPM (Transcripts Per Million) normalization takes into account an eventual biased estimation of long transcript numbers, by dividing the number of mapped reads by the transcript's length.