Acanthopleuribacter pedis KCTC 12899^T isolated from a chiton (Acanthopleura japonica) in Japan was obtained from the Korean Collection for Type Cultures (KCTC). This strain grows well on both Marine Agar 2216 (BD) and MA/10 [1/10 marine broth powder (BD), 3/4 natural sea water, agar powder 1.2 g/L].

Strain M133^T was isolated using serial 1/10 dilution plating on MA/10 from coral Porites lutea collected from the Weizhou Island (109°08′35″ E, 21°03’42″ N), China. This strain grows well on MA/10, but not full-strength Marine Agar 2216(BD).

The whole bacterial biomass was collected from MA/10 or MB/10 [1/10 strength marine broth (BD), 3/4 natural sea water] if not indicated.

For the genome analysis, strain M133^T and A. pedis KCTC 12899^T were incubated in MB/10 (1/10 strength marine broth [BD], 3/4 natural sea water) at 25°C with shaking. After 5 days of incubation, the biomass was collected and used for DNA extraction. Genomic DNA was extracted using the SDS method (Lim et al., 2016). The extracted DNA was visualized using agarose gel electrophoresis and quantified using the Qubit^® 2.0 Fluorometer (Thermo Scientific). The whole genome of strain M133 was sequenced using the PacBio Sequel platform and Illumina NovaSeq PE150, while the genome of A. pedis KCTC 12899^T was sequenced using the Illumina NovaSeq PE150. Genomes sequencing was performed by the Beijing Novogene Bioinformatics Technology Co., Ltd. Libraries for single-molecule real-time (SMRT) sequencing were constructed with an insert size of 10 kb using the SMRT bell TM Template kit (v1.0). Meanwhile, libraries for Illumina sequencing were generated using the NEBNext^® Ultra™ DNA Library Prep Kit (NEB, United States) following the manufacturer’s recommendations. For Illumina NovaSeq PE150 sequencing, raw data were filtered using readfq (v10) to obtain clean data. The specific processing steps were as follows: (1) removal of reads containing low-quality bases (mass value ≤ 20) over 40%; (2) removal of reads containing ≥ 10% N; and (3) removal of reads overlapping with adapters ≥ 15 bp. Clean data were assembled using the SOAP de novo (Luo et al., 2012), SPAdes (Bankevich et al., 2012), and Abyss (Simpson et al., 2009), and the results of these three software were integrated with CISA (Lin and Liao, 2013). For the PacBio Sequel platform, low-quality reads (less than 500 bp) were filtered, and clean reads were assembled using SMRT Link (Ardui et al., 2018) and corrected with Illumina data. Finally, the sequencing depth of A. pedis KCTC 12899^T was 55 × with a Q20 of 97%. For strain M133, the PacBio sequencing depth was 236 × , the Illumina sequencing depth was 73 × with a Q20 of 98%.

Gene search and annotation were performed using GeneMarkS (Besemer et al., 2001) and RASTtk (Aziz et al., 2008) in Kbase (Arkin et al., 2018). Metabolism pathways were annotated using the KASS server (Kanehisa et al., 2004). Annotation of COGs was performed using WebMGA (Wu et al., 2011). Peptidases were annotated using the Hotpep-protease (Bush, 2020) based on the Merops database. Carbohydrate-active enzyme annotation was conducted using the dbCAN meta server (Zhang et al., 2018). GIs and prophages were detected using the IslandViewer 4 webserver (Bertelli et al., 2017). Phaster (Arndt et al., 2016), respectively. The insert sequences were detected using IsSaga (Varani et al., 2011). The results of the detected prophages and insert sequences were manually checked with the GenBank annotations. Meanwhile, the pan-genome was calculated using OrthoMCL (2.0) with default parameters in the KBase server (Arkin et al., 2018). Duplicated genes for each genome were calculated from the pan-genome orthologous results. Intra-genome collinearity analysis was performed using MCScanX (Wang et al., 2012). Secondary metabolite biosynthetic gene clusters (BGCs) were analyzed using the online secondary metabolite search tool antiSMASH (Blin et al., 2021). Gene or gene cluster was visualized using Circos (Krzywinski et al., 2009). The average nucleotide identity was calculated using the online ANI calculator (Yoon et al., 2017). Lastly, the average amino acid identity was calculated using the EzAAI (Kim et al., 2021).

H₂S and N₂ production pathways were checked using transcriptome analysis. For transcriptome analysis, strain M133^T was firstly incubated in shaking MB/10 at 25°C. After 3 days of incubation, 50 ml of inoculum was inoculated into two bottles of 200 ml MB, respectively. Collected one bottle of M133^T using centrifugation after 1 day of shaking incubation at 25°C. Then another bottle of M133^T was changed to static incubation after 1 day of shaking incubation at 25°C until the emergence of black particles. Collected the black-stained biomass using centrifugation. Biomass was frozen using liquid nitrogen and stored at –70°C until RNA extraction. RNA was extracted using the RNAprep Pure Cell/Bacteria Kit (Tiangen, Beijing). Sequencing libraries were generated using NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (NEB, United States) following the manufacturer’s recommendations. Sequencing was performed using the Illumina Novaseq platform at the Beijing Novogene Bioinformatics Technology Co., Ltd. Clean reads were obtained after removing adapter and low quality reads. Reads containing N were removed; Reads containing low-quality bases (mass value ≤ 20) over 50% were removed. Sequencing depth was 95×, and quality score Q20 was 98%. Bowtie2-2.2.3 was used in mapping reads to reference genome (Langmead and Salzberg, 2012). Rockhopper was used to identify novel genes, operon, TSS, TTS and Cis-natural antisense transcripts (McClure et al., 2013). HTSeq was used to count the reads numbers mapped to each gene (Anders et al., 2015), and then FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. When FPKM < 1, no expression is considered. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package (Robinson et al., 2010). Differential expression analysis of two conditions was performed using the DEGSeq R package (Wang et al., 2010).

