Journal of Bioscience and Bioengineering VOL. 116 No. 3, 309e312, 2013 www.elsevier.com/locate/jbiosc
Genomic sequencing-based detection of large deletions in Rhodococcus rhodochrous strain B-276 Seikoh Saitoh,1, y Hiroaki Aoyama,1, 2, y, * Masako Akutsu,1 Kazuma Nakano,1 Naoya Shinzato,1 and Toru Matsui1 Tropical Biosphere Research Center, University of the Ryukyus, 1 Senbaru, Nishihara, Okinawa 903-0213, Japan1 and Laboratory of Cell and Functional Biology, Faculty of Science, University of the Ryukyus, 1 Senbaru, Nishihara, Okinawa 903-0213, Japan2 Received 1 November 2012; accepted 2 March 2013 Available online 6 April 2013
Bacteria of the genus Rhodococcus (Actinomycetes) have the ability to catabolize various organic compounds and are therefore considered potential genetic resources for applications such as bioremediation. We investigated a next-generation sequencing-based procedure to rapidly identify candidate functional gene(s) from rhodococci on the basis of their frequent genome recombination. The Rhodococcus rhodochrous strain B-276 and its alkene monooxygenase (AMO) gene cluster were the focus of our investigation. Firstly, 2 types of cultures of the R. rhodochrous strain B-276 were prepared, one of which was supplied with propene, which requires AMO genes for its assimilation, whereas the other was supplied with glucose as the sole energy source. The latter culture was anticipated to have a lower gene frequency of AMO genes because of their deletion during cultivation. We then conducted whole genome shotgun sequencing of the genomic DNA extracted from both cultures. Next, all sequence data were pooled and assembled into contiguous sequences (contigs). Finally, the abundance of each contig was quantified in order to detect contigs that were highly biased between the 2 cultures. We identified contigs that were overrepresented by 2 orders of magnitude in the AMO-required culture and successfully identified an AMO gene cluster among these contigs. We propose this procedure as an efficient method for the rapid detection and sequencing of deleted region, which contributes to identification of functional genes in rhodococci. Ó 2013, The Society for Biotechnology, Japan. All rights reserved. [Key words: Rhodococcus; Pyrosequencer; 454 GS Junior; Deletion detection; Genome comparison]
Rhodococcus (Actinomycetes) are known for their ability to degrade a broad range of organic compounds, including short- and long-chain alkanes and aromatic compounds (1). Their catabolic activities have attracted interest for applications such as desulfurization of fossil fuel (2), enantioselective synthesis of compounds (3), and biodegradation of hazardous pollutants such as polychlorinated biphenyl (4), trichloroethene (5), and thiodiglycol, a chemical weapon hydrolyzate (6). These functions are often mediated by gene products encoded by large plasmids harbored in the cell (7e10). Microorganisms readily undergo loss of genes not required for survival under specific environmental conditions (11). Under such conditions, mutant strains of Rhodococcus often arise, which are deficient in certain gene(s) due to loss of an entire plasmid or due to partial deletion of the encoding region. For example, Rhodococcus rhodochrous (formerly Nocardia corallina) strain B-276 can grow on propene as the sole carbon and energy source with alkene monooxygenase (AMO), an enzyme responsible for its assimilation (3,5).
* Corresponding author at: Tropical Biosphere Research Center, University of the Ryukyus, 1 Senbaru, Nishihara, Okinawa 903-0213, Japan. Tel.: þ81 98 895 8943; fax: þ81 98 895 8944. E-mail address:
[email protected] (H. Aoyama). y The first two authors contributed equally to this work.
Saeki and coworkers (5) showed that AMO-deficient variants of this strain, obtained under non-selective condition, had undergone partial or complete loss of a 185-kb plasmid (named pNC30) encoding the AMO gene cluster. Identification of such plasmids and analysis of their sequences are therefore important steps for the development of useful Rhodococcus strains. Recent years have seen the growth in popularity of the inexpensive next-generation sequencing (NGS) Roche 454 GS Junior system, which allows small laboratories to avail of high throughput sequencing technology. While its sequencing throughput is an order of magnitude smaller than that of higher grade models such as the 454 GS FLX system, the much smaller cost per run (w15% of the FLX) and its sufficient capacity for analyzing microbial genome (35 Mbp per run) make it a useful tool for microbial genomics. Using contig data from whole genome shotgun (WGS) sequencing performed using this system, we describe here a time-saving procedure for highlighting genomic differences between a wild-type strain and a variant deficient in gene(s) lost by deletion of a coding region in a plasmid. Our procedure is based on the fact that, in a microbial culture growing under conditions in which a plasmid-encoded gene is required for growth, the population of cells harboring the plasmid will be overrepresented compared with another culture grown under conditions in which the gene is not required. This difference
1389-1723/$ e see front matter Ó 2013, The Society for Biotechnology, Japan. All rights reserved. http://dx.doi.org/10.1016/j.jbiosc.2013.03.002
310
SAITOH ET AL.
J. BIOSCI. BIOENG.,
can be detected by comparison of WGS data between the 2 cultures with respect to the number of reads corresponding to the plasmid (Fig. 1). We sequenced WGS libraries prepared from 2 cultures of the R. rhodochrous strain B-276 whose sole energy sources were either propene, which requires AMO gene for its assimilation, or glucose [hereafter referred to as the selective condition (SC) or the non-selective condition (NSC), respectively]. Sequence data from the 2 WGS libraries were pooled and assembled to construct consensus sequences (contigs). The number of reads comprising each contig was counted, and the difference between the 2 libraries was analyzed statistically to detect overrepresented contigs. Finally, the overrepresented contigs were analyzed for the presence of the AMO gene cluster. MATERIALS AND METHODS Bacteria and growth conditions Cultures of R. rhodochrous B-276 (ATCC31338) were grown on Nutrient Broth no. 2 (Oxoid, Basingstoke, Hampshire, UK) containing 1% (w/v) glucose for the NSC, which was the condition used by Saeki and coworkers (5). For the SC, propene was supplied with propene (10%, v/v) via the gas phase with incubation in minimal medium (12). Cultures were grown in shaking flasks at 30 C with aerobic shaking (180 rpm) under dark condition for a week. Cells
Bacterial Strain (harboring a functional plasmid)
(1) Culturing a bacterial strain under 2 different conditions
Culture under selective condition (i.e., propene)
Culture under nonselective condition (i.e., glucose)
(2) Sequencing both cultures by the 454 technology
Read sequences of the selective condition
Read sequences of the nonselective condition
TABLE 1. Estimated proportions of cells capable to grow on propene of cultures after incubation under SC and NSC. Culture
NSC SC
Number of colonies
Grown
20 22
Proportion (%) Estimated
95% CIa
0 100
0e17 85e100
0 22
Summary of growth test in which propene was supplied as the sole energy source was shown (see text for the procedure). a 95% confidential interval based on exact confidential interval (13). were harvested by centrifugation (8000 g for 10 min), washed twice with sterilized water, followed by re-centrifugation (8000 g for 10 min) to obtain cell pellets for DNA extraction. Cell pellets were stored at 80 C until use. Examination of proportion of cells capable to utilize propene in the cultures after incubation Culture after incubation under SC or NSC was examined in terms of proportion of cells capable to grow on propene as the sole energy source. The test was conducted using an air-tight glass vial with 3 mL of the minimal medium, 28.5 mL air, an inoculum from a colony of SC or NSC, and 20 mL of propene injected using a syringe. Each vial was incubated at 30 C with aerobic shaking. The criterion for the test was whether OD600 measured after incubation for 96 h exceeded 0.8 or not. A 95% confidential interval of the proportion was calculated using exact confidential interval (13). DNA extraction and shotgun pyrosequencing Bacterial genomic DNA was extracted using a DNA extraction kit, ISOPLANT I (Nippon Gene, Tokyo, Japan). After extraction, the DNA was quantitatively and qualitatively analyzed by a Nanodrop (Thermo Fisher Scientific, Wilmington, DE, USA) and Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). NSC and SC DNA samples were subjected to shotgun pyrosequencing by using a 454 GS Junior titanium chemistry (454 Life Sciences-a Roche Company, Branford, CT, USA), and DNA libraries were thus prepared. Pyrosequencing runs were carried out in duplicate for the NSC library and once for SC library. All operations were carried out according to the protocols provided by the manufacturers. In silico analyzes Sequence data were pooled and assembled using Roche Newbler software (ver. 2.6) (14) with default settings. The abundance of DNA corresponding to the contig sequence in the 2 libraries was derived from the number of reads. For standardization purposes, the number of reads per 1 kbp sequence per 1 M aligned reads (hereafter referred to as RPKM) was used to compare the abundance of sequences between the 2 libraries. This analysis was performed using an in-house script in the R statistical software (ver. 2.13, R Core Development team). The difference of abundance between the 2 libraries for each contig was examined using the Fisher exact test based on the total aligned number of reads. P-values were adjusted with Bonferroni correction. The Basic Local Alignment Search Tool (BLAST; ver. 2.2.25) (15) was used to identify contigs containing the AMO gene cluster (GenBank: D37875.1) (16). The bb.454contignet Perl script (17) was then used to illustrate paths between the contigs assembled by the Newbler software. A putative annotation of the constructed draft genome sequence was made using Microbial Genome Annotation Pipeline (MiGAP, an auto annotation pipeline of DDBJ, http://www.migap.org/) for discussion on similarity of the genome to other genomes of rhodococci sequenced ever. The sequence data from the present study have been submitted to the DDBJ Sequence Read Archive (DRA) under accession no. DRA000810.
RESULTS AND DISCUSSION Contig A
Contig B
Contig C
Contig D
(3) Pooling the 2 sequence data sets and assembling into contigs
In the present study, we cultured R. rhodochrous strain B-276 grown under 2 different conditions (NSC and SC). Result of growth
Reads of the selective condition
TABLE 2. Statistic of sequencing runs using 454 GS Junior. Sequencing run no. 1
Contig B
Reads of the nonselective condition
(4) Counting the number of reads for each contig FIG. 1. Conceptual scheme of the procedure to detect large deletions in a bacterial genome, as suggested in this study: 1, A bacterial strain is cultured under 2 different conditions (feeding propene or glucose as the energy and carbon sources); 2, genome sequence from both cultures is analyzed by whole genome shotgun sequencing by using the 454 pyrosequence technology; 3, the 2 sequence datasets are pooled and assembled to construct a set of contigs; 4, the number of reads for each contig is counted, and the contig(s) whose abundance of read sequences is highly biased between the 2 libraries would be detected.
Librarya Sequencing throughput Total number of reads Total number of bases Assembly resultb Number of aligned reads Number of aligned bases % Aligned base a
2
3
NSC
NSC
SC
127157 57886794
126816 59169882
121990 55981563
126048 57484861 99.31
125660 58743091 99.28
121167 55750732 99.59
The non-selective condition (NSC) library or the selective condition (SC) library. Data from all 3 runs were pooled and assembled by Roche Newbler (ver. 2.6) using default settings. b
VOL. 116, 2013
DELETION DETECTION IN RHODOCOCCUS BY GENOME SEQUENCING TABLE 3. Statistic of assembly using Newbler.
Peak depth Number of contigs Total length of contigs (bp) N50 contig size (bp)a % Bases of QV 40
21 159 6576242 144721 99.91
Data from all 3 runs were pooled and assembled by Roche Newbler (version 2.6) using default settings. Threshold of the minimal contig size was 100 bp. a N50 contig size means that half number of total aligned reads were assembled into this size of contigs or greater.
contig00061
contig00079 contig00109
1e+01
contig00127
1e-01
1e+00
contig00129 contig00126
1e-02
Ratio of read sequence abundance
1e+02
contig00062 contig00073
1e+02
1e+03
1e+04 1e+05 Length of contig (bp)
1e+06
FIG. 2. Scatter plot of the ratio of contig abundance (of SC to NSC library, normalized by RPKM) against contig length (bp). Contigs relatively abundant in the SC and NSC libraries are plotted in the upper and lower halves of the coordinate plane, respectively. Empty and filled symbols denote statistical significance (Fisher exact test with Bonferroni correction, Padj < 0.05) and no-significance (Padj 0.05), respectively. Triangles denote data whose zero-values in the abundance were substituted by 1 before calculation of the ratio of abundance for logarithmic plot or avoidance of division by zero (for actual values, see Table S1).
test showed that the 2 cultures highly differed in proportion of cells capable to grow on propene as the sole energy source (Table 1; p < 0.001, Fisher exact test). We prepared DNA libraries from the 2 cultures and sequenced using a Roche 454 GS Junior sequencer. The sequencer yielded 173 Mbp from 3 sequencing runs, comprising 2 runs from the NSC library and 1 run from the SC library (Table 2). Sequence data were pooled and assembled using the Newbler software. The assembler gave an output of 159 contigs (>100 bp) in all, with total length of approximately 6.6 Mbp, a peak depth of 21, and N50 of 144.7 kbp (Table 3).
311
Comparison of contig abundance between the 2 libraries based on RPKM showed that 6 contigs, ranging in length from 913 bp to 25 kbp, were highly abundant in the SC library relative to the NSC library (Fig. 2). Contigs overrepresented by a factor of 10 are listed in Table 4 (for the entire listing, see Table S1). Among these, we successfully identified a contig encoding an AMO gene cluster (contig00061), consistently with the result of growth test described earlier. Interestingly, 5 of the 6 contigs had paths between themselves (with the exception of contig00127) in the mathematical graph as shown in Fig. 3, indicating that they represent consecutive parts of the plasmid encoding the AMO gene cluster. The fact that these were separate contigs, rather than a single one, may be attributed to the presence of similar sequences within the plasmid and/or between the plasmid and other genomic components (i.e., the chromosome and/or the other plasmids), which prevented solving the paths in the mathematical graph constructed by the Newbler assembler and resulted in separate contigs in the output. This notion is supported by the identification of 2 of the contigs, as shown in Fig. 3 (contig00126 and contig00129), as transposonderived sequences, which are typically multicopy genes or pseudogenes in Rhodococcus genomes (1). In general, multicopy sequences larger than a single read sequence (approximately 500 bp in this study) cause ambiguity in assembly computation (18). Determining the entire structure of the genome requires additional data such as that generated by Sanger sequencing over multicopy genes, or NGS data from paired-end libraries. From the viewpoint of the evolutionary strategy of Rhodococcus (1), transposon-derived sequences enable the organism to rapidly adapt to a wide range of environmental conditions via frequent genomic homologous recombination. Saeki and coworkers (5) estimated the size of the pNC30 plasmid, which encoded the AMO gene cluster, at approximately 185 kbp. In the present study, however, the total length of highly abundant SC library contigs was only approximately 79 kbp. This anomaly suggests a deletion of 79 kbp in the plasmid, likely via a recombination at one of the transposon-derived sequences (contig00126 or contig00129), rather than loss of the entire plasmid, which might have occurred in cells grown under NSC. The draft sequence of genome of R. rhodochrous strain B-276 includes 6179 putative coding sequences (CDS). Of these, 2287, 96, and 55 CDSs were homologous to known genes from Rhodococcus, Nocardia, and Gordonia strains, respectively. A total of putative 71 monooxygenase genes and 132 dioxygenase genes whose substrates include benzoate, biphenyl, catechol, extradiol, naphthalene, phthalate and their some derivatives were identified. Such information from the genome sequence will be worthy for further researches on this strain (Aoyama et al., in preparation). In the present paper, we have demonstrated the effectiveness of comparative analysis of 2 genome sequence datasets from 2 cultures under different growth conditions, one of which required specific functional genes contained within a plasmid (i.e., the AMO
TABLE 4. Contigs highly abundant in the selective condition (SC) library. Contig
Contig00062 Contig00061 Contig00073 Contig00079 Contig00109 Contig00127
Length (bp)
23329 24974 15947 12888 1725 914
Abundancea
Number of reads NSC
SC
NSC
SC
13 12 13 13 6 0
1292 1144 859 696 213 8
2 2 3 4 14 0
457 378 445 446 1019 72
Ratiob
Padjc
Note
206.4 198.0 137.3 111.2 73.7 Inf
<0.001 <0.001 <0.001 <0.001 <0.001 0.020
Included in Fig. 3 Coding the AMO gene cluster Included in Fig. 3 Included in Fig. 3 Included in Fig. 3 Suspected as contaminant DNA
Contigs overrepresented in the SC library by a factor of 10 are listed. For information on all contigs, see Table S1. a Number of reads per 1 kbp sequence per 1 M aligned reads (RPKM). b Ratio of RPKM of SC to NSC library. c P-value of Fisher exact test adjusted with Bonferroni correction.
312
SAITOH ET AL.
J. BIOSCI. BIOENG.,
FIG. 3. Contig connections of de novo assembly of R. rhodochrous strain B-276 around the contigs coding the AMO gene cluster (contig00061, abbreviated to c00061) drawn by bb.454contignet Perl script (15) based on the output from the Roche Newbler assembler. Boxes denote contigs with the information of the contig ID, length of contig sequence, and sequencing coverage (depth). Lines between boxes denote paths, each of which indicates that there were reads overlapping both of the contigs linked by the line. The numbers accompanying the lines denote the number of reads linking the 2 contigs.
gene cluster on pNC30). De novo assembly and quantification of abundance of reads for each contig allows for successful identification of the region of the genome containing the gene in question. This strategy simultaneously provides information on the whole genome sequence of organism, allowing researchers to carry out genome-wide analyzes between the plasmid containing the gene of interest and other genomic components (e.g., chromosome DNA and other plasmids). Our procedure represents a cost-effective approach to analyzing the genomic structure of microorganisms such as Rhodococcus that undergo frequent genome rearrangements. Supplementary data related to this article can be found online at http://dx.doi.org/10.1016/j.jbiosc.2013.03.002. References 1. Larkin, M. J., Kulakov, L. A., and Allen, C. C. R.: Biodegradation and Rhodococcus e masters of catabolic versatility, Curr. Opin. Biotechnol., 16, 282e290 (2005). 2. Matsui, T., Onaka, T., Tanaka, Y., Tezuka, T., Suzuki, M., and Kurane, R.: Alkylated benzothiophene desulfurization by Rhodococcus sp. strain T09, Biosci. Biotechnol. Biochem., 64, 596e599 (2000). 3. Furuhashi, K., Taoka, A., Uchida, S., Karube, I., and Suzuki, S.: Production of 1,2-epoxyalkanes from 1-alkenes by Nocardia corallina B-276, Eur. J. Appl. Microbiol. Biotechnol., 12, 39e45 (1981). 4. Seto, M., Kimbara, K., Shimura, M., Hatta, T., Fukuda, M., and Yano, K.: A novel transformation of polychlorinated biphenyls by Rhodococcus sp. strain RHA1, Appl. Environ. Microbiol., 61, 3353e3358 (1995). 5. Saeki, H., Akira, M., Furuhashi, K., Averhoff, B., and Gottschalk, G.: Degradation of trichloroethene by a linear plasmid-encoded alkene monooxygenase in Rhodococcus corallinus (Nocardia corallina) B-276, Microbiology, 145, 1721e1730 (1999). 6. El Bassi, L., Shinzato, N., Namihira, T., and Matsui, T.: Biodegradation of thiodiglycol, a hydrolyzate of the chemical weapon Yperite, by benzothiophene-desulfurizing bacteria, J. Hazard. Mater., 167, 124e127 (2009). 7. Kesseler, M., Dabbs, E. R., Averhoffl, B., and Gottschalk, G.: Studies on the isopropylbenzene 2,3-dioxygenase and the 3-isopropylcatechol 2,3-
8.
9.
10.
11.
12.
13. 14.
15. 16.
17.
18.
dioxygenase genes encoded by the linear plasmid of Rhodococcus erythropolis BD2, Microbiology, 142, 3241e3251 (1996). Kosono, S., Maeda, M., Fuji, F., Arai, H., and Kudo, T.: Three of the seven bphC genes of Rhodococcus erythropolis TA421, isolated from a termite ecosystem, are located on an indigenous plasmid associated with biphenyl degradation, Appl. Environ. Microbiol., 63, 3282e3285 (1997). Masai, E., Sugiyama, K., Iwashita, N., Shimizu, S., Hauschild, J. E., Hatta, T., Kimbara, K., Yano, K., and Fukuda, M.: The bphDEF meta-cleavage pathway genes involved in biphenyl/polychlorinated biphenyl degradation are located on a linear plasmid and separated from the initial bphACB genes in Rhodococcus sp. strain RHA1, Gene, 187, 141e149 (1997). Shimizu, S., Kobayashi, H., Masai, E., and Fukuda, M.: Characterization of the 450-kb linear plasmid in a polychlorinated biphenyl degrader, Rhodococcus sp. strain RHA1, Appl. Environ. Microbiol., 67, 2021e2028 (2001). Barton, N. H., Briggs, D. E. G., Eisen, J. A., Goldstein, D. B., and Patel, N. H.: Evolution. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY (2007). Matsui, T. and Furuhashi, K.: Asymmetric oxidation of isopropyl moieties of aliphatic and aromatic hydrocarbons by Rhodococcus sp. 11B, Biosci. Biotechnol. Biochem., 59, 1342e1344 (1995). Clopper, C. J. and Pearson, E. S.: The use of confidence or fiducial limits illustrated in the case of the binomial, Biometrika, 26, 404e413 (1934). Margulies, M., Egholm, M., Altman, W. E., Attiya, S., Bader, J. S., Bemben, L. A., Berka, J., Braverman, M. S., Chen, Y. J., Chen, Z., and other 46 authors: Genome sequencing in microfabricated high-density picolitre reactors, Nature, 437, 376e380 (2005). Altschul, S., Gish, W., Miller, W., Myers, E., and Lipman, D.: Basic local alignment search tool, J. Mol. Biol., 215, 403e410 (1990). Saeki, H. and Furuhashi, K.: Cloning and characterization of a Nocardia corallina B-276 gene cluster encoding alkene monooxygenase, J. Ferment. Bioeng., 78, 399e406 (1994). Iorizzo, M., Senalik, D., Szklarczyk, M., Grzebelus, D., Spooner, D., and Simon, P.: De novo assembly of the carrot mitochondrial genome using next generation sequencing of whole genomic DNA provides first evidence of DNA transfer into an angiosperm plastid genome, BMC Plant Biol., 12, 61 (2012). Schatz, M. C., Delcher, A. L., and Salzberg, S. L.: Assembly of large genomes using second-generation sequencing, Genome Res., 20, 1165e1173 (2010).