Sample collection and culturing<\/p>\n
The biofilm samples were collected from a gas production well, KTG3, in the Southern Kanto gas field (Chiba Prefecture, Japan). The biofilm was located on a wall of the separator discharge outlet at a depth of 20\u201340\u2009cm below the water surface. Samples were taken after eruption of brine on 29 August 2019, 16 July 2020 and 11 May 2023 (Extended Data Fig. 1<\/a>). The gas field is located in an old marine forearc basin and is based on subsurface gas reservoirs in turbidite sandstones, making it geologically and geochemically relevant to past seawater and marine sediment66<\/a>. The brine from the well KTG3 is rich in organic matter, iodine and iron. Detailed site information has been described previously37<\/a> and geochemical characteristics of the brine in the well KTG3 are provided in Extended Data Fig. 1b<\/a>. The collected biofilm sample was placed in a 100\u2009ml glass bottle filled with brine (no head space) and brought back to our laboratory under refrigerated conditions. The samples were also inoculated in 50\u2009ml glass vials containing anaerobic medium onsite (Supplementary Note\u00a02<\/a>). Our iTag analysis showed that the biofilm contained abundant and diverse microorganisms, mostly affiliated with uncultured or cultivation-limited microbial groups, including Promethearchaeota. To isolate new Promethearchaeota strains from the samples, we used a strategic batch-type cultivation35<\/a> that incorporated previous knowledge and our expectations for this archaeal group (Supplementary Note\u00a01<\/a>). Detailed information about the cultivation, including the medium composition and purity checks for the HC1 and SC1 cultures, is provided in the Supplementary Methods<\/a>. Unless otherwise indicated, the HC1 and SC1 cultures were grown in casamino acids\u2013peptone\u2013yeast extract medium (each at 0.005%, w\/v) and casamino acids\u2013powdered milk\u2013polypeptone peptone\u2013yeast extract medium (each at 0.005%, 0.005%, 0.01% and 0.05%, w\/v), respectively, with an inoculum size of 16% (v\/v) for HC1 and 20% (v\/v) for SC1. To measure the dissolved oxygen concentrations in KTG3 culture medium lacking reducing agents, we used a fibre optic oxygen meter (Fibox 3 trace, Presence) in combination with an oxygen sensor (SP-PSt6-YAU). The SP-PSt6-YAU sensor was affixed to the inner wall of the 50\u2009ml serum glass bottles using oxygen-permeable silicone adhesive.<\/p>\n Growth monitoring using qPCR and iTag analysis<\/p>\n DNA extraction, PCR mixture preparation and quantitative PCR (qPCR) analysis were performed as described previously4<\/a>,67<\/a>. PCR primers for qPCR used in this study were designed using the probe design tool of ARB program68<\/a>. The specificity of new primers was confirmed using BLAST and the ARB-SILVA database (https:\/\/www.arb-silva.de\/)<\/a>69<\/a>. The primers and PCR conditions are shown in Supplementary Table 23<\/a>. The doubling times of cultured archaea were calculated based on the semilogarithmic plot of the qPCR data. For iTag analysis, a universal primer pair 530F\/907R, which contained overhang adapters at 5\u2032 ends, was used70<\/a>. The procedures used for library construction, sequencing and data analysis were described previously4<\/a>. All samples for data acquisition were checked for growth and microbial composition by qPCR and iTag before or after the experiments.<\/p>\n Evaluation of growth temperature and substrate use<\/p>\n To evaluate the temperature range for growth, a pure co-culture of HC1 and a tri-culture of SC1 were incubated at a range of 10\u2009\u00b0C to 45\u2009\u00b0C or 10\u2009\u00b0C to 50\u2009\u00b0C. All incubations were performed in triplicate. After a certain period of incubation, 16S rRNA gene copy numbers of each strain were evaluated by qPCR.<\/p>\n To examine substrate use of the cultured archaeal strains, the basal media were supplemented with each substrate sterilized either by autoclaving or filtration. We performed two successive transfers under the same cultivation conditions and conducted qPCR targeting the archaeal strains\u2019 16S rRNA genes in both the first and second subcultures. All incubations were performed in duplicates at the optimum growth temperature for each strain.<\/p>\n To assess the effect of non-proteinaceous substrates on growth of archaeal strains, we added individual substrates to culture medium containing casamino acids\u2013peptone\u2013yeast extract for HC1 and casamino acids\u2013powdered milk\u2013polypeptone peptone\u2013yeast extract for SC1. Each cultivation condition was performed in duplicate.<\/p>\n Epifluorescence microscopy and cell size measurements<\/p>\n An epifluorescence microscope (BX53, EVIDENT) equipped with a phase-contrast objective lens and a colour CCD camera system (DP80 or DP75, EVIDENT) was used for routine microscopy observation, with image acquisition and recording performed using cellSense v.4.3 (EVIDENT). When it was necessary to obtain images without cell movement due to Brownian motion, agarose-coated glass slides were used71<\/a>. Cell and protrusion diameters of HC1 were measured from phase-contrast images using Fiji, an image analysis software based on ImageJ72<\/a>.<\/p>\n SEM<\/p>\n For SEM analysis, archaeal cells were fixed by adding 2.5% (w\/v) glutaraldehyde and incubating for 2\u2009h under anoxic conditions at their optimum growth temperatures and then overnight at 4\u2009\u00b0C. The sample preparation procedure was performed as described previously3<\/a>. The samples were observed under an extreme-high-resolution focused ion beam (FIB)-SEM (Helios G4 UX, Thermo Fisher Scientific) system operated at 1\u2009kV.<\/p>\n To evaluate the volumes and numbers of HC1 cell bodies, protrusions and membrane vesicles under different yeast extract concentrations, we acquired wide-area SEM images using Maps v.2.5 (Thermo Fisher Scientific), which was installed in the FIB-SEM operation system. Wide-area images were obtained by stitching 100 images into a larger 10\u2009\u00d7\u200910 image, with each image taken with an acceleration voltage of 1\u2009kV, a beam current of 13 pA, a field of view of 36.3\u2009\u00b5m \u00d7 24.2\u2009\u00b5m and an image resolution of 23.6\u2009nm per pixel. The acquired wide-area SEM images were then analysed using AMIRA 3D Pro 2024.1 and 2024.2 software (Thermo Fisher Scientific). First, the images were denoised by applying curvature-driven diffusion73<\/a>. Next, HC1 cell bodies, protrusions, membrane vesicles and Methanocalculus cells were segmented with the aid of a deep-learning module based on PyTorch74<\/a> and Keras (https:\/\/github.com\/fchollet\/keras)<\/a>. The deep learning model created from 2,000\u2009\u00d7\u20092,000 pixels was adapted to the full images, followed by manual modification and watershed segmentation. Each segmented component, except for protrusions, was labelled and the area and perimeter were analysed. Protrusions were skeletonized using the Auto Skeleton module75<\/a>, and the number of fibre segments and their curved length were measured after excluding looped fibres or those shorter than 10 pixels.<\/p>\n TEM<\/p>\n For the preparation of ultrathin sections of archaeal cells for TEM, we followed a previously established protocol3<\/a>. The section samples were observed using a conventional TEM (Tecnai 20, FEI) system at an acceleration voltage of 120\u2009kV. Detailed information about observations of HC1\u2019s membrane vesicles and thin filaments on the surface of SC1 protrusions is described in the Supplementary Methods<\/a>.<\/p>\n QFDE-EM<\/p>\n For QFDE-EM observations, archaeal cultures (440\u2009ml of HC1 culture and 25\u2009ml of SC1 culture) were concentrated to approximately 5\u2009ml through filtration using a 0.22-\u00b5m-pore-size polyethersulfone (PES) filter unit (Corning) and transferred to a glass vial in an anaerobic chamber (95:5 (v\/v) N2:H2 atmosphere; COY Laboratory Products). To remove the H2, which is inhibitory to the strains, this vial\u2019s headspace was replaced with N2\/CO2 gas (80:20, v\/v) through flushing. The subsequent procedure was performed as described previously3<\/a>,76<\/a>.<\/p>\n Cryo-EM<\/p>\n Before EM observation, 2\u2009ml of culture liquid was taken from a culture bottle using a syringe, and the culture liquid was concentrated to approximately 50\u2009\u00b5l by centrifugation at 20,400g for 10\u2009min at 20\u2009\u00b0C. Then, 3\u2009\u00b5l of the concentrated liquid culture was applied onto a glow-discharged Quantifoil Mo grid R2\/1 (Quantifoil MicroTools) after mixing with 15\u2009nm gold colloids (Electron Microscopy Sciences) and plunge-frozen into liquid ethane using the Vitrobot Mark IV (Thermo Fisher Scientific) at 4\u2009\u00b0C and 95% humidity. The frozen grids were screened on a 626 side-entry cryoholder (Gatan) with a 200\u2009kV JEM2200FS electron microscope (JEOL). Suitable grids were processed for cryo-ET in the 300\u2009kV Titan Krios G4 electron microscope (Thermo Fisher Scientific) equipped with a bottom-mounted Selectris X energy filter (slit width, 5\u2009eV) and a Falcon 4i direct detector camera. Tile series of each specimen were collected at a nominal magnification of \u00d719,500, resulting in an image resolution of 6.19\u2009\u00c5 per pixel, and a dose symmetric tile scheme was performed between +60\u00b0 and \u221260\u00b0 at 2\u00b0 increments using Tomography software implemented on a TITAN Krios G4, resulting in a total electron dose of 91.5 electrons per \u00c52. The tilt series was aligned with gold fiducials, and tomograms were reconstructed using filtered back projection or SIRT in the IMOD software77<\/a> with an image binning of 5. Segmentations of HC1\u2019s cellular components were performed using AMIRA 3D Pro software. Cellular membranes and granular structures were segmented manually with the help of a deep-learning module, and actin filaments were segmented manually.<\/p>\n Cryo-EM images and tomograms shown in Extended Data Fig. 2i\u2013n<\/a> and Supplementary Video\u00a02<\/a> were acquired using the CRYO ARM 300 II (JEOL) system. Detailed procedures for sample preparation and imaging conditions are provided in the Supplementary Methods<\/a>.<\/p>\n Immunofluorescence staining using lokiactin-specific antibodies<\/p>\n The protocol for antibody preparation and observation by immunofluorescence for lokiactin generally followed that of a previous study5<\/a>. The lokiactin antibodies were made using peptides specific for the lokiactin sequences of HC1 and SC1 (Supplementary Fig. 3<\/a>). The regions of the peptides are the same region as the peptides for \u2018Ca. L. ossiferum\u20195<\/a>. We made three antibodies in total (that is, ab1 for HC1; ab1 for SC1; and ab2 for both strains), and all of them were created using the service of Cosmo Bio and validated by enzyme-linked immunosorbent assay. Detailed procedures of western blotting and immunofluorescence staining are provided in the Supplementary Methods<\/a>.<\/p>\n HCR-FISH<\/p>\n Cells were fixed under anoxic conditions by adding paraformaldehyde solution (methanol-free, formaldehyde solution, 16% (w\/v), Electron Microscopy Sciences) to a final concentration of 4% (v\/v) directly into the culture bottles and incubating for 2\u20133\u2009h at 4\u2009\u00b0C. After fixation, cells were washed with PBS (153.8\u2009mM NaCl, 2.7\u2009mM Na2HPO4, 1.5\u2009mM KH2PO4 (pH\u20097.2)). Subsequently, cells were stored in 50% ethanol with PBS at \u221230\u2009\u00b0C until use in subsequent experiments. The design of specific probes and hybridization chain reaction (HCR)-FISH were performed in accordance with previously described protocols4<\/a>,78<\/a>. A list of the 16S rRNA-targeted oligonucleotide probes used in this study is provided in Supplementary Table 24<\/a>. Cells were counterstained with DAPI in the ProLong Diamond Antifade Mountant with DAPI (Thermo Fisher Scientific). FISH images were acquired using a STELLARIS STED confocal microscope (Leica) with a \u00d7100\/1.4\u2009NA oil-immersion or \u00d793\/1.3\u2009NA glycerol-immersion objective lens at \u00d72 or \u00d76 zoom, a line-scanning frequency of 400 or 700\u2009Hz and a line average of 3. Fluorescence signals were detected in confocal mode at 405, 499 and 553\u2009nm for excitation, and some images were deconvolved using the adaptive Lightning strategy in LasX software (Leica) when necessary. The 3D video was constructed from 27 sections of 8.17\u2009\u03bcm thick using the LasX software. Line-intensity profiles were obtained using LasX software for each of the DAPI and Alexa Fluor 488 channels along a line through the centre of the target cell.<\/p>\n Chemical analyses<\/p>\n To measure peptides and proteins, culture liquid samples were filtered through a 0.2\u2009\u00b5m pore-size polytetrafluoroethylene filter unit (Merck Millipore) before sample preservation and analysis. The molecular mass distribution of peptides and proteins was determined using an UPLC system equipped with the ACQUITY UPLC Protein BEH SEV column (Waters) as described previously79<\/a>.<\/p>\n Methane concentrations in the headspace of culture bottles were measured by gas chromatography as described previously4<\/a>.<\/p>\n Measurement of haem B was performed according to the previously reported method80<\/a>. For this measurement, archaeal cultures (440\u2009ml of HC1, and 900\u2009ml of P. syntrophicum MK-D1) were first concentrated to about 5\u2009ml using a 0.22-\u00b5m-pore-size PES filter unit (Corning) on a clean bench. The concentrated cell suspensions were then further concentrated by pelleting the cells through centrifugation at 20,400g for 10\u2009min at 4\u2009\u00b0C. The collected cells were washed twice with their respective basal medium without energy substrates.<\/p>\n Nitrate and nitrite concentrations were measured by high-performance liquid chromatography (LC) using the TSKgel SAX column (Tosoh, 35\u2009mm\u2009\u00d7\u20094.6\u2009mm inner diameter) with 3% (w\/v) NaCl as the eluent. The column temperature was maintained at 40\u2009\u00b0C, and detection was performed using a UV-VIS detector (GL-7451, GL Science) as described previously81<\/a>.<\/p>\n Iron concentrations in the brine were determined by quadrupole inductively coupled plasma mass spectrometry (ICP-MS; iCAP Qc, Thermo Fisher Scientific). To prevent metal precipitation before measuring the total iron concentrations, brine samples were acidified with high-purity concentrated HNO3 (TAMAPURE AA-100, Tama Chemical) at a final concentration of 1.5\u2009M. For ICP-MS analysis, the water sample was diluted with 0.3\u2009M HNO3, and four internal standards (Be, Sc, Y and In) were added to correct instrument drift. Mass interferences were eliminated using the kinetic energy discrimination mode with H2 reaction gas for high-sensitivity 56Fe analysis.<\/p>\n Genome, transcriptome, proteome and glycoproteome analyses<\/p>\n Genomic DNA from HC1 and SC1 was extracted according to a previous protocol67<\/a>. Whole-genome sequencing libraries were prepared and sequenced using the Illumina MiSeq system, with additional mate-pair sequencing for HC1 and Oxford Nanopore PromethION sequencing for SC1. Detailed procedures of genome assembly, gene prediction and functional annotation are described in the Supplementary Methods<\/a>.<\/p>\n For RNA-based sequencing analysis, total RNA was extracted from HC1 and SC1 cells collected during their exponential growth phase. Specifically, cells were collected from a 25\u2009ml HC1 culture grown for 25\u2009days and a 50\u2009ml SC1 culture grown for 51\u2009days. To concentrate the cells, the cultures were first filtered through a 0.22-\u00b5m-pore-size PES filter unit (Corning) on a clean bench, reducing the volume to approximately 4\u2009ml, and then centrifuged at 20,400g for 10\u2009min at 4\u2009\u00b0C for further concentration. Derailed procedures of RNA extraction, RNA-seq library preparation, sequencing and quantification of gene expression are provided in the Supplementary Methods<\/a>.<\/p>\n For proteomic analysis, HC1, SC1 and P. syntrophicum MK-D1 were cultured under the conditions described in the Supplementary Methods<\/a>. After confirming growth by qPCR, cells were concentrated using a PES filter unit followed by centrifugation, as done for the samples used in western blotting analysis. The collected cell samples were washed once with the basal medium lacking energy substrates, then suspended in 200\u2009\u00b5l of 100\u2009mM triethylammonium bicarbonate (pH\u20098.6) containing 2\u2009mM phenylmethylsulfonyl fluoride. Cells were disrupted by sonication at 2\u2009\u00b0C using a Q700 sonicator (QSonica) operating at 40\u2009kHz. The protein concentration of the cell-free extract was determined by a Qubit fluorometer (Thermo Fisher Scientific). All subsequent procedures, including LC\u2013tandem MS (LC\u2013MS\/MS) data acquisitions and data analysis, were performed as described previously82<\/a>,83<\/a>.<\/p>\n Glycoproteome analysis was performed using the same culture samples as those used for proteome analysis. The procedure was identical up to the point of cell disruption by sonication, after which half of the protein extract samples were allocated for glycoproteome analysis. Subsequent sample preparation, LC\u2013MS\/MS analysis and data processing were performed as previously described46<\/a>,84<\/a>.<\/p>\n Annotation of ESPs<\/p>\n ESPs in Promethearchaeota genomes were annotated using HMMER v.3.4 (hmmsearch)85<\/a> and HMM profiles are available from InterPro86<\/a>. Hits were determined based on an e-value cut-off of 0.0001 and a bit score above the gathering threshold (GA) specified in each HMM profile. Distant hits were identified based on less stringent thresholds for e value (0.01) and bit score (GA-10). We annotated homologues of elongation factor 2 and DNA polymerase epsilon based on (1) distant homology with PFAM domains found in eukaryotic elongation factor 2 (PF00679, PF14492 and\/or PF03764) and DNA polymerase epsilon (PF00136 and PF03104); (2) alignment using MAFFT (v.7.525; default options)87<\/a> with representatives of bacterial, archaeal and eukaryotic homologues annotated with the above domains in UniProtKB (release 2026_01)88<\/a> and clustered with 60% sequence identity (MMseqs2 v.17.b804f)89<\/a>; (3) trimming with trimAl (v.1.5.0; -gt 0.1)90<\/a>; (4) tree construction with FastTree (v.2.1.11; -lg)91<\/a>; and (5) identification of Promethearchaeota homologues that form a sister group with eukaryotic sequences<\/p>\n Phylogenetic analyses<\/p>\n Promethearchaeota genomes were collected from the Genome Taxonomy Database (GTDB) release 22092<\/a> and GenBank. Genomes not included in GTDB r220 were classified into species-level lineages defined in GTDB or new species based on sequence distances estimated using MASH (v.2.3)93<\/a>. Coding sequences were predicted using Prokka (v.1.14.6)94<\/a> with the default parameters and annotated with eggnog-mapper (v.2.1.12)95<\/a> and the eggnog (v.5.0) database using the default parameters96<\/a>. For constructing a genome-based phylogeny of Promethearchaeota, the following single-copy marker genes defined by GTDB were used for phylogenetic tree construction (Supplementary Table 25<\/a>). Maximum-likelihood estimation of the phylogeny was performed using IQTREE v.2.3.5 with ultrafast bootstrapping97<\/a>, a concatenated alignment of these marker genes trimmed using BMGE v.1.12 (-m BLOSUM30 -g 0.67 -b 3)98<\/a> and a universal distribution model (UDM0064 LCLR) as previously described99<\/a>. Maximum-likelihood estimation of the phylogeny for the 16S rRNA gene was performed using a sequence alignment generated with MAFFT (–maxiterate 1000 –localpair)87<\/a>, the optimum model selected by ModelFinder (SYM\u2009+\u2009R6)100<\/a> and regular bootstrapping. Bootstrap values were recalculated using BOOSTER (v.0.1.2)101<\/a>.<\/p>\n For maximum-likelihood estimation of the phylogeny of target genes, we first collected relevant sequences by using amino acid sequences obtained from Promethearchaeota as the query and searching against genes belonging to the same protein family collected from GTDB r220 species representatives92<\/a>. For this, the MMseqs2 v.16.747c6 search command102<\/a> was used with three iterative searches with different sensitivities from sensitivity level 4 to 8.5 (–min-seq-id 0.3 -c 0.7). To select representative sequences for tree construction, we (1) clustered the collected sequences using the MMseqs2 cluster command102<\/a> (–min-seq-id 0.6 -c 0.7); (2) constructed a rough tree using FastTreeMP (v.2.1.11)91<\/a> with an alignment generated with MAFFT (v.7.525)87<\/a> and trimmed with trimAl v.1.5.0 (-gt 0.1)90<\/a>; and (3) subsampled phylogenetically representative sequences from the tree using Parnas (v.0.1.6)103<\/a>. An alignment for the representatives was generated using MAFFT-DASH (–localpair –maxiterate 100 –dash)104<\/a> and trimmed using trimAl (-gt 0.3). Finally, we used this as an input for maximum-likelihood estimation of the phylogeny using IQTREE v.3.1.1 with either 1,000 ultrafast bootstrap replicates or 100 bootstrap replicates97<\/a>, and used ModelFinder100<\/a> to identify the best-fitting model among C10, C20, C30, C40, C50 and C60 (LG\u2009+\u2009CXX\u2009+\u2009G\u2009+\u2009F). For bootstrap replicates, transfer bootstrap expectation values were also calculated using BOOSTER (v.0.1.2)101<\/a>.<\/p>\n For estimation of the phylogeny of the Cox complex, we used a concatenated alignment of subunits 1, 2 and 3. For all archaeal and bacterial species representatives in GTDB release 220, each subunit was identified with hmmsearch85<\/a> against the corresponding protein families (PF00115, PF00116 and PF00510) using a domain e-value cut-off of 0.001 and coverage of 50%. Subunits were concatenated if the target organism encoded (1) one copy of each subunit across its genome; or (2) multiple copies but organized one copy of each within 10 loci of each other. For organisms outside of Promethearchaeota, Cox1 was clustered with 50% sequence identity and 80% coverage using MMseqs2 easy-cluster to choose representative sets of Cox1, 2 and 3. Alignment for each subunit and tree construction of the concatenated alignment was performed as described above but with ultrafast bootstrap support (-bb 1000).<\/p>\n Metabolic reconstruction and codon analysis<\/p>\n Gene expressivity was calculated using cordon56<\/a> using ribosomal proteins as a reference and the default settings. Peptidotrophy was annotated based on the presence of extracellular peptidases, glutamate dehydrogenases and glycine dehydrogenases with high expressivity (within the top 25% of each genome), all consistently found in isolates specialized in peptide degradation (strains MK-D1, HC1, SC1, Coprothermobacter proteolyticus, Necrotrophus endoclepta, Porphyromonas gingivalis and Salinivirga cyanobacteriivorans). Anaerobic peptidotrophy was annotated if high-expressivity OFORs and FeFe or NiFe hydrogenases were also present. Those with NiFe hydrogenase-containing gene sets homologous to those highly expressed by syntrophic Promethearchaeota isolates or archaeal electron-bifurcating hydrogenase MvhADG\u2013HdrABC were further annotated as H2-evolving. Aerobic peptidotrophy was annotated based on the presence of a high-expressivity Cox (within top 50% of each genome) along with a complete complex III (cytochrome b\/b6 with a cytochrome c, Rieske FeS protein and\/or a copper-binding electron carrier like plastocyanin or halocyanin). Aerotolerant peptidotrophy was annotated based on the presence of a high-expressivity Cbd with type 1 or 2 NADH:quinone oxidoreductase but lacking the ability to transfer all electrons from peptide degradation to oxygen reduction. Different expressivity thresholds were used for peptide metabolism and terminal oxidases based on the observation that expressivity tends to be in the 25\u201350% range for archaea that are known to be aerobic (that is, Halobacteriales).<\/p>\n For prediction of aerobic capacity in Promethearchaeota genomes, classifiers were trained using a dataset of strains with annotated aerobic capacity derived from BacDive105<\/a>. Strains were categorized as obligate aerobes, diazotrophic aerobes, facultative (an)aerobes and obligate anaerobes (with and without Cox). For classification, obligate aerobes, diazotrophic aerobes and facultative (an)aerobes were grouped as organisms with aerobic capacity, whereas obligate anaerobes were grouped as lacking aerobic capacity. Host-associated and phototrophic organisms (including oxygenic and anoxygenic phototrophs such as purple and green bacteria) were excluded, as host-associated oxygen gradients and light-driven metabolism can decouple oxygen exposure from intrinsic respiratory strategies, thereby confounding the relationship between genomic features and metabolic lifestyle. BacDive annotations were manually curated based on literature information and genomic annotation.<\/p>\n For genomic annotation, strains annotated as obligate aerobes in BacDive but encoding pathways for anaerobic respiration (that is, dissimilatory nitrate reduction or sulfate reduction) were reclassified as facultative anaerobes. Dissimilatory nitrate reduction and sulfate reduction were annotated based on the presence of genes with the following PFAM domains: Molybdopterin, Molydop_binding, Nitr_red_alph_N or NapB for dissimilatory nitrate reductases, and APS_reductase_C for adenylyl-sulfate reductase. Diazotrophy was annotated for isolates harbouring genes containing the PFAM domain Oxidored_nitro.<\/p>\n Strains classified as obligate anaerobes without experimental evidence of growth under aerobic conditions were excluded. Likewise, strains classified as facultative anaerobes without experimental evidence of enhanced growth under aerobic conditions were removed as they may rather be oxygen-tolerant anaerobes.<\/p>\n For each genome, OFOR and ODH abundance was calculated as log((n\u2009+\u20091)\/genome size (Mb)), where n corresponds to the number of OFOR or ODH homologues. Genes encoding OFOR and ODH were annotated based on PFAM annotations generated by eggNOG-mapper (v.2.1.2)95<\/a>. OFORs were identified based on the presence of PFAM domain combinations \u2018POR,POR_N\u2019, \u2018PFOR_II,POR_N\u2019, \u2018POR_N,TPP_enzyme_C\u2019 or \u2018PFOR_II,TPP_enzyme_C\u2019. OFORs containing three domains were annotated using hmmsearch and a domain e-value cutoff of 0.00001 and domain coverage cut-off of 0.8. ODH abundance was estimated based on the number of annotated E1 (decarboxylase) subunits (PFAM E1_dh), which determines the substrate specificity and serves as a proxy for the presence of 2-oxoacid dehydrogenase complexes. Genomic frequencies were calculated for the five most discriminative amino acids reported previously61<\/a>, as these capture the majority of the discriminatory signal while reducing redundancy among correlated amino acid frequencies.<\/p>\n Classification of aerobic capacity based on OFOR and ODH abundance or amino acid frequency was performed using LR and LDA implemented in scikit-learn (v.1.8.0). Features were standardized using z-score normalization before model training. LR was performed with balanced class weights to account for class imbalance. Models were evaluated using stratified k-fold cross-validation (StratifiedKFold; k\u2009=\u20095), in which the dataset was partitioned into training and test sets while preserving the proportion of organisms with and without aerobic capacity in each split. Predictions were generated using cross_val_predict, ensuring that each sample was evaluated in a held-out fold. Classification performance was assessed using accuracy and balanced accuracy based on these out-of-fold predictions. For each genome, posterior probabilities were obtained from the fitted models, and the maximum posterior probability was used as a measure of classification confidence. For Promethearchaeota, only genomes with \u226585% completeness and \u22645% contamination were analysed, with two Hodarchaeales MAGs encoding Cox and Complex I (JBSWRZ01) or Cox, Complex I and Complex III (JBSWRG01) as exceptions.<\/p>\n To quantify the similarity between metabolic categories, the centroid of each category was calculated as the mean of the normalized OFOR abundance and normalized ODH abundance across all genomes in that category. Pairwise distances between category centroids were first computed using Euclidean distance in the two-dimensional feature space: dij\u2009=\u2009((xi\u2009\u2212\u2009xj)2\u2009+\u2009(yi\u2009\u2212\u2009yj)2)0.5, where x and y represent the mean values of normalized OFOR abundance and normalized ODH abundance for categories i and j. These distances quantify the separation between metabolic regimes in the feature space and were visualized as a centroid distance matrix.<\/p>\n Moreover, Mahalanobis distances were calculated to account for the covariance structure of the data. For each pair of category centroids i and j, the Mahalanobis distance was computed as: Dij\u2009=\u2009(mui \u2013 muj)T\u2009\u00d7\u2009\u03a3\u20131\u2009\u00d7\u2009(mui\u2009\u2013\u2009muj), where mui and muj are the centroid vectors of the two categories and \u03a3 is the pooled covariance matrix estimated from all genomes in the feature space.<\/p>\n To examine archaea specifically, ODH and OFOR abundances were also comparatively analysed for archaeal isolates. Although archaea capable of aerobic growth often use ferredoxin as an electron carrier, such organisms nearly universally conserve a unique OFOR annotatable with three PFAM domains (PFOR_II, POR_N and POR; domain e-value\u2009\u2264\u20090.0001 and domain coverage\u2009\u2265\u200980% based on hmmsearch) (334 out of 340 analysed isolates; Supplementary Table 18<\/a>). The abundances of ODH and other OFORs were calculated for each archaeal isolate as described above. Promethearchaeota were compared with mesophilic archaea as, for thermophiles, increased temperature decreases oxygen solubility and may therefore influence the interaction between metabolism and oxygen.<\/p>\n Reporting summary<\/p>\n