Taxonomic assignment of strain M133^T was following the standard procedure of prokaryotic microbe identification, namely phenotypic analysis, phylogenetic analysis (Felsenstein, 1981; Saitou and Nei, 1987; Swofford, 1993; Kumar et al., 2016), phylogenomic analysis (Na et al., 2018) and chemotaxonomic analysis as done in Wang et al. (2020).

Coral P. lutea was collected from Weizhou Island and incubated in an aquarium at Guangxi University. For the autotrophic experiments, coral fragments were incubated in axenic natural seawater under natural light at room temperature (approximately 25°C) without feeding. During autotrophic incubation, coral fragments were washed with axenic natural seawater once a week to replace the seawater used during incubation. For the preparation of axenic natural seawater, natural seawater was sterilized at 121°C for 20 min, followed by the injection of axenic air for 20 min to restore the carbonate system.

Coral holobiont DNA was extracted from 1 cm × 1 cm coral tissue using the TaKaRa MiniBEST Universal Genome DNA extract kit (v5.0). The bacterial community was analyzed using 16S rRNA gene libraries generated with the 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) primers following the procedure described by Chen et al. (2021).

Acidipila rosea DSM 103428^T (GCF_004339725), “Acidisarcina polymorpha” SBC82 (GCF_003330725), Acidobacterium capsulatum ATCC 51196^T (GCF_000022565), Bryocella elongata DSM 22489^T (GCF_900108185), Candidatus Solibacter usitatus Ellin6076 (GCF_000014905), Edaphobacter modestus DSM 18101^T (GCF_004217555), Granulicella pectinivorans DSM 21001^T (GCF_900114625), Occallatibacter savannae AB23^T (GCF_003131205), Silvibacterium bohemicum S15^T (GCF_001006305), Terracidiphilus gabretensis S55^T (GCF_0014 49115), Terriglobus roseus DSM 18391^T (GCF_000265425), Terriglobus saanensis SP1PR4^T (GCF_000179915), Bryobacter aggregatus MPL3^T (GCF_000702445), Paludibaculum fermentans P105^T (GCF_015277775), Pyrinomonas methylaliphatogenes K22^T (GCF_000820845), Chloracidobacterium thermophilum B^T (GCF_000226295), Holophaga foetida DSM 6591^T (GCF_000242615), Geothrix fermentans DSM 14018^T (GCF_000428885), Thermoanaerobaculum aquaticum MP-01^T (GCF_000687145), and Luteitalea pratensis DSM 100886^T (GCF_001618865).

The GenBank/EMBL/DDBJ accession numbers for the 16S rRNA gene sequence and the whole genome sequence of strain M133^T are MN908335 and CP071793, respectively. The GenBank/EMBL/DDBJ accession number for the whole genome sequence of strain Acanthopleuribacter pedis KCTC 12899^T is JAFREP000000000. The raw 16S rRNA gene reads of bacteria from the coral P. lutea are deposited in the NCBI Sequence Read Archive database (SRA) under BioProjects PRJNA786650 and PRJNA787388.

The complete genome sequence of strain M133^T is 11,786,365 bp long with a G + C content of 60.2 mol%. The genome contains three 16S-23S-5S rRNA operons and 57 tRNA genes. A total of 6,923 protein-coding genes were identified. Meanwhile, the whole genome sequence of A. pedis KCTC 12899^T consists of 10,848,621 bp genome sequences obtained from 98 scaffolds, with a G + C content of 57.3 mol%. Five 5S rRNA, one 16S rRNA, one 23S rRNA, and 65 tRNA genes were identified. In total, 6,716 protein-coding genes were identified. Other acidobacterial genome sequences of type or Candidatus species with genome sizes of 2.66 to 9.97 Mb were obtained from GenBank (Table 1).

Based on the available sequences in GenBank, acidobacterial genome sizes range from 0.55 to 11.88 Mb, of which the completed genome sizes range from 2.32 to 11.79 Mb. The largest acidobacterial genome before this study belonged to Candidatus S. usitatus Ellin6076 (9.97 Mb) that was proposed to arise via horizontal gene transfer, followed by widespread small-scale gene duplication (Challacombe et al., 2011; Challacombe and Kuske, 2012). Thus far, strain M133 has the largest genome (11.79 Mb) among the completely sequenced acidobacteria, followed by A. pedis KCTC 12899^T (10.85 Mb).

To explore the mechanism underlying large genomes, the pan-genome was calculated using 21 selected acidobacteria. Proteins from the 21 selected acidobacteria were clustered into 29,763 families, of which only 418 families belonged to the core genome, suggesting the high diversification of these acidobacteria. The genome of strain M133^T contains 1,355 unique protein families (species-specific; total length of encoding genes is approximately 1.69 Mb) and 4,737 shared protein families (orthologs). The core genome of strain M133^T and A. pedis KCTC 12899^T is approximately 7.44 Mb, encoding 4,246 families of proteins, which is larger than most of the selected acidobacterial genomes (Table 1), suggesting that the last common ancestor of strain M133^T and A. pedis KCTC 12899^T had already obtained a large genome.

