Genetic Analysis of Benzothiophene Biodesulfurization Pathway of Gordonia terrae Strain C-6

Sulfur can be removed from benzothiophene (BT) by some bacteria without breaking carbon-carbon bonds. However, a clear mechanism for BT desulfurization and its genetic components have not been reported in literatures so far. In this study, we used comparative transcriptomics to study differential expression of genes in Gordonia terrae C-6 cultured with BT or sodium sulfate as the sole source of sulfur. We found that 135 genes were up-regulated with BT relative to sodium sulfate as the sole sulfur source. Many of these genes encode flavin-dependent monooxygenases, alkane sulfonate monooxygenases and desulfinase, which perform similar functions to those involved in the 4S pathway of dibenzothiophene (DBT) biodesulfurization. Three of the genes were found to be located in the same operon, designated bdsABC. Cell extracts of pET28a-bdsABC transfected E. coli Rosetta (DE3) converted BT to a phenolic compound, identified as o-hydroxystyrene. These results advance our understanding of enzymes involved in the BT biodesulfurization pathway.


Introduction
Sulfur-oxide gases (mostly SO 2 ) are released during the combustion of sulfur compounds, which results not only in serious air pollution, but also poison metal catalysts [1]. Benzothiophene (BT), dibenzothiophene (DBT), and their alkylated derivatives, account for more than 50% of the total sulfur content of commercial diesel [2,3]. These aromatic thiophenes are recalcitrant organic sulfur compounds and more resistant to hydrodesulfurization (HDS) treatment than other sulfur compounds such as mercaptans and sulfides [4]. Fortunately, several bacterial species can efficiently desulfurize aromatic thiophenes such as BT, DBT and their derivatives under ambient temperature and pressure [5].
As DBT is a common organosulfur compound found in a variety of fuels and more resistant to HDS treatment than other thiophene sulfides, it is widely regarded as a model compound for the isolation and characterization of biodesulfurization bacteria [6]. Two major pathways of DBT biodesulfurization have been reported and designated the "Kodama" and "4S" pathways, respectively. In the former, DBT is converted to hydrophilic organosulfur compounds by a series of oxidations of one of the aromatic rings of DBT [7]. The latter is a nondestructive desulfurization pathway and has been widely studied. In the "4S" pathway, DBT is successively converted to sulfoxide (DBTO), sulfone (DBTO 2 ), sulfinate (HPBSi) and hydroxybiphenyl (HBP) without degradation of either aromatic ring [8]. These four steps catalytic reactions are performed by three enzymes (DszC, DszA, and DszB), which are encoded by a single operon (dszABC). DBT is initially oxidized by DszC, first to DBT-5-oxide (DBTO) and then to DBT-5, 5'-oxide (DBTO 2 ). DszA catalyzes the transformation of DBTO 2 to 2-(2'hydroxyphenyl) benzene sulfinate (HPBSi), which opens the thiophenic ring. HPBSi is then desulfinated by DszB to produce 2-hydroxybiphenyl (2-HBP) [9]. These enzymes have been isolated, cloned, mutated, overexpressed, and crystallized [10][11][12][13][14]. It was found that DszC and DszA are flavin-dependent monooxygenases that require FMNH 2 or FADH 2 for their activities. Therefore, an additional enzyme, flavin reductase, is required, which catalyzes the reduction of flavins, such as FMN and FAD, by NAD(P)H to form reduced flavins.
The majority of DBT biodesulfurization bacteria, such as Rhodococcus sp. IGTS8, Rhodococcus erythropolis SHT87 [15], Mycobacterium sp. ZD-19 [16], Pseudomonas stutzeri TCE3 [17] and Rhodococcus erythropolis LSSE8-1 [18] seem to be incapable of BT desulphurization despite the similarity in their structures. Only a few bacteria harboring dszABC or their homologs (Table S1 in File S1), such as Paenibacillus sp. strain A11-2 [19,20], Mycobacterium goodii X7B [21], Gordonia alkanivorans RIPI90A [22,23] are capable of the biodesulfurization of both DBT and BT. However, the participation of DszA, DszB and DszC in the metabolism of BT has not yet been confirmed. Gordonia sp. 213E (NCIMB 40816) [24] and Sinorhizobium sp. KT55 [25] are specific BT biodesulfurization bacteria, which are capable of extracting sulfur from BT without breaking carbon-carbon bonds. Based on the characterization of metabolites produced during the BT biodesulfurization process by strain 213E and strain KT55, two types of BT biodesulfurization pathways have been proposed ( Figure 1 [25]. Metabolites of BTO 2 and BcOTO were detected in these desulfurization pathways although their final desulfurization products are different. Some bacteria that desulfurize BT via these two routes have been isolated and characterized, such as Paenibacillus sp. A11-2, Mycobacterium goodii X7B could desulfurize BT to yield ohydroxystyrene [20,21]. These previous studies have improved our understanding of the BT biodesulfurization mechanism, although the biochemical and genetic bases of the BT biodesulfurization pathway remain to be elucidated. G. terrae strain C-6 desulfurizes BT to yield a phenolic compound as the final product. It is hypothesized that this strain contains enzymes that are similar in function to those involved in DBT biodesulfurization, for that the biodesulfurization pathway and metabolites are so analogous to those of DBT biodesulfurization. Based on this hypothesis, we compared the transcriptomic profiles of strain C-6 cultured with BT or sodium sulfate as the sole source of sulfur. Among the significantly up-regulated genes, the bds (BT biodesulfurization) operon involved in the BT biodesulfurization pathway was identified. This information improves our knowledge of the mechanisms of aromatic thiophene metabolism.

Results and Discussion
Profiling of cell growth, BT utilization, and product yield of strain C-6 cultured with BT as sole source of sulfur When strain C-6 was cultured in the BSM S-medium containing BT as the sole sulfur source, maximum growth was achieved at 48 h, and the turbidity at 600nm at this time point was 1.36. No growth was observed without the addition of an alternative carbon source. These results indicate that this strain utilized BT as a sulfur source but not as sole carbon source. The Gibbs' assay showed that a phenolic compound was accumulated in the medium during the culture process, indicating that strain C-6 desulfurized BT to yield HPEal or ohydroxystyrene according to the reported pathway [24,25]. The maximum phenolic compound formation, close to 50% of the initial BT concentration, was obtained after cultivation for 48 h. This indicated that the consumption of BT was more rapid than the production of the phenolic compound during the culture process ( Figure 2). It can be speculated that this phenomenon is caused by the sublimation of BT, which is a sublimate organic compound at room temperature, in the culture process. This hypothesis was confirmed by loss of BT in a control experiment without inoculating strain C-6 in the culture ( Figure  2).

Global gene expression patterns of strain C-6 cultured with different sources of sulfur
Most organisms cannot convert sulfur directly into biomass, especially the sulfur contained in heterocycles such as BT and DBT. It must be converted initially into SO x 2-by a series of enzymes involved in the aromatic thiophene biodesulfurization pathway, and then converted subsequently into biomass by another group enzymes involved in sulfur cycle, such as sulfite reductase [26,27]. In order to clarify the BT desulfurization mechanism of strain C-6, we investigated the differential expression of genes in strain C-6 cultured with BT or sodium sulfate as the sole source of sulfur. In this case, the upregulated genes would be mainly related to sulfate limitation and/or induced by BT as part of BT desulfurization. To obtain an overview of strain C-6 gene expression profiles associated with growth on different sources of sulfur, cDNA samples were prepared and sequenced on anIllumina HiSeq 2000 platform. A total of 2,158,985,178 and 2,239,588,440 bases from 13,327,069 and 13,824,620 read pairs with a mean read length 81 bp were obtained from bacterial samples cultured with BT and sodium sulfate, respectively. These raw data were assembled into 3,157 contigs, and 4,800 unigenes were finally generated. Among these, 3,182 proteins had homologs in the COG database. With the exception of R, S categories of COG which are "function unknown" or "general function prediction only", others were related to the normal physiological metabolism of cells, such as regulation, transport and cell processing ( Figure S1 in File S1). In comparison with sodium sulfate as the sulfur source, the expression levels of 151 genes were deferent when BT was used as the sulfur source. Among these genes, 135 were up-regulated and 16 were downregulated ( Figure 3). Most of the significantly up-regulated genes (fold change more than 20) were shown to encode FMNH 2 -dependent monooxygenase, alkanesulfonate monooxygenase, desulfinase and transport related proteins, which are similar in function to enzymes involved DBT biodesulfurization of Rhodococcus sp. strain IGTS8 [8,9]. On this basis, 37 possible BT biodesulfurization related genes were selected and further investigated by RT-qPCR. All these genes were reconfirmed to be up-regulated in accordance with the RNA-seq result.

In silico analysis of BT desulfurization pathway genes
Both types of proposed BT biodesulfurization pathways are initiated by the activity of a monooxygenase, which catalyzes the reaction BT →BTO →BTO 2 ( Figure 1) [20,24,25]. Five of the significantly up-regulated genes (Table 1) are predicted to encode FMNH2-dependent monooxygenases, suggesting that the catalytic mechanism is analogous to that of DBT monooxygenase (DszC) [28]. Five additional genes were predicted to encode alkanesulfonate monooxygenases, which contribute to the conversion of BTO 2 to benzo[e] [1,2]oxathiin S,S-dioxide (Table 1). A desulfinase encoding gene was also identified among the significantly up-regulated genes. Its product catalyzes the last step of the 4S desulfurization pathway, which involves the hydrolysis of sulfinate or sulfonate from benzene sulfinate or benzene sulfonate [29]. In addition to these genes proposed to be involved in the 4S desulfurization pathway, NAD(P)H-dependent flavin reductase encoding gene was also identified among the significantly up-regulated genes. It provides flavin-dependent monooxygenase with reduced flavin [28]. As a general rule, genes involved in the same metabolic pathway usually exist in the form of an operon on the chromosome or plasmid. Based on this hypothesis, we sequenced the genome of strain C-6 [30]. The genetic map of chromosome reveals an operon, designated bdsABC (accession number:KC831580), containing only a desulfinase gene (Gene ID GTC6_01605, designated bdsA), a FMNH 2dependent monooxygenase gene (Gene ID GTC6_01600, designated bdsB) and an alkanesulfonate monooxygenase gene (Gene ID GTC6_01595, designated bdsC), which encode enzymes functionally analogous to DszB, DszC and DszA respectively. Other up-regulated genes are dispersed throughout the chromosome of strain C-6. As the most likely operon involved in BT desulfurization, the bdsABC operon was selected for the subsequential analysis of desulfurization activity although enzymes encoded by this operon show low homology with that involved in the "4S" desulfurization pathway of DBT (Table S1 in File S1). Another group of significantly upregulated genes was found to be related to the transmembrane transport of molecules, which maybe involved in the transport of BT and its metabolites [29].

Involvement of the bdsABC operon in BT biodesulfurization
To verify the putative roles of the bdsABC operon in the BT biodesulfurization pathway, we cloned and overexpressed this operon in E. coli Rosetta (DE3) with the expression vector pET28a (+). As expected, a phenolic compound was accumulated (indicated by the positive result in the Gibbs′ assay) when BT was incubated with cell extracts of E. coli Rosetta (DE3) overexpressing pET28a-bdsABC. NADH, FMN and dszD, a NAD(P)H-dependent flavin reductase involved in the DBT biodesulfurization pathway [28], were essential for the activity of these cell extracts, indicating that some of these enzymes involved in the BT biodesulfurization pathway are also flavin-dependent. In order to identify the structure of this phenolic compound and intermediate metabolites, the reaction mixture was extracted directly with ethyl acetate after incubation for 24h at 30°C. Approximately 95% volume of the ethyl acetate extracts were volatilized at room temperature for concentrating BT biodesulfurization metabolites. The remaining extracts were analyzed by GC-MS. Compared with the control experiment in which cell extracts of E. coli Rosetta (DE3)-  pET28a-bdsABC were replaced with that of E. coli Rosetta (DE3), no other different peaks were identified by the GC profiles, with the exception of the peak with retention time of 3.67 min. MS spectra of this peak are identical to that of ohydroxystyrene ( Figure S2 in File S1). Nonetheless, the BT biodesulfurization pathway, which is analogous to the 4S pathway of DBT, is becoming clearer based on the genetic analysis of BT metabolism, the final metabolite of BT and the results reported previously (Figure 4) [20,24,25].

Conclusion
G. terrae strain C-6 desulfurizes BT to yield ohydroxystyrene by mechanism analogous to the 4S biodesulfurization pathway. Transcriptomic profiles of strain C-6 cultured with BT or sodium sulfate as the sole source of sulfur indicated that three genes were involved in the BT biodesulfurization. These genes were located in the same operon. Cell extracts of E. coli Rosetta (DE3) overexpressing pET28a-bdsABC also converted BT to o-hydroxystyrene, which supported the hypothesis that the bdsABC operon is involved in BT biodesulfurization.

Chemicals
BT, BTO 2 and FMN were purchased from Sigma-Aldrich (Missouri, USA). Kanamycin and NADH were purchased from Amresco (Ohio, USA). Isopropyl-β-D-1-thiogalactopyranoside (IPTG) was purchased from Merck (New Jersey, USA). Polymerase chain reaction (PCR) primers were synthetized by Sangon (Shanghai, China). Pfu DNA polymerase was obtained from Promega (Wisconsin, USA). T4 DNA ligase and restriction endonucleases were obtained from Fermantas (Maryland, USA). All other reagents were of chromatography grade and were obtained from various commercial sources.

Substrate and product analysis
Quantitative analysis of BT was performed by HPLC (Agilent 1100 series, California, USA) fitted with a ZORBAX SB-C18 column (4.6 mm i.d. × 250 mm length). The mobile phase consisting of 80% methanol and 20% water was pumped at a flow rate of 1 ml/min. The absorbance of BT was measured at 245 nm. Final biodesulfurization product of BT was quantified by Gibbs′ assay according to previously described methods [31]. Qualitative analysis of BT biodesulfurization products was performed according to previously described methods [25]. Before analysis, BT and its metabolites in the aqueous phase were extracted with an equal volume of ethyl acetate. The ethyl acetate layer was centrifuged (1,600 ×g, 10min), and the supernatant was analyzed by gas chromatography (GC) or gas chromatography-mass spectrometry (GC-MS) according to previously reported methods [25].

DNA and RNA extraction
Strain C-6was cultured with BT or sodium sulfate as the sole sulfur source. Cells in culture were harvested by centrifugation (1,200 ×g, 10min) at the end of the logarithmic growth phase. Following incubation of cells at 37°C for 10 min with 125 μg/ml lysostaphin (AMBI, New York, USA), genomic DNA extraction was carried out with a bacterial genomic DNA purification kit (EdgeBio, Maryland, USA) according to the manufacturer's instructions. Cells were resuspended with TRIzol and disrupted with Mini-Bead beater (Biospec, California, USA).Total RNA was isolated according to the manufacturer's instructions of TRIzol-chloroform extraction (Invitrogen, California, USA). Genomic DNA was removed with Turbo DNase. Ribosomal RNA (rRNA) was removed with the Ribo-Zero Magnetic kit (G +/G-Bacteria) (Epicentre, Wisconsin, USA).

RNA sequencing and data analysis
The cDNA library was constructed with extracted mRNA by the TruseqTM RNA sample prep kit (Illumina, California, USA). Sequencing was carried out on the Illumina HiSeq 2000 platformat Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China.
The RNA-seq data from both of the samples were assembled by Trinity [32], and unigenes were predicted with Metagene (prokaryotic gene identification from environmental genome shotgun sequences).The functions of unigenes were annotated by searching against the non-redundant protein database with Blastp with E values less than 1.0×10 −5 . The clusters of orthologous groups (COG) functional categories were assigned with the STRING database (String v9.1). Reads from each sample were mapped to unigenes by BowTie software [33]. The expression of each unigene was calculated with the numbers of reads mapping to it. Differences in gene expression were analyzed by the R package, edgeR [34]. Genes with Q value <0.1 and fold-change 2were assigned as differentially expressed. Sequences were deposited in the NCBI SRA database under the accession no. SRA092218.

Cloning and expression of bds operon
The bds operon was amplified from the genomic DNA of strain C-6 by PCR with the following primers: 5'-ACACCATGGCTGAGGACGAAACCCCGATGACC-3' and 5'-GTGTGGAAGCTTGTTCTGTGGCAGGGGCTTCAG-3'. Sequence was deposited in the Genbank under the accession no. KC831580.The amplified product was ligated into the NcoI and HindIII sites of pET28a (+) to yield pET28a-bdsABC. This recombinant plasmid was introduced into E. coli Rosetta (DE3) by electroporation with a Micropulser electroporator (Bio-Rad, California, USA) as described previously to yield E. coli Rosetta (DE3) overexpressing pET28a-bdsABC [35].

Enzymatic assays
E. coli Rosetta (DE3) overexpressing pET28a-bdsABC was cultured in BSM S-with 0.3mM sodium sulfate as the sole source of sulfur (34 mg/L Kanamycin ) at 37°C to an optical density (OD) of 0.8 at 600 nm. The culture was then induced with 0.06 mM IPTG to express enzymes encoded by bdsABC. In order to avoid the formation of inclusion bodies, the cultivation temperature was reduced to 16°C. The cells were harvested by centrifugation (12,000×g, 5 min) at 4°C and washed with 0.85% (w/v) sodium chloride solution after incubation for 48h. Cells extracts were prepared as previously described [28]. Enzymatic assays were performed with 0.5ml cell extract supernatants, 0.5mM BT, 20 μM FMN, 1mM NADH, 10nM DszD purified from E. coli (DE3)-pET28a-dszD in100 mM pH 7.2 KPi buffer (total of 5 ml) at 30°C. Metabolites of this reaction system were analyzed with GC-MS according to previously reported methods [25].

Quantitative real-time PCR (RT-qPCR) verification
37 genes proposed involving BT biodesulfurization were chosen for the confirmation of RNA-seq data by qPCR with the SYBR Premix Ex Taq kit (Takara, Japan) according to the manufacturer's instructions with a real-time thermal cycler (Bio-Rad, Hercules, CA). Templates used for qPCR were cDNAs that was inverse transcribed from the extracted mRNA using M-MLV Reverse Transcriptase (Takara, Japan) with random primers. The primers used for RT-qPCR detection of selected genes are listed in Table S3 in File S1. The relative gene expression data were analyzed using the 2 -ΔΔCt method as described previously [36]. The results were analyzed using a one-way analysis of variance (ANOVA) statistical test. All quantitative PCR were repeated in three biological and three technical replications.

Supporting Information
File S1. Figure S1, Histogram of clusters of orthologous groups (COG) classification. A, RNA processing and modification; B, Chromatin structure and dynamics; C, Energy production and conversion; D, Cell cycle control, cell division, chromosome partitioning; E, Amino acid transport and metabolism; F, Nucleotide transport and metabolism; G, Carbohydrate transport and metabolism; H, Coenzyme transport and metabolism; I, Lipid transport and metabolism; J, Translation, ribosomal structure and biogenesis; K, Transcription; L, Replication, recombination and repair; M, Cell wall/membrane/envelope biogenesis; N, Cell motility; O, Posttranslational modification, protein turnover, chaperones; P, Inorganic ion transport and metabolism; Q, Secondary metabolites biosynthesis, transport and catabolism; R, General function prediction only; S, Function unknown; T, Signal transduction mechanisms; U, Intracellular trafficking, secretion, and vesicular transport; V, Defense mechanisms; W; Extracellular structures; Y, Nuclear structure; Z, Cytoskeleton. Figure S2, MS spectra of the phenolic compound produced by cell extracts of E. coli Rosetta (DE3) conceived with pET28a-bdsABC during BT biodesulfurization. Metabolites of BT in the cell extracts reaction system was extracted directly with ethyl acetate without adjusting its pH to 2.0, and analyzed by GC-MS according to the previously described method [25]. o-hydroxystyrene was identified by MS spectra corresponding to the peak with retention time of 3.67 in the GC profile. The molecular weight of o-hydroxystyrene decreased by 1 (the hydroxyl of o-hydroxystyrene exists in the form of negative ions) due to the extracting procedure. Table S1, The amino acid homology of BT/DBT biodesulfurization enzymes. Table S2, in File S1 Sulfur resource specificity of Gordonia terrae C-6. Table S3, in File S1 Primers used for RT-qPCR.