Insights into the genome structure of four acetogenic bacteria with specific reference to the Wood–Ljungdahl pathway

Abstract Acetogenic bacteria are obligate anaerobes with the ability of converting carbon dioxide and other one‐carbon substrates into acetate through the Wood–Ljungdahl (WL) pathway. These substrates are becoming increasingly important feedstock in industrial microbiology. The main potential industrial application of acetogenic bacteria is the production of metabolites that constitute renewable energy sources (biofuel); such bacteria are of particular interest for this purpose thanks to their low energy requirements for large‐scale cultivation. Here, we report new genome sequences for four species, three of them are reported for the first time, namely Acetobacterium paludosum DSM 8237, Acetobacterium tundrae DSM 917, Acetobacterium bakii DSM 8239, and Alkalibaculum bacchi DSM 221123. We performed a comparative genomic analysis focused on the WL pathway's genes and their encoded proteins, using Acetobacterium woodii as a reference genome. The Average Nucleotide Identity (ANI) values ranged from 70% to 95% over an alignment length of 5.4–6.5 Mbp. The core genome consisted of 363 genes, whereas the number of unique genes in a single genome ranged from 486 in A. tundrae to 2360 in A.bacchi. No significant rearrangements were detected in the gene order for the Wood–Ljungdahl pathway however, two species showed variations in genes involved in formate metabolism: A. paludosum harbor two copies of fhs1, and A. bakii a truncated fdhF1. The analysis of protein networks highlighted the expansion of protein orthologues in A. woodii compared to A. bacchi, whereas protein networks involved in the WL pathway were more conserved. This study has increased our understanding on the evolution of the WL pathway in acetogenic bacteria.