Horizontal gene transfer (HGT) events are important in acidobacteria for carbohydrate utilization (Naumoff and Dedysh, 2012; Naumoff, 2016), heavy metal resistance, iron uptake, secondary metabolism, and antibiotic resistance (Gonçalves and Santana, 2021). To explore HGT events, genomic islands (GIs) were explored in selected complete acidobacterial genomes (data from draft genomes were listed as references, which may underestimate the actual GIs, prophage, and insert sequence level). The quantities of GIs were highly variable among the genomes, while no habitant or taxonomic dependence was found (Table 1). The total lengths of GIs in an individual genome range from 0.08 to 0.87 Mb, accounting for 3.01–12.34% of the genome. Mobile genetic elements (MGEs), agents that affect DNA movement in HGT (Frost et al., 2005), were also analyzed. A total of 44 prophages were detected in the selected genomes, wherein 1–4 prophages were detected in an individual genome (Table 1). Interestingly, insert sequences, which were detected in all 21 genomes, were more abundant than prophages in acidobacterial genomes, ranging from 2–231 in an individual genome (Table 1). After manual verification, 44 GIs, 4 prophages, and 87 insert sequences were detected in the genome of strain M133^T. The total length of 44 GIs accounted for approximately 7.38% of the genome. Structurally, all four prophages had corresponding GIs in the genome. A few insert sequences had no corresponding GIs (Figure 1), indicating that these insert sequences may be species-specific. Some prophages and insert sequences co-localized in the same genomic region (Figure 1), suggesting gene transfer and/or recombination events between mobile genetic elements (Eichorst et al., 2018).

Gene duplications are common and highly variable among the selected acidobacterial genomes. Duplicated genes range from 218 to 2,076 for an individual genome, with a total length of 0.28–3.12 Mb, which accounts for 10.33–30.19% of the genome (Table 1). The total length of duplicated genes in the genome of strain M133^T is 3.11 Mb, which accounts for 26.38% of the genome. The duplicated genes in strain M133^T were detected genome-wide (Figure 1), mainly functioning in membrane transport, secondary metabolite biosynthesis, MGEs, and others (Supplementary Table 1). Furthermore, these duplicated genes suggest that there were also widespread gene duplications after HGT events in strain M133^T similar to Can. S. usitatus Ellin6076 (Challacombe et al., 2011). In addition to MGEs, some species-specific genes were duplicated, such as the rRNA genes (Supplementary Figure 1). In total, three rRNA operons were detected in the genome of strain M133, with one orthologous operon with tRNA-Ile-tRNA-Ala inserted between the 16S and 23S rRNA genes, which is consistent with most of the acidobacteria except for Bryobacteraceae (tRNA-Ala-tRNA-Ile). The other two operons have lost the tRNA-Ile-tRNA-Ala insertion and are located in different strands. The 5S rRNA end of the reverse-strand operon is linked to a site-specific tyrosine recombinase, indicating that the duplicated rRNA operons were from site-specific recombination.

Biosynthetic gene clusters (BGCs) usually result from the combination of HGT, duplication, and rearrangement events during complex evolutionary processes (Lopez, 2003; Gallagher and Jensen, 2015; Boutanaev and Osbourn, 2018). The results indicate that BGCs are as large as GIs in strain M133^T (Figure 1) and A. pedis KCTC 12899^T, accounting for 22.6% and 20% of the genome, respectively (Table 1). For other acidobacteria, BGCs account for less than 6.5%. Some prophages and/or insert sequences were co-localized with BGCs in the same genome regions (Figure 1), suggesting that MGEs have affected the evolution of these BGCs.

The total length of GIs, duplicated genes, and BGCs (plus overlaps) is 5.12 Mb in the genome of strain M133^T, which accounts for 43.43% of the genome size. Meanwhile, the total length of these in Can. S. usitatus Ellin6076 and P. fermentans P105^T are 3.70 and 3.21 Mb, accounting for 37.11% and 33.86% of the genome, respectively (Table 1). Therefore, HGT and gene duplication have impacted the long-term evolution of these large genome acidobacteria.

The quantity of clusters of orthologous groups (COGs) for acidobacteria ranged from 1260 to 1962, and larger genomes encoded more COGs (Supplementary Table 2). According to N class COGs of A. capsulatum ATCC 51196^T (Kishimoto et al., 1991), 17 of the 21 selected acidobacteria should be motile (Supplementary Table 2), although some of them were reported to be non-motile (Liesack et al., 1994; Eichorst et al., 2007; Kulichevskaya et al., 2010, 2014; Pankratov and Dedysh, 2010; Okamura et al., 2011; Dedysh et al., 2012; Crowe et al., 2014; Vieira et al., 2017). C. thermophilum B^T, G. fermentans DSM 14018^T, and T. aquaticum MP-01^T contain 2–7 families of N-class COGs, which are consistent with the non-motile reports (Coates et al., 1999; Losey et al., 2013; Tank and Bryant, 2015).

Meanwhile, the quantities of G-class COGs varied among habitats (Supplementary Table 2), which decreased in the sequence of soil acidobacteria, animal-associated acidobacteria (A. pedis KCTC 12899^T and strain M133^T), and microaerophilic or anaerobic acidobacteria (H. foetida DSM 6591^T, C. thermophilum B^T, G. fermentans DSM 14018^T, and T. aquaticum MP-01^T), indicating the adaptive carbohydrate utilization abilities of these strains (Ward et al., 2009; de Castro et al., 2013; Belova et al., 2018; Eichorst et al., 2018; Kalam et al., 2020). The Q-class of COGs (Supplementary Table 2) indicated that strain M133^T and A. pedis KCTC 12899^T harbors genes involved in secondary metabolite biosynthesis, transport, and catabolism compared with other acidobacteria (details are displayed in the section “Secondary Metabolism”).

