Experimental evolution reveals nitrate tolerance mechanisms in Desulfovibrio vulgaris

Abstract

Elevated nitrate in the environment inhibits sulfate reduction by important microorganisms of sulfate-reducing bacteria (SRB). Several SRB may respire nitrate to survive under elevated nitrate, but how SRB that lack nitrate reductase survive to elevated nitrate remains elusive. To understand nitrate adaptation mechanisms, we evolved 12 populations of a model SRB (i.e., Desulfovibrio vulgaris Hildenborough, DvH) under elevated NaNO₃ for 1000 generations, analyzed growth and acquired mutations, and linked their genotypes with phenotypes. Nitrate-evolved (EN) populations significantly (p < 0.05) increased nitrate tolerance, and whole-genome resequencing identified 119 new mutations in 44 genes of 12 EN populations, among which six functional gene groups were discovered with high mutation frequencies at the population level. We observed a high frequency of nonsense or frameshift mutations in nitrosative stress response genes (NSR: DVU2543, DVU2547, and DVU2548), nitrogen regulatory protein C family genes (NRC: DVU2394-2396, DVU2402, and DVU2405), and nitrate cluster (DVU0246-0249 and DVU0251). Mutagenesis analysis confirmed that loss-of-functions of NRC and NSR increased nitrate tolerance. Also, functional gene groups involved in fatty acid synthesis, iron regulation, and two-component system (LytR/LytS) known to be responsive to multiple stresses, had a high frequency of missense mutations. Mutations in those gene groups could increase nitrate tolerance through regulating energy metabolism, barring entry of nitrate into cells, altering cell membrane characteristics, or conferring growth advantages at the stationary phase. This study advances our understanding of nitrate tolerance mechanisms and has important implications for linking genotypes with phenotypes in DvH.

Introduction

Sulfate-reducing bacteria (SRB) carry out dissimilatory sulfate reduction, and play an important role in the biogeochemical cycling of sulfur, carbon, nitrogen, and metals [1, 2]. They are ubiquitous in anaerobic environments such as marine and lacustrine sediments [3], where nitrate and sulfate are commonly found [4–6]. Nitrate concentrations in natural environments are generally <2 mM [7, 8], while in contaminated sites, nitrate concentrations can be greater than 100 mM [9]. Elevated nitrate could inhibit sulfate reduction performed by SRB due to its competition for organic electron donors [6], which has a special importance in contaminated or engineered systems since it mitigates sulfide production in the environment [10–12]. To confer an advantage under elevated nitrate, several SRB can use nitrate as an alternative electron acceptor [13]. However, the model SRB, Desulfovibrio vulgaris Hildenborough (DvH) lacks a functional nitrate reductase [14, 15], and how it responds to elevated nitrate and thrives in such environments for a long time exposure remains elusive [16]. Therefore, it is important to understand nitrate tolerance mechanisms, and we used DvH to advance our knowledge about SRB responses and functions in the environment in this study.

Previous studies showed complex responses of DvH to elevated nitrate by different approaches, such as iTRAQ peptide tags [17], microarrays [4], transposon liquid enrichment sequencing (TnLE-seq) [18], and mutagenesis [19]. For example, transcriptomic analysis of DvH showed that elevated nitrate up-regulated key functional genes involved in osmotic, ionic, and nitrite stress responses, and few common patterns of gene expression were observed in response to elevated nitrate, nitrite, or salt [4]. Further analyses indicated that a specific nitrate cluster (DVU0245-DVU0251) was related to increased nitrate tolerance, but not for nitrite tolerance in DvH [18], and elevated nitrate and nitrite had distinct effects on inhibition of DvH by mutagenesis [19]. Although these studies have provided new insights into our understanding of the response of DvH to short-term exposure to elevated nitrate at the protein, transcriptional and gene levels, more studies are needed to reveal the genetic basis of nitrate tolerance when DvH is exposed to elevated nitrate that mimics an environmental stress, especially for a long time.

