dataesr/scientific-paragraphs-categorization huggingface dataset:
In order to find additional SNPs, 6,329 contigs (Supplementary Table  S2 ) similar to the sixteen complete small microviruses were retrieved from 2,946 assembled viromes and metagenomes/metatranscriptomes from the IMG database (tBLASTn of the contigs against all ORFs of the sixteen newly identified microviruses, e-value < 0. 001). To further evaluate selective pressure on ORFs, the pN/pS ratio was computed for each of them, dividing the number of non-synonymous mutations (pN) by the number of synonymous ones (pS). The identified contigs, mostly partial microvirus genomes, were then mapped to the sixteen genomes of interest, and SNPs were assessed using the methodology described previously. All measured parameters, namely LRS%, SNPs, pN/pS, and conservation, are summarized for all ORFs in Supplementary Table  S5.
To test whether the gene conservation, composition, and SNP analysis could help finding bona fide genes, we performed these analyses on the bullavirus phiX174 whose gene content is experimentally validated. For the gene conservation analysis, fifteen microvirus genomes ranging from 0 to 100 per cent intergenomic identity with phiX174 were randomly chosen from GenBank sequences and our in-house virome contigs. For these fifteen genomes and phiX174, ORFs were predicted and our two-step clustering strategy was applied on them. For the SNP analysis, twenty-five phiX-like genomes were selected from Gen-Bank (intergenomic identity > 95%; Supplementary Table  S4 ) and mapped to the phiX174 genome to detect SNPs. These SNPs were used as references to calculate the pN/pS ratio on the experimentally verified genes and the predicted ORFs of phiX174.
Finally, the percentage of high-degeneracy amino acids (LRS%) was computed for each ORF and gene of phiX174. For each of the sixteen genomes and phiX174, a genomic map was drawn and colored according to ORF conservation, LRS%, and pN/pS values using the Python library DNA Features Viewer  (Zulkower and Rosser 2020). SNPs were plotted with ggplot  (Wickham 2009)  as the sum of ratios of SNP depths over total depths across a sliding window of 15 nt, divided by 15 to normalize the window size. Additionally, the number of tBLASTn hits was evaluated as a heatmap using a 15-nt window, also plotted with DNA Features Viewer, in order to emphasize the conserved zones of each ORF. All these results were used to manually predict bona fide genes and refine their 5 ′ ends. Once all genes were defined, secondary structures of intergenic regions were predicted using RNAfold v2.
The sixteen newly detected genomes are significantly shorter than any known phages Using the HMM profiles of the Rep and Cap proteins from thirtythree reference microviruses, a total of 12,495 circular contigs and genomes similar to microvirus ones were identified in sequence data from GenBank, in IMG microbial metagenomes and metatranscriptomes  (Nowinski et al. 2019) , and in 2,946 assembled viral metagenomes. Although many contigs (1,314) were shorter than the smallest genome of a cultivated microvirus (4. 248 kb), only a handful of contigs stood out as extremely small (Supplementary Fig. S1 ). Eleven contigs smaller than 3. 5 kb were thus recovered, and their Cap proteins were compared to the other 12,484 sequences in order to identify related ones (BLASTp > 50 on bitscore). Five additional contigs, themselves small (<3.
692 kb), had a significant similarity with one of the former eleven smallest ones, leading to a total of sixteen sequences. These sixteen circular contigs (Supplementary Table  S3 ) had a region similar to at least one of the two HMMs, all being similar to the Rep profile (probability > 99%) and eleven similar to the Cap profile (probability > 99%). Two genomes were found in the GenBank database (assembled from viromes), three in microbiomes, one in a metatranscriptome, and ten in the set of public viromes we assembled. These sixteen genomes originated from fifteen different data sets derived from marine water (n = 7), sediments (n = 1), freshwater river (n = 1), wastewater systems (n = 4), and samples associated with cellular aquatic hosts, i. e. marine sponge (n = 1) and freshwater fishes (minnow, n = 2).
These five proteins were compared and were indeed similar to each other, although sometimes distant with a minimum BLASTp bit-score of 124. For the search of Cap proteins to be more sensitive, an HMM profile was built out of the multiple alignment of these five proteins and was compared to the HMM profile built from the alignment of the thirty-three reference Cap proteins and that of the eleven Cap proteins detected earlier. Although only distantly related, these two HMM profiles were significantly similar (e-value of 10 -5 ). In addition, protein structure prediction confirmed the presence of an SJR domain in these candidate Cap proteins. Thus, the two core proteins of microviruses were detected in all sixteen newly detected genomes. As a final verification, these sixteen Cap and Rep proteins were compared to all viral proteins from the KEGG (Kyoto Encyclopedia of Genes and Genomes;  Kanehisa and Goto, 2000)  database (BLASTp, e-value <10 -5 ) and all proteins identified as similar were proteins from microviruses.
e. SNPs) observed for each ORF was also used to further help discriminate between coding and non-coding ORFs. Different software exists to perform similar tasks such as Phy-loCSF  (Lin, Jungreis, and Kellis 2011)  or RNAcode  (Washietl et al. 2011) , and these methods are based on the analysis of a multiple alignment of at least ten similar genomes. Unfortunately, the sixteen genomes analyzed here do not have similar genomes at the nucleotidic level in any database, and the VIRIDIC analysis showed that these sixteen genomes are not even similar to each other at the nucleotidic level. Yet, high levels of microdiversity were shown for some viral species in a given environment and even in a given metagenome  (Martinez-Hernandez et al. 2017; Gregory et al. 2019 ). Thus, metagenome reads and contigs, even those representing only partial genomes, carry some microdiversity that could be used for detecting SNPs and then distinguish coding from non-coding ORFs.
As phiX174's genome was not assembled from metagenomic data, its SNPs were identified using twenty-five genomes from GenBank that are very similar to the phiX174 genome (at least 95% genome identity; Supplementary Table  S4 ). All non-coding ORFs had a number of negative SNPs higher than the number of synonymous ones. This result is expected, as negative SNPs are more frequent randomly than synonymous ones (Supplementary Fig. S8A ). Indeed, 168 negative and 138 synonymous SNPs are found when considering all codons and their nine possible SNPs. Conversely, the ten validated genes had more synonymous than negative SNPs. This result reflects a purifying selection on these coding regions at the codon level, with the selective removal of amino acid substitutions that are deleterious. Considering positive SNPs that are the most frequent randomly (n = 174), most genes here again favored synonymous mutations over them.
(E) For all ORFS, SNPs were detected among GenBank genomes very similar to phiX174 and indicated in blue if they lead to synonymous mutations, and green, orange, and red if they lead to positive, negative, and strongly negative mutations, respectively, according to the BLOSUM62 matrix. A bar chart displaying the abundance of each type of SNPs is plotted, with grey backgrounds delimiting the predicted ORFs. (1) on the left, the distance between the studied genome and all other similar genome (one line per genome) using VIRIDIC intergenomic distance and amino acid identity percentage for Cap and Rep; and (2) on the right, ORFs are displayed as columns in the order of the genome. (C) For each gene, the LRS% is computed and displayed. (D) Considering the SNPs detected (see below), ORFs are colored according to pN/pS values.
Rep was missing in one genome and was replaced by four genes in the reverse strand. Finally, the OG63 overprinted in Rep in reekeekeeviruses was not predicted by any in silico automated software. Overall, PHANOTATE predicted small genes in all intergenic regions longer than 100 nt, leading to a total of eighty genes instead of the sixty-nine genes. The five overprinted genes were not detected by automated software and were only recovered using our custom approach.
Tentative annotation of each conserved gene The custom gene prediction strategy resulted in the identification of seven proteins conserved in all or part of the sixteen genomes and three singletons found in the smallest genome reeSedAu. Among these proteins, only Rep and Cap were similar to any Pfam domains and to proteins of known function encoded by members of the Microviridae family. For all these OGs, especially the ones for which the function remained unknown, their secondary structure was tentatively predicted using ColabFold  (Mirdita et al. 2022 ) and  then compared to the known structures of phiX174 proteins (Protein Data Bank ids: 2BPA, 1CD3, 4JPP;  McKenna et al. 1992; Dokland et al. 1999; Berman et al. 2000; Sun et al. 2014). As previously described, the Cap proteins of the five reekeekeeviruses were so distantly related to the one of the known microviruses that they were difficult to identify even using remote homology tools.
DNA secondary structures In phiX174, the primosome, i. e. the protein complex containing a primase and helicase from the bacterial host, initiator of the complementary negative-strand synthesis, recognizes a hairpin region ahead of the spike coding gene G (Shlomai and Kornberg 1980). Figure  6. Protein structure prediction for the final inferred genes. PhiX174 is on top with structures retrieved from the Protein Data Bank (PDB): Cap, spike, and DNA-binding proteins were retrieved from the PDB entry 2BPA, while scaffolding proteins are solved in 1CD3, and the pilot protein is described in entry 4JPP. When not available, a prediction was run with ColabFold and marked with a star on the figure. Then, one representative genome was chosen for reekeekee-and roodoodooviruses, alongside the smallest assembled genome, to predict the structure of each of the proteins from the final curated genomes.
Quantification of these new clades in the environment The abundance of these two new clades was evaluated in the fifteen data sets in which the sixteen genomes were assembled, as well as in thirty-five additional viromes described in the same publications. For each metagenome, reads were compared using BLASTx to an in-house database of 13,390 Microviridae genomes that included the sixteen genomes. Unsurprisingly, the proportion of microviruses in viromes (Supplementary Fig. S13 right ) was found to be highly variable (between 0. 01 and 49. 91%). Among them, Gokushovirinae and Alpavirinae were almost always the most abundant groups (Supplementary Fig. S13 , left), followed by Pichovirinae and Group D. Even among the viromes from which they were assembled and despite the use of multiple displacement amplification in seven out of fifteen data sets, the reekeekee-and roodoodooviruses represent a small fraction of all microviruses (only up to 2 × 10 -4 and 13%, respectively).
Following library prep, all libraries were pooled in equal volume by combining 2 μL of each library for a final bead clean up with 0. 85X AmpPure beads (Beckman Coulter Life Sciences, Indianapolis, IN, USA). This QC pool was then sequenced on an MiniSeq MO flow cell (Illumina). The resulting index distribution was used to re-pool the libraries for an SP flow cell (Illumina) sequencing run with sample LR1 pooled at maximum volume available. Sequencing data generated in this study were deposited in NCBI Sequence Read Archive. Data from Movile Cave is available under bioproject accession number PRJNA673084. The genome of Thiovulum sp. from Frasassi is available under bioproject accession number PRJNA846597. The genomes of Ca. T. imperiosus, and Ca. T. karukerense are available under bioproject accession number PRJNA830902.
Thiovulum sp. genome assembly from public databases To obtain genomic information from additional cave-dwelling Thiovulaceae we downloaded all available metagenomic libraries from the Frasassi caves in Italy (SRR10997432, SRR1559028, SRR1559230, SRR1559353, SRR1560064, SRR1560266, SRR1560848, SRR1560849, SRR1560850, SRR8191123, SRR8194889, SRR8197024, SRR8200784, SRR8202337, SRR8203764). The short read libraries were quality-trimmed using Trimmomatic  [38]  and scanned for the presence of Thiovulum 16 S rRNA gene using PhyloFlash  [28] , revealing that library SRR1560850 contained >170,000 reads of Thiovulum sp. 16 S rRNA gene. A metagenomic assembly of library SRR1560850 was therefore conducted using Megahit  [39] , after which the assembly was binned using Metabat2  [40]. The obtained bins were taxonomically annotated using the GTDB-TK tool  [41]  resulting in one Thiovulum bin.
Genome quality assessment, annotation, and comparison The completeness of the Thiovulum genomes obtained was assessed using CheckM  [42]. The continuity of the Movile Thiovulum genome was evaluated using the unicycler-check module in Unicycler v0. 4. 9b. Annotation of all genomes was performed using the command-line Prokka  [43]  and DRAM  [44]  tools as well as the KEGG  [45] , EggNOG 5. 0  [46] , PATRIC  [47, 48] , and RAST  [49, 50]  annotation servers. A COG  [51]  analysis was done using the ANVIO tool  [52]. OperonMapper  [53]  was used to inspect the organization of genes into operons in the circular Movile Cave genome. CRISPRs where identified using CRISPR finder tool  [54]. Metabolic models of the annotated genomes were calculated using PathwayTools (V25. 5) [55]  using the RAST annotation of each genome as the model basis.
S rRNA gene phylogenetic trees A maximum-likelihood phylogenetic tree was calculated, including only 16 S rRNA gene sequences longer than 1200 nt, using FastTree 2  [60]  using all Thiovulum sequences in the SILVA REF (V138. 1) database (n = 71), three 16 S rRNA gene sequences obtained from the assembled genome of the Movile Cave Thiovulum (see below) and sequences of several different Sulfurimonas species as an outgroup. The tree was calculated using the generalized timereversible model and the discrete gamma model with 20 rate categories. Local support values were calculated with the Shimodaira-Hasegawa test.
Transcriptomic analysis The six libraries containing cDNA sequences (three from Air Bell 2 and three from Lake Room) were quality-trimmed using trimommatic  [38]  and mapped against the complete genome of the Movile Cave Thiovulum sp. using Salmon (version 1. 6)  [61]. Ribosomal RNA data were removed from the mapping results and TPM (Transcripts Per Kilobase Million) values were recalculated to reflect mRNA expression. Further removal of tRNA sequences did not alter the results. The RNA data was analyzed using the iDEP (v. 0. 95) online tool  [62]  that provides an online graphical user interface for the DeSEQ2  [63]  and Limma  [64]  packages for RNAseq analysis. Differential expression was considered significant with a 2-fold difference and a false discovery rate smaller than 0. 1. Taxonomic composition of the active community was obtained by analyzing the 16 S rRNA gene from the transcriptomic read libraries using PhyloFlash as above (V 3.
Samples were sputter coated with 10 nm gold and examined for SEM on a JEOL JSM 5510 LV microscope (Jeol, Japan), and for BFSTEM on a Hitachi SU8230 (Tokyo, Japan). Energy-dispersive X-ray spectroscopy (EDX) analysis was performed using an EDX analyzer (Oxford Instruments, Abingdon, UK) and with the INCA 300 software.
Genome analysis Of the four genomes assembled in this study, only the assembly of metagenomic data from Movile Cave resulted in a closed circular genomic sequence (Fig. S2 ). The size of the assembled genomes ranged between 1. 72-2. 20 MBP and GC content varied between 28-37% (Table  2 ). Genome completeness was estimated using CheckM  [42]  ranging between 90-96% and evalG (as incorporated into the PATRIC server (now BV-BRC);  [72] ) ranging between 95. 8%   and 100% with the exception of Thiovulum sp. Tui Malila which was evaluated at 66% completeness. Using a Campylobacteraceae specific set of marker genes did not improve the completeness prediction. The marker genes that could not be identified are listed in Supplementary Table  2. Of these 9 genes/domains (PF00308. 13, PF02464. 12, PF04085.
9, PF08299. 6, PF08459. 6, PF12344. 3, TIGR00615, TIGR02124, TIGR03423) could not be found in any of the 6 genomes and additional 4 genes/domains were missing from both freshwater genomes (TIGR00196, TIGR00539, TIGR00121, PF03853. 10). All genomes were analyzed using different tools with the results summarized in Supplementary Datasets 1-5. Coarse and Fine consistency as evaluated by EvalG  [72]  sulfur cycling, including two versions of the sulfide:quinone oxidoredutase (Type VI and IV) that oxidizes sulfide to polysulfide, and the polysulfide reductase gene nrfD that carries out the reverse process. In the Movile strain, nrfD was found in a 3-gene potential operon together with the large subunit of the assimilatory nitrate reductase (narB) and the ttrB; tetrathionate reductase subunit B.
Rhodanese sulfur transferase proteins were identified in all genomes. In the Movile strain these were located in a 6-gene operon containing two other subunits of a nitrate reductase (narH and narG). As for the two other genes in this operon, one is related to cytochrome c and the other could not be annotated (G2Y-810). The tauE sulfite exporter was identified in all genomes as well, in the Movile strain genome as part of a 5-gene operon involved in the transport of molybdate (modCABD). Sulfite dehydrogenase, dissimilatory sulfite reductase (dsrAB), the sox genes or adenylyl sulfate reductase (aprAB) that carry out the sulfide oxidation to SO 4 2-could not be found by any of the annotation tools nor by manual BLAST against all sequences available for each of those proteins in the UniProt database.
Differential expression analysis (Fig. 7C ; Supplementary Dataset 1) revealed that 222 genes were more expressed in the Lake Room compared to Air Bell 2, while the opposite comparison resulted in 42 genes. In contrast, the TPM normalized expression of the housekeeping genes rpoB, gyrB, and bipA was similar in both compartments, therefore, confirming that the differential expression of genes is not a methodological bias. Over half of the genes more expressed in the Lake Room encoded for hypothetical proteins to which no function could be assigned. Retron-type reverse transcriptases were the most dominant group of genes (n = 15) also exhibiting some of the highest transcription level with TPM values up to 19,000. Genes over expressed in samples from Air Bell 2 were related to energy generation including cytochromes c and b as well as F-type ATP synthase.
DISCUSSION We investigated the morphological, phylogenetic, and genomic aspects of a planktonic Thiovulum species that belongs to a newly recognized clade of freshwater Thiovulum. We were able to identify, and assemble, another genome of a freshwater Thiovulum from public metagenomic data from the Frasassi caves. We further compared these two genomes from freshwater clade to the genomes of four marine Thiovulum spp. : T. sp. ES  [11]  and T. sp. GCA_012963895. 1 (Tui Malila)  [71] , retrieved from public databases, and the genomes of Ca. T. imperiosus  [2]  and Ca. T. karukerense  [26]  whose genomes are presented here for the first time. The main aspects of this comparison are depicted in Fig. 6  and in more details in Supplementary Dataset 7. Supplementary Figs. 3 4 5 6 7 8 offer a detailed, Pathway-Tools automatic-reconstruction of the pathways detected in each of the six genomes.
Nevertheless, while genomic information is suggestive of the presence or absence of specific enzymes and pathways, additional experiments and gene expression data are required to determine which of the genes are utilized and under which environmental conditions. Our transcriptomic analysis points out that at the time of sampling the nap and εhao genes were minimally expressed as compared to other sulfur and nitrogen metabolism genes (Fig. 6 , Supplementary Dataset 1), suggesting that the DNRA pathway was not very active in the Thiovulum community at the time of sampling. In contrast, the high expression of both copies of rhodanase genes as well as the sulfite exporter (tauE) suggest that elemental sulfur may have been converted to sulfite and excreted. We propose that, if O 2 is available, sulfide is oxidized to elemental sulfur with oxygen as an electron acceptor (as may have been the case for part of the community at the time of sampling given the high expression of cytochrome c oxidase cbb3; Fig.
The presence and high activity of Thiovulum at the surface of both the Lake Room and Air Bell 2, environments that differ significantly in the overlaying atmosphere (normoxic vs. hypoxic, respectively), suggest that either the Ca. T. stygium strain Movile is not sensitive to these differences in O 2 concentration, or it mostly aggregates at different locations in the stagnant, oxic surface (<1 mm) of the different cave chambers. The expression profiles differed significantly between the Lake Room and Air Bell 2, and it is evident (Fig. 6 ) that most genes recognizable as involved in cell metabolism had higher expression levels in Air Bell 2, though not all at statistically significant levels (Supplementary Dataset 1). More than half of the genes overexpressed by Thiovulum in the Lake Room could not be assigned any function, making it impossible to understand it's specific metabolic activity in that compartment of the cave.
ACKNOWLEDGEMENTS The authors thank  GESS  team for logistics with sampling in the cave. We would also like to thank  Pheobe Laaguiby  at the  University of Vermont  for performing the outstanding  Oxford Nanopore  sequencing and  Bo Barker Jørgensen ,  Emily Fleming ,  Carl Wirsen ,  Tom Fenchel  and the anonymous reviewers for their valuable suggestions that improved the quality of this manuscript. We would also like to thank the  Extreme Microbiome Project (XMP)  for providing the DNA extraction reagents and methods as well as  Laura Gray  and  Mehdi Keddache at Illumina Corp  for providing partial sequencing reagents through its partnership with XMP. TB was supported by a grant of  Ministry of Research and Innovation , project number  PN-III-P4-ID-PCCF-2016-0016 (DARKFOOD ). SMS was supported by a grant of  Ministry of Research and Innovation (UEFISCDI)  projects number  PN-III-P4-ID-PCE-2020-2843  ( EVO-DEVO-CAVE ).
JWA acknowledges the support from a grant from the  User Support Program Space Research  (grant  ALW-GO/13-09 ) of the  Netherlands Organization for Scientific Research (NWO). MB was funded through the  German Research Foundation (DFG)   Eigene Stelle  project  BI 1987/2-1. The computational resources for the assembly of the Thiovulum genome were provided to J. -F. Flot by the  Consortium des Équipements de Calcul Intensif (CÉCI)  funded by the  Fonds de la Recherche Scientifique de Belgique (F. R. S. -FNRS)  under Grant No. 2. 5020. 11;DI was funded through the  German Research Foundation (DFG)   Eigene Stelle  project  IO 98/ 3-1. JM Volland was supported by  Gordon and Betty Moore Foundation  grant  GBMF7617. The work (proposal:  10. 46936/10. 25585/60001074 )conducted by the  U. S. Department of Energy Joint Genome Institute  (https://ror.
DATA AVAILABILITY All sequenced data from Movile Cave is available under bioproject accession number PRJNA673084. The genome of Thiovulum sp. from Frasassi is available under bioproject accession number PRJNA846597. The genomes of Ca. T. imperiosus, and Ca. T. karukerense are available under bioproject accession number PRJNA830902. All other data are provided as part of this manuscript and its supplementary material, including annotation results of all genomes that are provided as supplementary datasets.