Carbohydrate-active enzymes (CAZymes) prediction results (Supplementary Table 3) revealed that facultative anaerobic P. fermentans P105^T and Can. S. usitatus Ellin6076 have the most abundant families of CAZymes (111–112), followed by aerobic soil acidobacteria (75–103), animal-associated acidobacteria (56–63), and anaerobic bacteria (23–35). Glycoside hydrolases, which are involved in hydrolysis and/or rearrangement of glycosidic bonds, contributed to the differences in CAZyme abundance among the genomes (Supplementary Table 3). Soil acidobacteria have 40–71 families of glycoside hydrolases, while animal-associated and anaerobic acidobacteria have very few glycoside hydrolase families (23–24 and 7–11, respectively). These results further confirmed the outstanding carbohydrate utilization ability of soil acidobacteria (Ward et al., 2009; de Castro et al., 2013; Belova et al., 2018; Eichorst et al., 2018; Kalam et al., 2020). Interestingly, only 48 CAZymes were detected in the aerobe Bryobacter aggregatus MPL3^T isolated from acidic Sphagnum peat bogs, of which 18 families belonged to glycoside hydrolases.

In this study, Acidobacteria are observed to be versatile in carbohydrate degradation. Genes encoding for chitinases (GH18, 23) involved in chitin degradation; α-amylases (GH57), glucoamylases (GH15), and glycogen/starch phosphorylases (GT35) in starch degradation; endoglucanases (GH51, 9, and 10) and β-glucosidases (GH3) in cellulose degradation, and xylanases (GH10, 51) in xylan degradation were detected in most of the acidobacterial genomes (Supplementary Table 3). Previously, some acidobacteria have been shown to degrade a range of carbohydrates. For example, soil bacteria Silvibacterium bohemicum S15^T (Lladó et al., 2016) and “Acidisarcina polymorpha” SBC82 (Belova et al., 2018) have been reported to degrade starch, cellulose, chitin, xylan, and pectin. Genes encoding carbohydrate hydrolases indicated that strain M133^T and A. pedis KCTC 12899^T can degrade chitin (mainly from the cell wall of fungi and the exoskeleton of arthropods) using N-acetylglucosaminidase (GH20, Meekrathok et al., 2021), but not pectin (mainly from plants) using polygalacturonases (GH28) and α-rhamnosidases (CBM67, GH106).

Furthermore, acidobacteria are versatile in polysaccharide biosynthesis. Starch is synthesized by the starch synthase and 1,4-alpha-glucan branching enzyme. Genes encoding starch synthases (GT5) were detected in 20 of the selected acidobacterial genomes, excluding Geothrix fermentans DSM 14018^T, while those encoding the 1,4-alpha-glucan branching enzyme (CBM48 + GH13_9) were detected in only seven genomes (Supplementary Table 3). These observations indicate that most of the acidobacteria can synthesize amylose, while only a few (including strains M133 and A. pedis KCTC 12899^T) can produce glycogen. In particular, Acidobacteriaceae cannot produce glycogen.

Cellulose is important in biofilm formation and cellular water holding and is synthesized by various cellulose synthases. To date, only the genomes of A. capsulatum and Terriglobus saanensis have been reported to biosynthesize cellulose using bcs operons (Kalam et al., 2020). Genomic detection using the Kbase platform and conserved domains search in GenBank indicated that seven acidobacteria in the family Acidobacteriaceae, namely A. capsulatum ATCC 51196^T, T. saanensis SP1PR4^T, A. rosea DSM 103428^T, B. elongate DSM 22489^T, G. pectinivorans DSM 21001^T, T. gabretensis S55^T, and “A. polymorpha” SBC82, contain bcs operons. The distribution of the bcs operon was variable in the genus Edaphobacter, and no bcs operon was detected in the genome of E. modestus DSM 18101^T. Cellulose synthase was not detected in strain M133^T and A. pedis KCTC 12899^T.

According to the catalog of MEROPS peptidase, the genomes of selected acidobacteria contain 45–95 families of peptidases (Supplementary Table 4). Species from the same family have similar numbers of peptidase families. The anaerobic Holophagaceae foetida DSM 6591^T has the lowest number of peptidase family (45). Acidbacteriaceae had 57–77 peptidase families, followed by animal-associated Acanthopleuribacteraceae (79–84), and Bryobacteriaceae has the highest number of peptidase families (79–95). Other acidobacteria have 62–87 families of peptidases.

Bacterial collagenase H (M09B), imelysin (M75), and MtfA peptidase (M90) were only detected in strain M133 and A. pedis KCTC 12899^T (Supplementary Table 4). Moreover, collagenase (U32), subtilisin (S08A), pappalysin-1(M438), astacin (M12A), adamalysin (M12B), serralysin (M10B), fungalysin (M36), lysostaphin (M23B), immune inhibitor A peptidase (M06), thermolysin (M04), and bacteriocin-processing peptidase (C39) were detected in strain M133^T and A. pedis KCTC 12899^T. In particular, a high abundance of subtilisin (11–12) and thermolysin (8–11) coding sequences were detected in strain M133 and A. pedis KCTC 12899^T (Supplementary Table 4). Membrane dipeptidase (M19), carnosine dipeptidase II (M20F), glutamate carboxypeptidase (M28B), aminopeptidase M29 (M29), streptogrisin A (S01E), carboxypeptidase Y (S10), SpoIVB peptidase (S55), Xaa-Pro dipeptidyl peptidase (S15), and sedolisin (S53) were absent in strain M133^T and A. pedis KCTC 12899^T, whereas they are widely detected in other acidobacteria.