Experimental evolution of asexual microbial populations coupled with whole-genome resequencing provides a powerful approach to identify adaptive mutations and elucidate genotype-phenotype relationships [20, 21]. Beneficial mutations typically emerge after an adaptive evolution of microorganisms [20, 22], which has been used to discover adaptive changes under specific conditions for prokaryotes such as Escherichia coli, Pseudomonas aeruginosa, and Desulfovibrio vulgaris [22–26], phages ϕ6 and T4 [27, 28], and eukaryotes such as Saccharomyces cerevisiae and Chlamydomonas reinhardtii [29, 30]. For example, a previous study indicated a rapid selective sweep of mutations occurred in D. vulgaris populations under NaCl stress within the first 100 generations, and four mutations contributed to shortened lag phases, increased growth rates and biomass yield [23]. Also, identical or similar genetic changes (e.g., mutations occurred in same/similar genes with similar phenotypic changes) in replicate populations are considered as a strong indicator of beneficial mutations resulting from stress adaptation in the experimental evolution [31–34]. Parallel evolutionary losses of ribose catabolic function occurred in 12 E. coli populations during 2000 generations in a glucose minimal medium to gain beneficial fitness in glucose-limited environments [35], and such genetic parallelism at the population level was also observed in the gene Dde_2265 encoding sulfate adenylythransferase (sat) in Desulfovibrio alaskensis evolved under perchlorate stress for about 117 generations, leading to increased perchlorate tolerance [36]. Therefore, it is expected that experimental evolution and resequencing of whole genomes provide powerful approaches to link genotypes to phenotypes, and explore nitrate tolerance mechanisms in DvH.

In this study, we aimed to understand nitrate tolerance mechanisms in DvH through experimental evolution, which could provide a comprehensive understanding of genetic basis of the adaptation evolution in DvH. We hypothesized that a set of beneficial genetic changes would emerge in DvH populations after exposure to elevated nitrate (NaNO₃) for 1000 generations [22, 23], and that some of these mutations would be parallel among evolved populations [20, 32, 34, 36]. We evolved 12 populations of DvH under elevated nitrate for 1000 generations, compared genotypic and phenotypic characteristics of ancestral (AN) populations and EN populations, linked their genotypes with phenotypes by mutagenesis, and explored its nitrate tolerance mechanisms. We found that six functional groups with a high frequency of mutations, which were directly or indirectly involved in nitrate tolerance. Most importantly, we showed new evidence that nonsense and frameshift mutations in key nitrate associated genes conferred nitrate tolerance and genetic parallelism in nitrate-evolved populations. This study advances our understanding of nitrate tolerance mechanisms in DvH.

Materials and methods

Strains, cultivation conditions, and propagation

We used 12 AN populations (AN1-12) of D. vulgaris Hildenborough (ATCC 29579) from a previous study [37], and a clone of DvH was the founder of 12 populations. The first subculture was non-stressed for all 12 populations, and they were frozen to serve as 12 AN populations, from which we randomly chose three AN populations (AN2, AN8, and AN11) to determine a suitable NaNO₃ concentration for experimental evolution, which allowed AN populations to reach a stationary phase within 48 h (Fig. S1). AN populations were grown at 37 °C and propagated in the defined LS4D medium [38] amended with 10 mM NaNO₃ for 500 generations, then with 30 mM NaNO₃ for another 500 generations to keep ongoing nitrate selective pressure for DvH populations. Cultures were propagated by transferring 1% (100 µL) of the final volume (10 mL) into a fresh medium every 48 h, and populations were archived every 100 generations at −80 °C for later studies. Every 48 h, the OD₆₀₀ was recorded for each population before transfer (Fig. S2), and the statistical significance of biomass differences (OD₆₀₀) of 6.7 generations (after the first transfer) vs. 1000 generations (after the last transfer) was based on ANOVA tests.

Growth phenotypes of ancestral and evolved populations

The growth of both AN and nitrate-evolved (EN) populations in the LS4D medium with 30 mM or 100 mM NaNO₃ or without additional NaNO₃ was measured by a spectrophotometer (Thermo Spectronic 20D+, Waltham, Massachusetts, USA) with three replicates for each population, and OD₆₀₀ values of all populations were recorded every 2–3 h over a period of 70 h. Lag phase was defined as the time after-inoculation to OD₆₀₀ = 0.2; growth rate was 2.303× the slope of the linear portion of growth curve by plotting Ln (OD₆₀₀) as the y-axis and time as the x-axis; biomass yield was the maximum OD₆₀₀ as previously described [37]. Statistical significance of multiple pairwise-comparison was based on the ANOVA test. In addition, extracellular nitrate concentration of AN and EN populations was measured by ion chromatography (ICS-90A, Dionex, USA).

Genomic DNA extraction and whole-genome resequencing of AN and EN populations

Genomic DNA (gDNA) of all AN and EN populations was extracted with a GenElute Bacterial Genomic DNA kit (Sigma) following the manufacturer’s instructions. The quality of gDNA was assessed with a Nanodrop ND-1000 (Thermo Fisher Scientific, Inc., Wilmington, DE, USA). Illumina HiSeq sequencing of gDNA samples was conducted at the Oklahoma Medical Research Foundation (Oklahoma City, OK, USA) following standard library preparation protocols and Illumine sequencing protocols. The sequencing data have been deposited into NCBI Sequence Read Archives and are accessible through Bioproject PRJNA630554.