| INTRODUC TI ON
Acetogenic bacteria, or acetogens, are obligate anaerobes converting one-carbon substrates, such as carbon dioxide, formate, methyl groups, or carbon monoxide into acetate using molecular hydrogen as electron donor through the Wood-Ljungdahl (WL) pathway, a process known as acetogenesis (Ragsdale & Pierce, 2008). Acetogenesis was first described in the early '30 and has been extensively studied in Clostridia (Drake, 1994). The WL pathway was considered for a long time to be a specific trait of species belonging primarily to the Firmicutes (Ragsdale & Pierce, 2008), but a number of recent studies have shown that this pathway is far more spread in the microbial tree of life than previously thought (Adam, Borrel, & Gribaldo, 2018;Borrel, Adam, & Gribaldo, 2016;Graber & Breznak, 2004;Hug et al., 2013;Strous et al., 2006).
Acetogenic species have been found in the archaeal kingdom, although most Archaea produce methane instead of acetate as end product (Borrel et al., 2016), in Chloroflexi (Hug et al., 2013), Spirochetes (Graber & Breznak, 2004), and Planctomycetes (Berg, 2011;Strous et al., 2006). Due to its low ATP requirement, the WL pathway can be found in prokaryotes adapted to conditions that approach the thermodynamic limits of life (Schuchmann and Mueller, 2014). In addition, comparative genomic analyses of extant microbial taxa revealed that the predicted last common universal ancestor possessed the WL pathway (Adam et al., 2018;Weiss et al., 2016). It is thus conceivable that the WL pathway represented an efficient way to produce energy in the early Earth environment before the great oxidation event, that is the enrichment of oxygen in the early earth atmosphere as a consequence of the emergence of organisms able to perform oxygenic photosynthesis (Poehlein et al., 2012;Weiss et al., 2016). The main advantages of the WL pathway include the following: its versatility; it can be coupled to methanogenesis or to energy conservation via generation of electrochemical gradients; its modularity, since some species utilize partial WL pathways to channel electrons produced during fermentation to CO 2 ; its flexibility, as several organisms use different coenzymes and/or electron carriers, and in some cases the WL pathway is reversed (e.g., it generates molecular hydrogen and carbon dioxide from acetate for energy production (Schuchmann & Mueller, 2016).
There is a growing interest toward acetogens, as they can be used as biocatalyst for the conversion of synthesis gas (a mixture of H 2 and CO and/or CO 2 ) into fuels or chemicals with low energy supply (Bengelsdorf et al., 2016;Cavicchioli et al., 2011;Shin et al., 2018).
The genome structure and encoded functions of the members of the genus Acetobacterium (Balch, Schoberth, Tanner, & Wolfe, 1977), are still not very well understood. The genes involved in the WL pathway of Acetobacterium woodi are divided into three clusters (Poehlein et al., 2012). Each of them consists of 6 to 10 syntenic genes, with their products orchestrating a specific phase of the WL pathway ( Figure 1). Cluster I consists of 7 genes encoding formate dehydrogenase and accessory enzymes catalyzing the reduction of carbon dioxide to formate. Cluster II contains 6 genes, underpinning the four steps leading from formate to acetyl-CoA. Cluster III encodes the enzymes involved in carbon fixation and production of acetate from acetyl-CoA (Poehlein et al., 2012). Here, we report new genome sequences of four acetogenic bacteria and perform a comparative genomic analysis focused on the gene clusters and protein networks of the WL pathway.
The three Acetobacterium species were grown in DSM 614 medium amended with fructose at a temperature of 22°C, while Alkalibaculum bacchi was grown in DSM 545 medium at a temperature of 37°C.

| DNA extraction, library preparation, and sequencing
Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue kit (Hilden, Germany), according to the manufacturer's protocol for gram-positive bacteria. Bacterial cells were harvested by centrifugation at 10,000g for 15 min and kept at 37°C for 1 hr with the enzymatic lysis buffer provided by the supplier. Cells were then placed at 56°C for 30 min and treated with RNase A. After column purification, DNA was eluted with 100 ml 10 mmol/L Tris/HCl, pH 8.0. Genomic DNA purity and integrity were assessed by measuring the absorbance at 260 nm (A260) and the ratio of the absorbance at 260 and 280 nm (A260/A280) with a NanoDrop ND-1000 spectrophotometer (Thermo Scientific). Genomic DNA concentration was measured by using the Qubit fluorometer (Thermo Fisher). Libraries were prepared using the Nextera XT DNA library preparation kit (Illumina, USA) with default settings, and sequenced on an Illumina MiSeq platform.

| Genome assembly and annotation
The quality of the reads was checked using the software fastqc (Andrews, 2010), and adaptor sequences were removed using trim_galore (Krueger, 2016). The assembly was performed with the software SPAdes version 3.8.0 (Bankevich et al., 2012), using all default parameters and the option "-careful." After assembly, contigs shorter than 500 bp and/or with a coverage below 3 were removed.
Pairwise Average Nucleotide Identity (ANI) values were calculated among the five sequenced genomes and the reference genome of A. woodii using the software pyani (Pritchard, Glover, Humphris, Elphinstone, & Toth, 2016). The output was visualized using the inhouse developed software DiMHepy, publicly available at https :// github.com/lucaT ribol i/DiMHepy. Genomes were annotated using Prokka (Seemann, 2014), using an ad hoc database created starting from the genome of A. woodii. Amino acidic sequences predicted by Prokka were used as input for EggNOG mapper for prediction of functional features (Huerta-Cepas et al., 2017). The outputs of Prokka were imported in R (R Core Team, 2012) for graphical depiction of genomic maps using the R-package GenoPlotR (Guy, Kultima, Andersson, & Quackenbush, 2011), based on the coordinates found by Prokka. To infer the number of shared genes among the five genomes we used Roary (Page et al., 2015), leaving all default settings beside the blastp identity parameter, that was set to 60 because the comparative analysis included a species from another genus (i.e., Alkalibaculum bacchi). Venn diagrams, based on presence/absence of homologous genes as inferred by Roary, were drawn using the web tool of the Bioinformatics and Evolutionary Genomics Department of the University of Gent (http://bioin forma tics.psb.ugent.be/webto ols/ Venn/).
To identify biosynthetic gene clusters for secondary metabolites, the genome sequences for each of the strains were uploaded in fasta format to the antibiotics and Secondary Metabolites Analysis SHell (antiSMASH) web server (Blin et al., 2017).

| Prediction of orthologues and paralogues
The protein sequences for the five species were predicted by Prokka, and all-versus-all sequence similarity searches between the protein set of each pair of the five considered species were performed independently using the BLASTp program of the BLAST package (Camacho et al., 2009). As proposed by Rosenfeld and DeSalle (2012), a paralogy analysis may consider an E-value threshold that maximizes the number of detectable protein families (Rosenfeld & DeSalle, 2012). Therefore, all similarity searches were initially carried out using an E-value cutoff of 10 −3 . In order to identify orthologues, we used a python software developed by Ambrosino et al. (2018). The software accepts the output of the BLAST similarity searches as input, implementing a Bidirectional Best Hit (BBH) approach (Hughes, 2005;Huynen & Bork, 1998;Overbeek, Fonstein, D'Souza, Pusch, & Maltsev, 1999;Tatusov, Koonin, & Lipman, 1997).
Such approach establishes that proteins a i and b i , from species A and B, respectively, are the best orthologues if a i is the best scored hit of b i , with b i being the best scored hit of a i , in all-versus-all BLAST similarity searches (Hughes, 2005). For paralogy prediction, all-versus-all similarity searches were performed for each species using the BLASTp program.

| Protein similarity networks
Networks of proteins based on the inferred similarity relationships were built. The network construction procedure extracted all the connected components into different separated undirected graphs by using NetworkX package (Hagberg, Schult, & Swart, 2008). Each node in the network represents a protein and each edge represents an orthology or paralogy relationship. A filtering step was introduced to select for each species only the E-value cutoff that maximized the number of paralogue networks. The selected E-values were e -10 for TA B L E 2 Genome annotation statistics, including number of CDS predicted by Prokka, antiSMASH gene clusters analysis and protein family annotation by eggNOG mapper (for A. woodii the analysis was done on the reference strain with acc.no. CP002987) for Alkalibaculum bacchi. Cytoscape software (Shannon et al., 2003) was used for the graphical visualization of the networks.

| Genome-wide analyses reveal close similarity between A. tundrae and A. paludosum
The number of reads per genome was on average 814.008 ± 251.751; the assembly resulted in an average number of contigs of 53 ± 9 (  Figure 2). It should be pointed out that A. bakii DSM 8239 was sequenced in another study (Hwang, Song, & Cho, 2015). We compared the previously sequenced genome of A. bakii with our data and found an ANI value of 99.76% over an alignment length of 4.12 Mb.
The ANI analysis confirms the evolutionary relationships between these species (Simankova et al., 2000), with A. paludosum and A. tundrae being most closely related within the genus Acetobacterium with an ANI of 95% over an alignment length of 6.4 Mbp. Alkalibaculum bacchi branched outside of the Acetobacterium F I G U R E 4 Organization of the three gene clusters in the four Acetobacterium genomes. Orthologues are connected with purple shades group, and displayed an ANI value of 70%, over an alignment length of 5.4 Mbp.
The annotation using Prokka found on average 3,343 ± 393 coding sequences. Proteins were assigned using EggNOG mapper to 2,460 ± 221 protein families ( Table 2).
The number of gene clusters involved in the production of secondary metabolites identified by the antiSMASH analysis was 12, 16, 15, and 18 in A. bacchi, A. bakii, A. paludosum, and A. tundrae, respectively (Table 2). A single cluster of genes for fatty acid biosynthesis per genome was found by the ClusterFinder algorithm, and this cluster was in all cases homologous to a cluster of 10 genes in The pangenome consisted of 9,262 genes, with a core genome of 363 genes (whose annotation is provided in Table A1), the number of core genes Acetobacterium spp. was 1,241. The number of unique genes into a single genome ranged from 486 to 2,360, in A. tundrae and A. bacchi, respectively (Figure 3).

| Gene cluster organization of the WL pathway is well conserved in Acetobacterium spp
As mentioned above, the WL pathway in A. woodii is encoded by three gene clusters. We examined the organization of those genes in three newly sequenced Acetobacterium species. The gene order was perfectly conserved (syntenic), compared with the reference strain Acetobacterium woodii, in the three clusters. A. bakii showed a truncated version of the formate dehydrogenase gene (fdhF1), whereas the other genes in this cluster were conserved (Figure 4). To confirm this observation, we searched the homologue of fdhF1 in the genome of A. bakii deposited in NCBI, which could not be identified. Consistently, a truncated version of fdhF1 in A. bakii was also found by Shin et al. (2018). In the genomes of A. tundrae and A. paludosum, the gene encoding formyl-tetrahydrofolate synthetase (fhs1, from cluster II), was duplicated ( Figure 4). One possible explanation for this feature could be the duplication of this specific gene as an adaptive trait. Examples of gene duplication are frequently connected to environmental adaptation (Tatusov et al., 1997), often through gene dosage (Bratlie et al., 2010;Kondrashov, 2012).

Gene cluster III presented no rearrangements in any of the four
Acetobacterium genomes (Figure 4). Conversely, in Alkalibaculum bacchi, genes of the WL pathway were organized in a different way compared to the Acetobacterium genus, as none of the three clusters was found to be complete. Genes appeared instead to be scattered all over the bacterial chromosome (Table A2). Only the formate dehydrogenase genes (and not the accessory proteins) of cluster I were found on two separate contigs. All genes of cluster II were found, although they were split between two contigs. All but two genes of cluster III were found on the same contig, although the gene order was not maintained (Table A2).

| Protein network analysis reveals gene expansion dynamics for WL pathway proteins
The comparative analysis performed on all considered species led to the construction of networks of protein orthologues and paralogues.
Prediction of orthologues between the five species was performed using a Bidirectional Best Hit (BBH) approach. Overall, 20,712 BBHs were detected. Paralogues were detected by all-against-all sequence similarity searches. Using as an input the predicted 20,712 orthology relationships, we considered the associated paralogues in all species, which led to the identification of a total of 2,135 distinct networks ( Figure 5). A general overview of the generated networks indicates that a consistent core of networks (922) contained proteins present in all considered species, while only 9, 21, 5, 7, and 48 networks contained proteins exclusively found in A. woodii, A. paludosum,A. tundrae,A. bakii,and A. bacchi,respectively (Figure 5).
We then inferred gene conservation or divergence between species pairs, calculating the number of proteins per species for each network ( Figure 6). We defined duplicated proteins starting exclusively from the previously detected orthologue pairs. Specifically, we defined 455 two-protein networks connected by a single orthology relationship, 1,424 networks including 3-9 proteins, and 256 networks containing 10 or more proteins (Figure 6a). The networks distributed along a hypothetical bisector (Figure 6b), which represent the protein families that did not undergo significant changes in the number of members between species pairs. In contrast, networks that are distant from the bisector represent expansions or reductions in the number of proteins of related protein families in A.
woodii compared to the other species. Furthermore, it is possible to infer the most conserved protein families between A. woodii and the other species by considering the networks with the highest number of orthologues (large circles in Figure 6).
We then selected the A. woodii proteins encoded by the genes of the WL pathway, identifying them within the generated networks. The proteins encoded by the gene clusters I, II, and III led to the discovery identification of 13 distinct networks ( Figure A1). At least one protein F I G U R E 7 Selected networks displaying different amplification patterns in genes involved in the Wood-Ljungdahl pathway. An extended version of this figure including all proteins of the WL pathway is presented in Figure A1 per cluster presented cliques of one orthologue per genome (Figure 7), this is the case for FdhD in cluster I, FolD in cluster II and AcsD in clus- bacchi proteins are completely missing ( Figure A1). This confirms the divergence highlighted in the previous comparative analyses.

| CON CLUS IONS
We obtained draft genome sequences for three Acetobacterium species and a acetogenic bacterium, Alkalibaculum bacchi. This study emphasizes the degree of genomic divergence and conservation of protein families within the genus. Having a closer look at the gene clusters involved in WL pathway, we revealed rearrangements and homology patterns that expands our understanding regarding the evolution of this metabolic pathway in the Acetobacterium genus with the perspective of future exploitation of these bacteria for industrial applications.

ACK N OWLED G M ENTS
The study was financed in part by the Autonomous Province of Trento (ENAM project) in cooperation with the Italian National Research Council (CNR). The authors thank Matthias Kirschberg for providing useful edits on the manuscript.

CO N FLI C T O F I NTE R E S T S
None declared.