Genes encoding enzymes that involved in peptide assembly, regulation, resistance, and synthesis of secondary metabolites are usually physically clustered into biosynthetic gene clusters (BGCs, Medema et al., 2015). More than 40 BGCs were detected in genomes of strain M133^T and A. pedis KCTC 12899^T, respectively (Table 1). According to the antiSMASH prediction, these BGCs may function in non-ribosomal peptides, polyketides, terpenes, thioamitides, arylpolyenes, and ranthipeptides biosynthesis. While other acidobacterial genomes contain only a few BGCs (0–13) (Table 1). Genes of non-ribosomal peptide synthetases and type I polyketide synthetases accounted for more than 40% of the BGCs in strain M133^T and A. pedis KCTC 12899^T. The exact secondary metabolites of strain M133^T and A. pedis KCTC 12899^T need further study.

Bacterial long-chain polyunsaturated fatty acids (PUFAs) are synthesized by multi-domain protein complexes akin to type I iterative fatty acids and polyketide synthases (Shulse and Allen, 2011). The typical PUFA-synthesizing gene cluster was discovered in marine gammaproteobacteria, in which five genes, pfaABCDE, participated in the de novo biosynthesis of PUFAs (Shulse and Allen, 2011). EPA was detected in both strain M133^T and A. pedis KCTC 12899^T using MIDI gas chromatography, with EPA contents of 6.2% and 12.2%, respectively (Supplementary Table 5). Among the five genes in the EPA-synthesizing gene cluster, four genes (pfaB and pfaC were annotated as one gene: pfaB/C) were detected in strain M133^T and A. pedis KCTC 12899^T, wherein three genes formed the cluster pfaDAB/C, while pfaE encoding for phosphopantetheinyl transferase was located far from this cluster (Supplementary Figure 2). This EPA biosynthesis cluster is much similar to those in deltaproteobacteria “Desulfococcus oleovorans” Hxd3 and Sorangium cellulosum So ce 56 (S. bulgaricum Soce 321^T was reported to produce EPA, Mohr et al., 2018), but not the EPA-producing gammaproteobacteria Shewanella pealeana ATCC 700345^T (Shulse and Allen, 2011). None of the pfa genes were detected in the other acidobacterial genomes.

There are several types of respiration electron acceptors in acidobacteria, such as oxygen, nitrate, Fe (III), Mn (IV), and sulfur. Most acidobacteria, mainly Acidobacteriaceae, are aerobic or microaerobic. The A-type cytochrome c oxidase was ubiquitous in Acidobacteriaceae, while the high-oxygen affinity cbb3-type cytochrome c oxidase and cytochrome bd ubiquinol oxidase were detected in only half of the selected genomes (Supplementary Table 6). This illustrates that some of Acidobacteriaceae can respond to microoxic concentrations (Morris and Schmidt, 2013; Eichorst et al., 2018). Meanwhile, only the cbb3-type terminal oxidase has been detected in the genome of the microaerophile C. thermophilum (Tank and Bryant, 2015). Both strain M133 and A. pedis KCTC 12899^T contain the A-type cytochrome c oxidase, cbb3-type cytochrome c oxidase, and cytochrome bd ubiquinol oxidase, suggesting that these two strains can survive under variable oxygen conditions.

Nitrate is an energy-efficient oxygen substitute under anaerobic conditions. However, it seems that acidobacteria are not adept at nitrate reduction. Genes for assimilatory nitrate reduction were detected in five of the 21 genomes, and those for nitrite reduction were detected in six genomes. Dissimilatory nitrate reduction genes were detected in four genomes, namely G. ferrnentans DSM 14018^T, T. aquaticum MP-01^T, A. pedis KCTC 12899^T, and strain M133, while those for nitrite reduction were found in five genomes. These results may be due to a few assimilatory nitrate reductases involved in dissimilatory nitrogen metabolism, as suggested by Morozkina and Zvyagilskaya (2007); Eichorst et al. (2018). Denitrification genes nirS and norBC were detected in G. ferrnentans DSM 14018^T, T. aquaticum MP-01^T, P. fermentans P105^T, Luteitalea pratensis DSM 100886^T, and A. pedis KCTC 12899^T, whereas the complete denitrification pathways, which included napA, nirS, norBC, and nosZ were detected only in strain M133^T. Transcriptome analysis confirmed the expression of genes nirS (FPKM 2311: 1.2), norB (1364: 0.6), norC (856: 0), and nosZ (723: 1.3) under anaerobic conditions in strain M133^T. API 20NE results also suggest that strain M133^T can reduce nitrate to nitrogen.