Mutation calling

Mutations in all AN and EN populations were called by our in-house pipeline (http://zhoulab5.rccc.ou.edu:8080/) with the DvH genome sequence (NC_002937.3) as the reference. Btrim was used for quality control with filter parameter settings, including windows size as 5, average quality score as 30 and minimal insert size as 50 [39]. The alignment was generated by Bowtie2, and the default settings of bowtiw2-build and bowtie2-x were used for database building and alignments [40]. Samtools was used to call mutations including single nucleotide polymorphisms (SNPs) and insertions/deletions (INDELs); Samtools-mplieup was used to filter variants, probabilistic realignment for computation of base alignment quality was disabled to reduce false SNPs caused by misalignments, minimum mapping quality for alignment was 0, and minimum base quality for a base was 13 [41]. The output variant call format (VCF) file was generated by Bcftools [42]. Phred scale quality score was used as the filter score. Breseq was also used to call mutations with the default setting [43], and new junctions/insertions in output/index.html and output/marginal.html (allele frequency >10%) were added into the mutation data analysis. Compared to the original NCBI reference DvH strain, the ancestral DvH strain used in this study already has 22 genomic differences [23], which were removed from our mutation analysis. The allele frequency of mutations was calculated based on DP4 values: forward reference alleles, reverse reference alleles, forward alter alleles, and reverse alter alleles by Samtools [41].

Mutation validation

To validate mutations called in this study, a subset of regions bearing mutations with different filter values were chosen (Table S1). Genomic DNA from AN or EN populations was extracted and amplified by PCR, and PCR amplicons were sequenced by a Sanger platform at the Oklahoma Medical Research Foundation to verify those mutations identified by the Illumina sequencing platform. All the primers used in this study are listed in Table S1.

Mutation data analysis and prediction of mutated protein structure

Mutations were assigned to a gene if they occurred within a coding region or within a 100-bp upstream of open reading frames [44]. We combined output VCF files of whole-genome sequence data from Bcftools [42] to generate SNPs/INDELs tables for AN or EN populations (Tables S2 and S3). Briefly, mutation data analysis included five major steps: (i) all mutation information was extracted from output VCF files of AN and EN populations, including gene name, position for the mutation, mutated bases and reference bases; (ii) proteins encoded by these genes and Clusters of Orthologous Groups (COGs) of the proteins [45] were added into the table manually; (iii) as mutations or polymorphisms were detected in the ancestor [23], we only analyzed newly acquired mutations in the EN population (defined as newly acquired mutations), including new junctions occurred in coding genes that predicted by breseq (Table S4); (iv) to focus on mutated functional genes with newly acquired mutations in EN populations, a gene identified with newly acquired mutations was labeled with the highest allele frequency of the mutation occurring in the gene for each EN population (Table S5); and (v) gene co-expression prediction was based on information in operon and regulon prediction in Microbesonline (http://www.microbesonline.org/), and experimental evidence was based on publications about mutated genes, especially genes that were co-mentioned in same publication identified by String [46].

Types of detected mutations in this study include: (i) nonsense mutation: one-base change in a DNA sequence that results in a premature stop codon; (ii) frameshift mutation: an insertion or deletion in a DNA sequence that often results in a non-functional protein (e.g., truncated proteins) [47]; (iii) missense mutation: one-base change in a DNA sequence that results in a different amino acid in the mature protein; and (iv) silent mutation: one-base substitution in a DNA sequence that does not involve amino acid replacement in the coded protein.

We also used Phyre2 [48] to predict changes in the tertiary structure of HcpR encoded by DVU2547, one of the most frequently mutated genes in EN populations with the highest hits (>99% confidence) for predicted models (Fig. S3).

Mutagenesis and growth phenotypes of deletion mutants

DvH strain JZ001 (Δupp) was generated using the founder strain of AN populations, and this strain was used as the parental strain to generate marker exchange (ME) mutants of genes by a marker replacement approach [49]. The growth of parental strain and ME mutants of NSR was tested in LS4D, LS4D + 30 mM NaNO₃ and LS4D + 0.15 mM NaNO₂. The growth phenotypes of the parental strain and ME mutants of NRC genes were tested in LS4D, LS4D + 30 mM NaNO₃, and LS4D + 100 mM NaCl. Genetic complementation was done for all the mutants (ME2394, ME2395, ME2396, ME2405, ME2543, ME2547, and ME2548) and EN7 variant complemented with the plasmid containing a native promoter, ribosome binding site and wild type gene (Table 1). Further reverse-transcription PCR was done for the mutants and their complemented strains at the exponential phase to confirm the expression of complemented genes. The growth of complemented strains was tested in LS4D and LS4D + 30 mM NaNO₃. Statistical significance of growth parameters was based on the T-test. All DvH strains used in this study are listed in Table 1.

Results

Nitrate tolerance increased in nitrate-evolved populations

To choose a suitable nitrate concentration for our experimental evolution, three ancestral (AN) populations (AN2, AN8, and AN11) were randomly selected and grown under different nitrate concentrations. Their lag phase increased while growth rates and biomass gradually decreased as nitrate concentrations increased from 0 mM to 30 mM (Fig. S1). We evolved 12 populations under elevated nitrate (10 or 30 mM NaNO₃), which is higher than nitrate concentrations in natural environments for 1000 generations, and all 12 nitrate-evolved (EN) populations reached a higher biomass compared to AN populations (Fig. S2).

EN populations had a shorter lag phase (Fig. 1a) and a faster growth rate (Fig. 1b) compared to AN populations under normal growth conditions and showed significantly (p < 0.05) higher nitrate tolerance compared to AN populations (with 30 mM NaNO₃) (Fig. 1). With 100 mM additional NaNO₃, the growth of EN populations was repressed but they could still grow relatively well while AN populations barely grew within 70 h (Fig. 1). The nitrate concentration (30 mM) did not change at the end of the cultivation period for AN or EN populations, indicating that a change in nitrate concentration was not detected after DvH was evolved for 1000 generations under elevated nitrate.

Increased NaNO₃ tolerance in EN populations compared to AN populations.

Genes and functional gene groups with new mutations in EN populations

Sequencing reads of AN and EN populations covered about 99% of the DvH reference genome sequence (NC_002937.3) with an average sequencing depth of 88 × for AN populations, and 92× for EN populations. A total of 58 and 167 mutations were respectively detected in AN and EN populations, and the number of mutations ranged from 8 to 26 for each AN population (Table S2), while the number of mutations ranged from 13 to 33 for each EN population (Table S3). We focused on newly acquired mutations and associated genes in the EN populations, including 119 newly acquired mutations (containing nine large deletions and one insertion predicted using breseq) in 44 genes associated with energy, amino acid, carbohydrate, nucleotide, lipid, and inorganic ion metabolism, cellular process and signal, and information storage and processing (Table S4).

Among those 44 genes with newly acquired mutations, we focused on genes with high mutation frequency at the population level (>33%, refs to genes bearing mutations in more than four of 12 EN populations), and six functional gene groups of 19 genes were found based on gene co-expression prediction or experimental evidence from previous studies (Tables S4 and S5): (a) nitrosative stress response genes (NSR: DVU2543, DVU2547, and DVU2548): DVU2547 (hybrid cluster protein regulator, hcpR) harboring nine mutations in seven of 12 EN populations, regulates the expression of DVU2543 (hcp, identical mutation in two EN populations) [14, 50] and exists in the same operon with DVU2548 (acpD, five mutations in three EN populations) [51]; (b) nitrogen regulatory protein C family genes (NRC: DVU2394-2396, DVU2402, and DVU2405): a two component system DVU2394 (ntrC)/DVU2395 (response regulator/histidine kinase) in the operon DVU2394-DVU2396 (with 18 mutations in 11 EN populations) regulates the operon DVU2397-2405 [52], in which DVU2402 mutated in one EN population and DVU2405 harboring three mutations in two EN populations; (c) nitrate cluster (DVU0246-0249 and DVU0251) [18] with 17 mutations in all 12 EN populations; (d) fatty acid synthesis genes (DVU1204 and DVU1208): DVU1204 (3-oxoacyl-(acyl-carrier-protein (ACP) synthase II, fabF) harboring six mutations in seven EN populations and DVU1208 (phosphate acyltransferase, plsX) harboring an identical mutation in two EN populations are predicted in one operon involved in fatty acid synthesis [23, 51]; (e) iron regulatory genes (DVU0942 and DVU2571): DVU0942 (fur) harboring four mutations in six EN populations is a regulator of ferrous iron transport protein B gene DVU2571 (feoB) [53], which mutated in the other six EN populations; and (f) a two component system DVU0596/DVU0597 (lytR/lytS) harboring five mutations in six EN populations regulates the expression of carbon starvation protein genes DVU0598 (cstA) and DVU0599 (cstB) [52] (Table 2).

A large proportion (84/119) of mutations occurred in those 19 genes of six functional groups in EN populations (Tables 2 and S4). For example, seven EN populations had mutations in DVU2547, including an identical mutation occurred in EN7 and EN8; nine EN populations had mutations in the two component system DVU2394/DVU2395; all 12 EN populations harbored mutations in nitrate cluster. For the other 27 mutated genes with 35 newly acquired mutations, we found seven genes with 17 new mutations acquired in at least two EN populations (Table S6), and 18 genes acquired new mutations in only one EN population (Table S5). In the following analysis, we focused on those highly mutated six functional groups (with mutation frequency ≥ 50% at the populations level) to link their genotypes with phenotypes towards understanding nitrate tolerance mechanisms in DvH.

Six functional gene groups were highly mutated in EN populations

We further analyzed potential effects of mutations on functions of their associated genes/proteins in the six functional groups. First, nonsense and frameshift mutations were identified in the first three functional gene groups, presumably resulting in premature or truncated proteins. (i) Nonsense and frameshift mutations were identified in the NSR group of three genes with allele frequencies ranged from 17 to 100% (an average of 72%). Seven of 17 mutations in this group were fixed or nearly fixed (>90%), and interestingly, five of them were nonsense and frameshift mutations. Also, the same mutation was identified in this group, such as a nonsense mutation in coordinate 2654808 in DVU2543 (Hybrid cluster protein, hcp) in EN4 and EN5, and in position 2660403 of DVU2547 (Hcp regulator, hcpR) in EN7 and EN8 (Table 3). (ii) Among 22 mutations that occurred in the NRC group, 11 nonsense and frameshift mutations were identified in this group with four of them fixed in all EN populations, and the allele frequency of mutations ranged from 11 to 100% with an average of 58% (Table 4). (iii) Nonsense and frameshift mutations were also observed in three genes (DVU0246, DVU0247, and DVU0248) of nitrate cluster with a mutation allele frequency ranged from 13 to 77% (an average of 36%). All the mutations occurred in DVU0246 (Pyruvate phosphate dikinase, ppdk) and DVU0247 (Response regulator, cheY like/RR) could result in truncated proteins. Identical frameshift mutation was identified at the coordinate 279997 of DVU0246 in EN1 and EN8, and the same missense mutation occurred at the coordinate 285948 of DVU0251 (Transmembrane protein, tauE like) in EN 7 and EN 12 (Table 5). Interestingly, eight of ten mutations detected using breseq were located in these three functional groups (NSR, NRC, and nitrate cluster), which were nonsense and frameshift mutations (Table S4).

Second, missense mutations occurred in the other three gene groups, potentially resulting in functional changes in their encoding proteins. Mutations in fatty acid synthesis and iron regulatory genes had higher allele frequencies (43–100% and 88–100% with an average of 92% and 99%, respectively) compared to other gene groups (Table S4). For fatty acid synthesis genes, mutations in DVU1204 (fabF) were observed in seven out of 12 EN populations, which were fixed or nearly fixed, and the same SNP was detected in DVU1204 in EN4 and EN5, while a same silent mutation was found in DVU1208 (plsX) in EN8 and EN10 (Table S7). Also, all EN populations harbored mutations in iron regulatory genes, and the same mutation occurred in DVU0942 (fur) in EN4, EN5, and EN9 (Table S8). In addition, six mutations were detected in lytR/lytS with the lowest allele frequency (11–40% with an average of 19%) among those six gene groups, but all of the mutations were missense except an intergenic mutation at the upstream of DVU0597 start codon in EN3 (Table S9).

The results showed that nonsense, frameshift and missense mutations occurred in those six gene groups, and they may result in prematurely terminated proteins, truncated proteins, or function-changed proteins, thus affecting functions of DvH in response to elevated nitrate. Furthermore, as those mutations appeared to parallelly occur in EN populations, it is expected that such observed mutations may be related to increased nitrate tolerance in DvH.

Deletion mutants of NSR and NRC genes resulted in improved nitrate tolerance

To determine if nonsense and frameshift mutations are involved in nitrate tolerance in EN populations and to understand possible nitrate tolerance mechanisms in DvH, we generated representative deletion mutants of NRC and NSR genes, and analyzed their growth under elevated nitrate, NaCl, or nitrite. The NSR group includes three genes: hcp (hybrid cluster protein, DVU2543), hcpR (hcp regulator, DVU2547) and acpD (DVU2548) (Table 2). DVU2547 encoding HcpR was one of the most frequently mutated genes in seven out of 12 EN populations (Table 3). Also, these mutations were identified in key domains of hcpR (DVU2547), including DNA-binding domain (EN1-1, EN1-2), cAMP-binding domain (EN1-2), and dimer interface (EN3-1, EN3-2, EN6, EN9) (Fig. S3), and the same nonsense SNP was fixed in EN7 and nearly fixed in EN8. Analysis of those mutants indicated that the growth of parental strain JZ001 and NSR gene deletion mutants was quite similar without nitrate or nitrite addition (Fig. 2a), but these mutants grew faster than JZ001 with 0.15 mM nitrite added although a lower final biomass was obtained for the mutants (Fig. 2b). Especially, when 30 mM was added, all three mutants had shorter lag phases, higher growth rates, and higher biomass yield than JZ001 (p < 0.05), indicating increased nitrate tolerance in those mutants (Fig. 2c). In addition, a similar growth was observed for JZ001P (JZ001 with empty vector pMO719) and the complemented strains (ME2543 and ME2547) without nitrate (data not shown), while decreased nitrate tolerance was observed in the complemented strains (Fig. 2d) and EN7 variant (Fig. S4) under elevated nitrate, possibly due to the overexpression of complemented genes. It is noted that we could not obtain a complemented strain for ME2548 after several attempts. The results demonstrated that nonsense and frameshift mutations in the NSR group (Table 3) could improve nitrate tolerance in EN populations.

Growth phenotype of DvH parental strain (JZ001) and marker exchange mutants (ME2543, ME2547 and ME2548) of nitrosative stress response genes (NSR).

NRC group contains a two component system DVU2394 (Response Regulator, RR)/DVU2395 (Histidine Kinase, HK), which regulates the expression of operon DVU2397-2405 encoding alcohol dehydrogenase and heterodisulfide reductase (Table 2). Marker exchange (ME) mutants of DVU2394, DVU2395, DVU2396, and DVU2405 in the NRC group were generated. While no significant growth differences were observed between ME mutants and JZ001 in LS4D (Fig. 3a) or LS4D with 100 mM NaCl (Fig. 3b), all NRC mutants had higher growth rates (p < 0.05) in LS4D amended with 30 mM nitrate except ME2405, which had the shortest lag phase compared to both JZ001 and other ME mutants (Fig. 3c). As expected, the parental strain (with empty plasmid pMO719)-like phenotypes were observed for all complemented strains of NRC mutants (Fig. 3d). The results confirmed that nonsense and frameshift mutations in the NRC gene group (Table 4) likely led to increased nitrate tolerance in EN populations.

Growth phenotype of DvH parental strain (JZ001) and marker exchange mutants (ME2394, ME2395, ME2396, and ME2405) of nitrogen regulatory protein family C genes (NRC).

Discussion

Experimental evolution coupled with whole-genome resequencing enable us to link genotype with phenotype, understand nitrate tolerance mechanisms, and facilitate the advancement of our knowledge about SRB responses to elevated nitrate in the environment. In this study, we found that most newly acquired mutations occurred in six gene groups of DvH through experimental evolution under elevated nitrate, and deletion mutants of NSR and NRC genes confirmed their roles in nitrate tolerance, indicating that those nonsense and frameshift mutations were beneficial. Also, identical or similar mutations were observed in EN populations evolved independently, suggesting genetic parallelism of adaptive selection at the population level.

Beneficial mutations in stress-evolved populations are considered as drivers of adaptation [21, 32]. Previous studies of adaptive evolution in Desulfovibrio species showed beneficial mutations emerged along the experimental evolution and conferred a fitness advantage under stress conditions [23, 36]. In this study, we found a set of high-frequency mutations in six gene groups of EN populations, and they conferred nitrate tolerance in DvH, appearing to be beneficial. However, native functions of those gene groups are not likely involved with nitrate metabolism but largely by gene regulations, such as three two-component systems (DVU2394/DVU2395, DVU0596/DVU0597, and DVU0246/DVU0247) and two global transcriptional regulators/regulons (DVU2547/DVU2543 and DVU0942/DVU2571) [14, 18, 52–54]. Those regulatory systems may perceive nitrate stress and control downstream responses, and based on functions of mutated genes and the fact that no nitrate was consumed, our results suggest that adaptive mechanisms of DvH under elevated nitrate may include (i) shifting of energy metabolism, (ii) barring entry of nitrate into the cell, (iii) altering cell membrane characteristics, and (iv) conferring growth advantages at the stationary (GASP) phase.

First, a shift of energy metabolism could provide physiological advantages for DvH to survive nitrate stress [4], thus increasing nitrate tolerance after experimental evolution. A previous study indicated that a deletion mutant of DVU0916 (rex), a repressor of sulfate adenylyl transferase increased nitrate tolerance [18]. We found that a large number of mutations were identified in NSR and NRC genes involved in energy metabolism, and functional loss of those genes conferred nitrate tolerance. DVU2547 (hcpR) of NSR was proposed to be involved in energy metabolism through repressing sulfate reduction pathways and activating expression of DVU2543 (hcp), which responds to reactive nitrogen species and was up-regulated under elevated nitrite or nitrate [4, 55–57]. A possible explanation for growth advantages of an hcpR mutation was suggested to be derepression of sulfate-reducing genes resulting from the loss of the repressor HcpR (Fig. 4a). Also, we found NRC genes conferred nitrate tolerance, which has not been reported before. In this gene group, DVU2405 (alcohol dehydrogenase, adh) regulated by a two-component system (DVU2394/DVU2395) was shown to be one of the most highly expressed genes contributing to energy metabolism with different electron donors in DvH [58]. Thus functional loss of genes in the NRC group may lead to a shift to methyl/S-adenosyl-methionine metabolism due to a lack of alcohol oxidation pathways catalyzed by alcohol dehydrogenase encoded by DVU2405 [4, 58] (Fig. 4b). Our results suggested that those frameshift and nonsense mutations in the NSR and NRC gene groups could relieve nitrate stress through shifting energy metabolisms, and those new findings advance our understanding the role of NSR and NRC groups in nitrate tolerance in DvH.

A schematic representation of possible nitrate tolerance mechanisms in DvH.

Second, blocking nitrate entry into the cell could confer nitrate tolerance in DvH. The nitrate cluster genes were proposed to allow non-specific nitrate transport, while transposon mutants of nitrate cluster genes were found to have a growth advantage under nitrate stress [18]. Also, DVU0249 encodes an outer membrane-associated homodimer [59], and DVU0251 is a TauE like transmembrane protein [18]. Therefore, we speculate that those nitrate cluster genes may be involved in nitrate uptake or transport, and the functional loss of those genes may prevent nitrate entry into the DvH cell [18], which is consistent with our study, showing that nonsense and frameshift mutations occurred in nitrate cluster (Fig. 4c).

Third, altering membrane lipid composition is crucial for bacterial survival and adaptation under environmental stresses [60]. Membrane fatty acid synthesis was shown to be involved in heat stress, salt stress, and oxidative stress adaptation in E.coli [61]. Previous studies also indicated that the percentage of unsaturated fatty acid in cell membrane was affected under salt stress, while the formation of acyl-accepter was catalyzed by DVU1204 (fabF), a possible essential gene in DvH [23, 62], and it was up-regulated under heat stress [63]. We found that the gene fabF was high-frequently mutated in EN populations, which could increase unsaturated fatty acid percentages in cell membrane, thus increasing nitrate tolerance possibly by relieving osmotic stress (Fig. 4d).

Fourth, GASP could be beneficial for bacteria to adapt to different stresses due to competitive ability conferred by GASP mutations (e.g., in the gene rpoS) emerged during a long-term stationary phase [64]. A previous study with E. coli indicated the transcriptional regulator RpoS exhibited GASP mutations and responded to multiple stresses [65]. However, DvH does not have an annotated rpoS gene. Alternatively, iron homeostasis maintained by DVU2571 (feoB) and up-regulation of DVU0598 and DVU0599 encoding carbon starvation proteins regulated by lytR/lytS, which suggested being involved in stresses response as DvH cells transitioned into the stationary phase [66]. Furthermore, DVU2571 and lytR/lytS were also mutated in serial transfers with or without salt stress [23]. Such an experimental evolution may emerge mutations in those genes up-regulated at the stationary phase, as evolved DvH populations were kept at the stationary phase for ~24 h in every transfer (48 h). In this case, it is reasonable to assume that newly acquired mutations in those GASP-related genes may increase nitrate tolerance through relieving possible environmental stresses like nutrient depletion (Fig. 4e, f).

Notably, some beneficial mutations occurred in evolved populations may just increase the overall growth (e.g., shortened lag phase, increased growth rate, increased biomass) without direct effects for a specific stress during an experimental evolution [67]. For example, our previous study showed that a mutation in DVU0597 was related to increased growth and biomass as well as shortened lag phase while another mutation in DVU2472 was related to increased biomass and shortened lag phase in salt-evolved populations of DvH [23]. In this study, we found shortened lag phase and increased growth rate for EN populations with or without nitrate stress, which is similar to a previous study, showing that both control-evolved populations (EC) and salt-evolved populations (ES) had a shorter lag phase with or without salt stress compared to their ancestral populations (without evolution) in DvH [37]. Also, a few shared mutated genes (DVU1204, DVU2571, and DVU0596/DVU0597) were found in those populations (EN, EC, and ES), which may confer an increase of general stress tolerance [23, 60, 62–64, 66]. Furthermore, the fact that high concentrations of iron could be regulated by DVU2571 during the lag phase emphasized the role of iron accumulation in the early stage of bacterial growth [68]. Therefore, our results highlighted potential roles of general stress response genes, especially GASP associated genes (DVU0597 and DVU2571) in promoting the overall growth in EN populations.

In addition, parallel genetic changes were observed at the mutation, gene, and gene group levels in independent populations. At the mutation level (identical mutations among EN populations), identical mutations were observed in six gene groups and other less defined genes, and they were nearly fixed or fixed. For example, identical nonsense and frameshift mutations occurred in the NSR gene or the nitrate cluster of independent EN populations, resulting in truncated proteins and increased nitrate resistance. Identical fixed or nearly fixed missense mutations were also observed at the coordinate 1034711 in DVU0924 in EN4, EN5, and EN9 (Tables S3 and S8). At the gene level (different mutations in the same gene across EN populations), nonsense and frameshift mutations were observed for genes in the NSR, NRC, and nitrate cluster, which may result in nonfunctional proteins and similar phenotypic changes that conferred nitrate tolerance. Moreover, parallel genetic changes were also reflected at the gene group level (genetic changes occurred in functional gene groups of independent populations), leading to similar phenotypic shifts, especially in NSR, NRC, and nitrate cluster. For example in the NRC group, loss-of-function mutations were identified in DVU2394 (RR) in EN5 and EN9, in DVU2395 (HK) in EN5, EN9, EN10, and EN12, in DVU2396 in EN1, EN2, and EN6, indicating that these genetic changes were similar and contributed to nitrate tolerance in independent populations. Notably, in this study, parallel evolution was more likely to happen via loss-of-function mutations. Therefore, such genetic parallelism involved in SNP/INDELs coupled with their high allele frequencies under elevated nitrate may be largely beneficial mutations in EN populations.

In summary, this study reveals nitrate tolerance mechanisms largely associated with six gene groups by experimental evolution in DvH. When DvH was exposed to elevated nitrate for 1000 generations, a number of beneficial mutations occurred. Nonsense and frameshift mutations were identified in NSR, NRC and nitrate clusters, resulting in increased nitrate tolerance through regulating energy metabolism and barring entry of nitrate into DvH cells, while missense mutations in fatty acid synthesis genes, iron regulatory genes, and lytR/lytS may confer general stress tolerance by experimental evolution under elevated nitrate through altering cell membrane characteristics and conferring GASP. Also, genetic parallelism is reflected at the mutation, gene, and gene group levels, allowing us to further understand DvH responses to elevated nitrate. This study unravels genetic basis of evolutionary changes in nitrate-evolved populations and provides a comprehensive understanding of nitrate tolerance mechanisms, which has important implications for linking genotypes with phenotypes in DvH. Future studies may further understand possible indirect nitrate tolerance mechanisms using a coinciding evolved control without additional nitrate and explore the contribution of important individual mutations and other highly mutated genes to nitrate tolerance in EN populations.

Acknowledgements

The experimental work of this study (October 2013 to June 2017) was largely supported by ENIGMA-Ecosystems and Networks Integrated with Genes and Molecular Assemblies (http://enigma.lbl.gov), a Science Focus Area Program at Lawrence Berkeley National Laboratory is based upon work supported by the U.S. Department of Energy, Office of Science, and Office of Biological & Environmental Research under contract number DE-AC02-05CH11231, the Office of the Vice President for Research at the University of Oklahoma (JZ), and the National Program on Key Basic Research Project of China (973 Program, No. 2015CB150505 to SC); the subsequent integration work (October 2018 to March 2020) of this study, including supplementary experiments, data analysis and manuscript writing, was mainly supported by the National Natural Science Foundation of China (Grant numbers 31770539 and 91951207 to ZH, grant number 31672262 to QY), and the Special Funds for Scientific and Technological Innovation in Guangdong Province (Grant number 2018A030310302 to BW).

Compliance with ethical standards

Conflict of interest

The authors declare that they have no conflict of interest.

References

Supplementary information The online version of this article (https://doi.org/10.1038/s41396-020-00753-5) contains supplementary material, which is available to authorized users.