Fe (III) or Mn (IV) reduction is an important anaerobic respiration pathway in acidobacteria. The outer membrane c-type cytochromes OmcA and MtrC (also known as OmcB) have been previously demonstrated to play an important role in the reduction of Fe (III) and Mn (III/IV) (Beliaev et al., 2001; Myers and Myers, 2001). OmcA and MtrC were detected in T. aquaticum MP-01^T, P. fermentans P105^T, and Can. S. usitatus Ellin6076. Both T. aquaticum MP-01^T (Losey et al., 2013) and P. fermentans P105^T (Kulichevskaya et al., 2014) have been reported to anaerobically reduce Fe (III) or Mn (IV). Fe (III) respiration in Can S. usitatus Ellin6076 needs to be confirmed using growth-based tests. The gene for MtrF, a homolog of MtrC, was detected in the genome of G. ferrnentans DSM 14018^T, whereas no gene for OmcA was detected, although G. ferrnentans DSM 14018^T was confirmed to use Fe (III) as an electron acceptor (Coates et al., 1999).

Oxidized sulfur compounds are relatively low-efficient electron acceptors during respiration. Genes encoding assimilatory sulfate reductases were detected in 13 of the 21 acidobacterial genomes, and those encoding dissimilatory sulfate reductases were detected in two genomes (Supplementary Table 6). Genes encoding for sulfite reductase were detected in 16 genomes, of which only one (H. foetida DSM 6591^T) was dissimilatory sulfite reductase (Supplementary Table 6). Only assimilatory sulfate and sulfite reduction genes were detected in strain M133 and A. pedis KCTC 12899^T.

The bifunctional enzyme sulfhydrogenase was first reported in the anaerobic archaeon Pyrococcus furiosus, and this enzyme produces sulfide hydrogen from elemental sulfur or polysulfide under anaerobic conditions (Ma et al., 1993). Complete sulfhydrogenase genes hydBGDA (Pedroni et al., 1995) and accessory genes hypCDEF were also detected in the genome of strain M133. Therefore, sulfhydrogenase may function as a secondary anaerobic respiration pathway in this strain. High nutrition, especially tryptone or other protein derivatives, promotes hydrogen sulfide production. However, the actual substrate of sulfhydrogenase in strain M133 remains unclear. In addition, the amino acid sequence of the M133 sulfhydrogenase β subunit (homolog to sulfite reductase) was searched using NCBI BLAST, and both sulfite reductase and 4Fe-4S dicluster-containing proteins were retrieved. After manually checking the RASTtk re-annotated genomes in Kbase, approximately 31 genomes were confirmed to harbor complete sulfhydrogenase genes (αβγδ) (Supplementary Table 7). Phylogenetically, these sulfhydrogenase-containing strains belonged to approximately 10 bacterial phyla, including Acidobacteria, Cyanobacteria, “Chloroflexi,” “Deinococcus-Thermus,” “Nitrospinae,” Proteobacteria (αβγδ), Bacteroidetes, Actinobacteria, Planctomycetes, and Verrucomicrobia, and two archaeon phyla. Although hydrogenase and sulfite reductase were detected in Firmicutes (Clostridium), no sulfhydrogenase was found. Three archaeon strains from Euryarchaeota and Thaumarchaeota, including Pyrococcus furiosus DSM 3638, were also found to harbor complete sulfhydrogenase genes. This universally and randomly taxonomic distribution pattern suggests that HGT may have contributed to the dissemination of sulfhydrogenase. To examine the possibility of HGT, the G + C content of hydBGDA and the corresponding genome (where the hydBGDA were extracted) were analyzed. Results indicated that the G + C content of hydBGDA is highly variable among genomes (32.84–74.08 mol%), while the G + C content difference between hydBGDA and the whole genome is usually 2–4 mol% (Supplementary Table 7). This difference between a common ancestor and the individual genomes is a long-term evolutionary result. Therefore, the existence of sulfhydrogenase in strain M133 is likely the evolutionary remnant and not a result of recent HGT.

Several strains have been reported to undergo fermentative growth (T. aquaticum MP-01, G. fermentans DSM14018, H. foetida TMBS4, and A. ailaaui PMMR2) (Liesack et al., 1994; Coates et al., 1999; Losey et al., 2013; Myers and King, 2016). However, the function of the possible fermentation enzymes remains to be confirmed (Eichorst et al., 2018). Genes involved in fermentation, such as lactate dehydrogenase, fumarate reductase/succinate dehydrogenase, and phosphate acetyltransferase, were also detected in strain M133 and A. pedis KCTC 12899^T; however, no fermentation growth was detected in the API 50 CH sheet.

To date, species closely related to strain M133^T (>92% 16S rRNA gene identity) have seldom been detected in the environment. To explore these allied species, GenBank databases were blasted using the 16S rRNA gene sequences of strain M133^T and A. pedis KCTC 12899^T. Only two operational taxonomic units (OTUs) (GU319135 and MT858037) were retrieved from the standard nucleotide database. SRA datasets for seawater and algae (Mei et al., 2019), sediment (Zhang et al., 2017; Ghate et al., 2021), coral and coral reefs (Liang et al., 2017; Pearman et al., 2019; Chen et al., 2021), and sponge (Baquiran and Conaco, 2018; Wu et al., 2018) were blasted, and the results indicated that strain M133^T allied reads were occasionally detected in the datasets of the Red Sea coral reef (PRJNA479721) and the Weizhou Island coral P. lutea (PRJNA786650). This occasional appearance of strain M133^T allied reads in high throughput sequencing datasets indicated that these bacteria do exist in related environments, but are usually not detected. Technically, DNA template preparation and PCR amplification can influence the outcome of certain reads. It is possible that the allied species exist in very low abundance in related habitats and/or cells are resistant to universal DNA extraction methods. Another reason may be the bias from the universal primer (Takahashi et al., 2014) and the barcode of indexed-primer (O’Donnell et al., 2016) used in PCR amplification. The exact cause of the low detection of strain M133^T allied species requires further study.

To date, the hosts or habitats of the allied species of strain M133^T have been limited to animals. Strain M133^T was isolated from the coral P. lutea, and five reads (100% identities to strain M133^T) were detected in another coral P. lutea 1 year later (PRJNA786650). Interestingly, the bacterium clone A3M_UNP0_21 (GU319135, with 98% 16S rRNA sequence identity to strain M133^T) was also detected in a coral (Acropora eurystoma in Israel) (Meron et al., 2011). Meanwhile, A. pedis KCTC 12899^T was isolated from the chiton Acanthopleura japonica in Japan (Fukunaga et al., 2008), and its closest bacterial clone OTU20 (MT858037, 100% identity) was detected in the gut of fish Seriola rivoliana in Mexico. Although the exact niches of the allied reads of strain M133^T (92–99% to strain M133^T or A. pedis KCTC 12899^T) are unclear in the PVC plates deployed in the Red Sea coral reefs (Pearman et al., 2019), it is evident that benthic animal species could settle on similar autonomous reef monitoring structures (Leray and Knowlton, 2015), raising the possibility that allied reads were also from animals. In addition, some distant Acanthopleuribacteraceae reads also displayed an intimate association with coral. After three months of autotrophic incubation and weekly axenic washing, the abundance of Acanthopleuribacteraceae reads (85.9–91.7% identities to strain M133^T, PRJNA787388) increased evidently in coral P. lutea, coinciding with the significant bacterial community composition shift (Supplementary Figures 3, 4). As the survivors of long-term autotrophy and repeated axenic wash, these distant reads should also be closely correlated.

In this study, it can be observed that strain M133^T and A. pedis KCTC 12899^T have evolved to associate with animals. Both strain M133^T and A. pedis KCTC 12899^T possess the ability to initiate and maintain symbiotic interactions with animals. Microbial motility and chemotaxis are pivotal for the onset and maintenance of symbiotic interactions (Raina et al., 2019). Moreover, both acidobacteria are mobile by using the flagella (Supplementary Figure 5; Fukunaga et al., 2008). Ankyrins, which was also detected in strain M133^T and A. pedis KCTC 12899^T, can aid bacteria in evading the eukaryotic immune system (Jahn et al., 2019) in bacteria-animal symbiosis. For intracellular symbiotic bacteria, there is usually a dramatic genome reduction (McCutcheon and Moran, 2011; Jäckle et al., 2019), so strain M133^T and A. pedis KCTC 12899^T are not likely to form an intracellular symbiosis with animals. The discovery of OTU20 (MT858037) in the fish gut suggests that strain M133^T and A. pedis KCTC 12899^T may also exist in the animal gut environment. Sulfatases can digest the highly sulfated glycans in the gut, whereas the activation of sulfatases requires post-translational modification catalyzed by radical S-adenosyl-L-methionine proteins (Benjdia et al., 2007, 2011). Both genes encoding for sulfatase and radical S-adenosyl-L-methionine proteins were detected in the genomes of strain M133^T and A. pedis KCTC 12899^T. The low oxygen and nitrate respiration capabilities of these bacteria further facilitate their survival and function in animal gut-like environments. The general metabolic pattern indicates that there is a positive selection resulting from long-term animal associations. Glycoside hydrolases related to plant polysaccharide degradation, such as alpha-L-rhamnosidase (GH106), endo-1,4-beta-xylanase (GH10), and alpha-galactosidase (GH27), are usually absent from strains M133^T and A. pedis KCTC 12899^T. In contrast, animal disease-related virulence factors, such as collagenase (Harrington, 1996; Penttinen et al., 2016), astacin (Ricard-Blum and Vallet, 2016), adamalysin (Ricard-Blum and Vallet, 2016), subtilisin (Imamura et al., 2017; Łagowski et al., 2021), and thermolysin (Kong et al., 2015; Tsaplina et al., 2020) are present in strain M133^T and A. pedis KCTC 12899^T.

Meanwhile, the allied species of strain M133^T may benefit animal hosts in dealing with environmental fluctuations. Several bacteria have been shown to protect their animal hosts from pathogens and heat stress, such as Pseudoalteromonas sp., which inhibits coral pathogens (Vibrio sp.) using antibiotics and peptidase (Offret et al., 2016; Richards et al., 2017; Rosado et al., 2019). Ruegeria sp. inhibits the coral pathogen Vibrio coralliilyticus (Miura et al., 2019), whereas Muricauda sp. protects coral endosymbionts from thermal stress by producing zeaxanthin (Motone et al., 2020). Although the exact products of BGCs are still unclear, the large secondary metabolite biosynthesis capacities of strain M133 and A. pedis KCTC 12899^T are still encouraging. According to antiSMASH known-cluster-blast, BGCs of strain M133 are very similar to the biosynthesis gene clusters of jerangolid A, oocydin A, nematophin, and microsclerodermin M. These compounds have been reported to have cytotoxic activities against fungi and gram-positive bacteria (Gerth et al., 1996; Strobel et al., 1999; Melikhova et al., 2016; Cai et al., 2017). Meanwhile, EPA accounts for 1% of the total fatty acids of the whole P. lutea holobionts; therefore, strain M133, an EPA-producer, may also benefit P. lutea holobionts in dealing with temperature and pressure variation (Valentine and Valentine, 2004).

The closest taxonomic neighbor of strain M133^T, based on the 16S rRNA gene sequence similarity, is A. pedis KCTC 12899^T, they share 92.4% identity (recommended genus demarcation is 95%, Yarza et al., 2014). The average nucleotide identity between strains M133^T and KCTC 12899^T is 70% (recommended genus demarcation is 74%, Barco et al., 2020). The average amino-acid identity between strains M133^T and KCTC 12899^T is 63.97%, which is lower than the recommended genus demarcation of 68% (Konstantinidis and Tiedje, 2005). These indices indicated that strain M133^T and A. pedis KCTC 12899^T belong to different genera. Phylogenetic analysis based on the 16S rRNA gene sequences using maximum-likelihood, neighbor-joining, and maximum-parsimony algorithms (Supplementary Figures 4–6) and phylogenomic analysis using 92 concatenated sequence (Figure 2) indicated that strain M133^T forms a distinct branch beside A. pedis KCTC 12899^T in the family Acanthopleuribacteraceae.

Cells of strain M133^T are Gram-negative, non-spore-forming, motile, aerobic rods (Supplementary Figure 7). This strain can be distinguished from A. pedis KCTC 12899^T by the ability to reduce nitrate to nitrogen and the production of hydrogen sulfide. Chemotaxonomic features can also differentiate these two strains. The major cellular fatty acids (> 5%) of strain M133^T are isoC_15:0 (35.6%), C_16:0 (29.8%), C_15:0 (8.1%), and C_20:5ω3c (6.2%), while those of A. pedis KCTC 12899^T are C_16:0 (29.6%), isoC_15:0 (17.2%), isoC_17:0 (12.4%), C_20:5ω3c (12.2%), and C_16:0N alcohol (5.8%) (Supplementary Table 5). The major polar lipids of strain M133^T are phosphatidylethanolamine, phosphatidylglycerol, and diphosphatidylglycerol, which are similar to those of A. pedis KCTC 12899^T (Supplementary Figure 8). Other features of strain M133^T are listed in Table 2 and the new taxon description.

In summary, phylogenomic, phylogenetic, chemotaxonomic and phenotypic differences to the closest type strain A. pedis KCTC 12899^T indicated that strain M133^T represents a new species in a new genus in the family Acanthopleuribacteraceae, for which Sulfidibacter corallicola gen. nov., sp. nov. is proposed.

Sulfidibacter (sul.fi.di.bac’ter. N.L. neut. n. sulfidum, sulfide; N.L. masc. n. bacter, rod; N. L. masc. n. Sulfidibacter, sulfide producing rod).

Cells are Gram-stain negative, non-spore-forming, motile, aerobic rods. Catalase and oxidase are positive. Nitrate reduction is positive. Major respiratory quinone are menaquinone 6 and 7 (MK6 & 7). Major cellular fatty acids are isoC_15:0, C_16:0, C_15:0 and C_20:5ω3c. Major polar lipids are phosphatidylethanolamine, phosphatidylglycerol, and diphosphatidylglycerol.

Type species is Sulfidibacter corallicola.

Sulfidibacter corallicola (co.ral.li′co.la. L. neut. n. corallum, coral; L. masc. suff. -cola, inhabitant dweller; N.L. n. corallicola, coral-dweller).

Cells have following features in addition to the genus description. Cells are usually 2.3–5.3 μm in length and 0.9–1.5 μm in width. Cells are motile by single polar flagellum. Colonies on marine agar are circular, smooth and yellow in color. H₂S can be produced when incubated in high strength of marine agar or marine broth. Cells can grow at 15–37°C, pH4-10 under 0–10% NaCl (w/v). Starch, tween20, 40 and 60 are degraded. Nitrate is reduced to nitrogen. Arginine dihydrolase, urease, β-glucosidase, β-galactosidase and protease are positive in API 20NE test. Alkaline phosphatase, esterase C4, esterase C8, lipase C14, leucine arylamidase, valine arylamidase, cystine arylamidase, trypsin, chymotrypsin, acid phosphatase, naphthol-AS-BI-phosphohydrolase and N-acetyl-β-glucosaminidase are positive in API ZYM test. D-maltose, D-trehalose, D-cellobiose, gentiobiose, sucrose, turanose, D-salicin, N-acetyl-β-glucosamine, N-acetyl-β-galactosamine, N-acetyl neuraminic acid, D-glucose, D-fucose, L-fucose, inosine, D-arabitol, myo-inositol, glycerol, glucose-6-PO₄, D-fructose-6-PO₄, gelatin, L-aspartic acid, L-glutamic acid, L-histidine, L-pyroglutamic acid, L-serine, D-galacturonic acid, D-gluconic acid, D-glucuronic acid, glucuronamide, quinic acid, D-saccharic acid, D-lactic acid methyl ester, L-lactic acid, citric acid, α-keto-glutaric acid, D-malic acid, Tween40, γ-amino-butyric acid, β-hydroxy-butyric acid, α-keto-butyric acid, acetic acid and formic acid are oxidized in Biolog GenIII microplate. The cellular fatty acids also contain C_14:0, isoC_17:0, isoC_11:0, C_14:1ω7c, isoC_13:0. The polar lipids also contain one phospholipid, five unidentified ninhydrin-positive lipids, and three unidentified polar lipids.

Type strain M133^T (=MCCC 1K03775^T = KCTC 72445^T), was isolated from coral Porites lutea, China. Genome accession in GenBank is CP071793, and the DNA G + C content is 60.2 mol%.