A Roadmap for Natural Product Discovery Based on Large-Scale Genomics and Metabolomics

Creation and analysis of a GCF network

We created a dataset comprised of NPGCs involved in the synthesis of NRPS, type I and type II PKS, NRPS-independent siderophores (NISs), lanthipeptides, and thiazole-oxazole modified microcins (TOMMs). These include NPGCs identified in 830 actinomycete genomes, 344 of which were sequenced for this project (Supplementary Data Set 1), as well as 412 gene clusters with established natural products from Genbank. Three distance metrics were calculated for every NPGC pair: (i) the number of homologous genes shared (Fig. 1a); (ii) the proportion of nucleotides involved in a pairwise alignment (Fig. 1b); and (iii) the amino acid sequence identity between the domains of repeated protein modules (Fig. 1c). When used alone, each of these metrics has limited ability to separate homologous from paralogous gene clusters. For example, BLAST-based ortholog detection and clustering as reported previously ¹⁶ is especially limited by the repeated, modular structure of large PKS and NRPS genes, which produce very low e-values for genes in different gene cluster families. However, a combined score that incorporates all three metrics produces coherent GCFs that include only highly related NPGCs, as assessed by manual inspections of gene cluster diagrams and GCF network visualizations (Supplementary Results, Supplementary Fig. 1; an interactive version of Supplementary Figure 1 is available at www.igb.illinois.edu/labs/metcalf/gcf). The all-v-all nature of this approach was very computationally expensive, taking several months on a server utilizing 40 processors, mostly dedicated to PROmer comparisons. The combined scoring metric was not systematically compared to the BLAST only clustering because the families produced by this metric clearly produced families containing unrelated gene clusters (data not shown). To further validate the approach, we examined the gene cluster families that contain multiple NPGCs with established natural products. (Note that characterized gene clusters that are the sole members of GCFs are not informative in this regard and were not included in the analysis.) The combined metric correctly grouped 103 characterized gene clusters into forty-one GCFs that direct synthesis of highly similar natural products (Supplementary Note and Supplementary Fig. 2). The single known NPGC that was not grouped with relatives making highly similar compounds is an unusual actinomycin gene cluster with a large internal duplication, previously reported to be divergent from other actinomycin gene clusters ¹⁷.

The GCF network created using these methods from the 830 actinobacterial genomes was comprised of 140,986 genes from 11,422 NPGCs grouped into 4,122 GCFs (Fig. 2a). There were 77 GCFs that contain at least one characterized NPGC along with related, as yet uncharacterized, gene clusters in the newly sequenced genomes. These “anchored” GCFs provided 1193 new NPGCs that are likely to direct the synthesis of novel derivatives of known compounds (Supplementary Data Set 2), each with a high expectation of exhibiting a targeted bioactivity. We are currently developing functionality that will allow users to submit their own genomic data into this framework.

The ability of the GCF approach to recognize whether one NPGC is the same or different from another allows the use of rarefaction analysis to estimate the total number of natural product scaffolds made by actinomycetes and the number of strains that would need to be screened to find them all. This statistical method examines the rate at which new GCFs accumulate as new genomes are randomly added to the data set. This process is somewhat complicated by the non-random nature of the data; over 40% of the genomes in our dataset were derived from a few medically relevant genera: Mycobacterium (142 genomes); Corynebacterium (68 genomes); Propionibacterium (89 genomes); Actinomyces (24 genomes); and Gordonia (20 genomes (Fig. 2a). To correct for these phylogenetic biases, we randomly sampled our data from phylogenetic bins (operational taxonomic units, or OTUs) for extrapolation analyses and to produce a rarefaction curve (Fig. 2b). Three different estimation methods, Chao1, ACE and extrapolation, implemented within the software EstimateS ¹⁸, were used to analyze the data; all gave similar results (Supplementary Table 1). Although the NRPS GCFs were the most abundant in our current data, the extrapolation predicts that type I PKS will become more abundant as more genomes are sampled. These classes were followed, in order, by type II PKSs, lanthipeptides, NRPS-independent siderophores (NISs), and TOMMs. Collectively, we estimated that the actinomycetes encode ~17,350 GCFs in these six classes. Because each GCF represents the biosynthetic capacity to produce a suite of related molecules, these GCFs are predicted to direct the synthesis of hundreds of thousands, but not millions, of natural products. The extrapolations predict that essentially all actinomycete NPGCs will be identified after ca. 15,000 genomes are sequenced from different OTUs, an achievable proposition given the decreasing cost of sequencing. Considerably fewer genome sequences will be required to access the full repertoire of the less abundant classes (Fig. 2b).

Our analysis of the abundance and distribution of the GCFs in the dataset provided an efficient path for selection and prioritization of actinobacterial strains for future natural product discovery efforts. Consistent with previous reports¹⁹, our data suggested that certain phylogenetic groups are gifted with an abundance and diversity of NPGCs (taxa between Streptomycetales to Pseudonocardiales, inclusive, in Fig. 2a). In total, these gifted groups encoded 73.6% of the GCFs present in our dataset (Table 1). Based on these finding, we asked whether inclusion of the non-gifted genomes biased the diversity estimates. However, when only the gifted taxa were analyzed, the predicted numbers of GCFs were essentially the same as presented in the previous section (Supplementary Table 2). Importantly, our data showed a strong correlation between phylogenetic distance and the GCF complement of individual strains, with the degree of conservation varying according to biosynthetic class and phylogenetic group (Fig. 3, Supplementary Fig. 3). Among the gifted taxa, related strains share 80% of their NRPS and 73% of their type I PKS GCFs at the 1% divergence level (based on a concatenated ribosomal protein phylogeny), These values drop sharply beyond the 1% cut-off, such that strains share only 6% and 3% of their respective NRPS and PKS I GCFs at a ribosomal protein divergence of 4-5%. We observed similar trends for most other natural product classes, although the rate of decline in GCF conservation across genetic distance for lanthipeptides and type II PKS clusters was not as pronounced (Supplementary Fig. 3). In contrast, NRPS-independent siderophores were conserved across a much larger phylogenetic distance: 84% and 66% conservation at 0-1% and 4-5% ribosomal protein divergence, respectively. Importantly, we observed the same correlations between phylogeny and GCF with the use of a simple 297 bp fragment of the rpoB gene suitable for amplicon sequencing (Fig. 3b). Significantly, 1% ribosomal protein divergence corresponds to a suggested “species-level” cutoff within Streptomyces (Supplementary Fig. 4), consistent with the idea that natural product biosynthetic capacity drives speciation within this genus²⁰.

Correlation of GCFs with secondary metabolite production

Our analysis of GCFs substantially reinforces the idea that Actinobacteria represent a major source of new natural products; however, to move beyond the concept, the molecules produced by each GCF must be identified. With this in mind, we developed a liquid chromatography-high resolution mass spectrometry (LC-HRMS) method to quantitatively measure the exported metabolome of the strains being screened. We tested this method on pooled extracts produced by 178 of the strains included in our GCF data set after growth in four different media. After background subtraction, we detected an average of 105 compounds produced by each strain, many of which were found in multiple strains/samples. In total, 2,521 unique compounds (intact masses) were identified in the MS¹ data. Among these, we identified 110 previously characterized natural products after automated search of a database consisting of 9,817 known actinomycete natural products. After we collapsed highly similar molecules (e.g. actinomycin A and actinomycin D) into groups, we experimentally linked 27 known natural products to nine GCFs in our network. In the samples we analyzed, these 27 compounds were observed and verified by high-resolution tandem MS 268 times (representative data shown in Supplementary Fig.5). Significantly, we identified the correct GCF in the genome of the source organism in 92% of these cases. Thus, the method has a relatively low rate of false positives. Interestingly, these same gene clusters could be found in many strains that did not produce the compounds in question, i.e. the gene clusters were cryptic in these cases. In total, known GCFs were cryptic in 77% of the strains we examined. This observation belies the idea that NPGCs are usually cryptic, with major implications for discovery efforts. Accordingly, one may need to screen on average only four isolates with a given NPGC to find one in which the corresponding natural product is produced at levels sufficient to allow mass spectrometric detection. We wish to emphasize that our assessment of cryptic NPGCs applies only to the known compounds that we identified. It remains to be seen whether this trend will apply to all NPGCs.

We were particularly interested in whether an automated method could be developed for linking specific gene clusters to specific molecules (i.e. specific exact masses from the MS data) based on the simple fact that a molecule cannot be produced without the biosynthetic genes also being present. To achieve this, we performed a binary correlation between subsets of each GCF - based on clustered domain sequences - present in each genome and the MS¹ intact masses found in each extract. Without manual intervention, the binary correlation of the GCF and MS data sets associated experimentally established biosynthetic gene clusters with production of oxytetracycline, benarthin, nonactin, and enterocin, with correlation scores at the extreme high end of the overall distribution (Fig. 4a, Table 2). Further, manual searches revealed the presence of fragmentary gene clusters (due to the draft nature of the genome sequences) that had been omitted from our automated GFC network in some of the strains that produced MS² verified natural products. This improved the correlation scores for actinomycin D and pyridomycin, as did the use of the full GCF for proferrioxamine D2, such that the manual correlation score was actually higher than the best automated correlation score. Application of the automated method to the full data set also revealed a strong correlation between a previously uncharacterized GCF and the related compounds desertomycin A, B, D, and E and oasamycin A, D, and E (Fig. 4b, Supplementary Fig. 6). This cluster encodes PKS modules consistent with the known molecule, as well as proteins with predicted functions that match the unique structural features of the molecule, specifically N-methylation, and mannosylation. Thus, we believe it highly likely that this previously unassigned GCF is responsible for the synthesis of both desertomycin and oasamycin. We also observed a highly significant correlation between the griseobactin GCF and benarthin. It had previously been speculated, but never proven, that the two molecules share a biosynthetic pathway²¹. Our correlation data provide additional support for this hypothesis, and show that this method can also be used to find biosynthetic intermediates.

Genome Sequencing

All strains were received directly from the Agricultural Research Service (NRRL) culture collection, Peoria, IL. ATCC® Medium #172 broth ²⁵ was used for liquid cultures (per liter: glucose, 10.0 g; soluble starch, 20.0 g; yeast extract, 5.0 g; N-Z amine type A, 5.0 g; CaCO₃, 1.0 g; agar, 15.0 g). Genomic DNA was prepared using the UltraClean Microbial DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA) and sequencing library preparation was performed with Nextera version 2 kits (Illumina, San Diego, CA). Genomes for the 343 strains reported here were sequenced at the University of Illinois at Urbana-Champaign Keck Sequencing Center and the University of Wisconsin Biotechnology Center. These genomes were sequenced in four separate batches. 87 were sequenced using version 1 Nextera library preparation kits. The 24 genomes handled at the University of Wisconsin Biotechnology Center were sequenced on the Illumina GAIIx platform. Sixty-four genomes, including one that had insufficient reads on the GAIIx, were sequenced on an Illumina HiSeq using v2 chemistry at the University of Illinois. For these 87 genomes, an initial set of genome assemblies was performed using Velvet version 1.0.15 ²⁶, SOAPdenovo v1.04 ²⁷, EULER-SR v1.1.2 ²⁸. Subsequently, faux reads created by breaking larger contigs from these assemblies into 1999 bp pieces and 400,000 reformatted Illumina paired-end reads were used as input for gsAssembler v2.5.3 ²⁹. The resulting contigs were then used as input, along with all Illumina reads for a final scaffolding step with SSPACE v1.1 ³⁰. For the remaining genomes, sequencing was performed on an Illumina HiSeq 2000 with version 3 chemistry and 24 samples pooled per lane. These genomes were assembled IDBA UD version 1.0.9 ³¹. Gene prediction was performed with Prodigal version 2.5.0 ³². The remaining 487 genomes were downloaded from NCBI. All available gene predictions were kept and not altered. NCBI taxonomy information was used in an attempt to include all Actinomycetales genomes in NCBI as of February 2013. A small number of these genomes were not carried all the way through the pipeline, primarily due to a lack of complete uniformity in locus tag and formatting conventions. All further processing was performed with Perl scripts interfaced with MySQL.

Bioinformatics

Gene clusters were identified as previously described ¹⁹. Orthologs were defined using OrthoMCL to analyze an initial set of 231 genomes ³³. Genes included within gene cluster boundaries for additional genomes were recruited to this initial data set by defining orthologs based on best BLAST hits to the original 231 genome data set below an e-value cutoff of 1e-10.

Gene cluster families were created using three similarity measures. The first is the mean proportion of orthologs shared by two gene clusters, with each gene only being counted once to prevent inflated scores from repeated genes (e.g. polyketide synthases and nonribosomal peptide synthetases). The cutoff used for this score was 0.5 (50%). The second uses PROmer alignments of every possible pairwise combination of gene clusters that were scored as the mean proportion of each gene cluster that is part of an alignment ³⁴. The cutoff used for this score was 0.5 (50%). The third is based on similarity of one type of domain or full length gene for each biosynthetic class as follows: type I and II polyketide synthases, ketosynthase domain; nonribosomal peptide synthetases, adenylation domain; NISiderophores, IucA-IucC domains; lanthipeptides, lanthionine-cyclase containing genes; and microcins, dehydratase genes. The program uclust was used to group domains at a 70% similarity cutoff ³⁵. The cutoff used for this score was that half of the total number of key domains or genes present in a pair of clusters must be in the same 70% similarity groups. To be placed into a gene cluster family, all three of these cutoffs must be passed. These scores were then weighted, counting the domain similarity twice, and converted to a distance metric as D = a + h + 2s, where D is the distance metric, a is the alignment score, h is the proportion of shared homologous genes and s is the highest clustering threshold in uclust that groups at least half of the domains in a pair of clusters together. These distance scores were used as input for density-based clustering with DBSCAN ³⁶ implemented in the R package fpc (http://cran.r-project.org/web/packages/fpc/) with the following parameters: eps, 0.3; and MinPts, 2. The clusters assignments produced using DBSCAN were used as input to color the gene cluster family networks in the program Cytoscape to aid manual auditing ³⁷. All GCFs were manually examined, with special attention paid to GCFs split into separate clusters by DBSCAN.

Lanthipeptides and TOMMs were subjected to additional scrutiny based on predicted precursor peptides. The lan-cyclase domain containing genes and the cyclodehydratase, previously known as the docking protein, genes were aligned with MUSCLE v3.7 ³⁸ and used as input for FastTree version 2.1.5 SSE3, OpenMP ³⁹ with option –gamma for rate optimization under the Gamma20 model. The resulting phylogenetic framework was used to examine precursor peptides, with special attention paid to likely secondary structure for lanthipeptides, which is tied to resulting function ⁴⁰.

All GCFs were manually checked for similarity based on gene cluster diagrams highlighting orthologous genes. This visualization was performed with in-house Perl scripts interfaced with MySQL and output as HTML with JavaScript. The visualizations are provided in Database S1. The annotations presented of the website are based on the most frequently used annotation for each set of homologous genes, as found within the set of annotated genomes in Genbank. This generic method allows automated annotation based on existing data. We strongly recommend additional bioinformatics analyses for genes of interest.

Ribosomal proteins used to create the phylogenetic tree in Figure 2a were detected using pHMMs created with HMMER using the files developed in a previous study ⁴¹. The proteins were aligned using ClustalW ⁴² and the tree was created with FastTreeMP ³⁹ and visualized with the Python library ETE2 ⁴³. Genomic data was presented with Circos ⁴⁴.

All sequence distance comparisons, for both amino acid and nucleotide sequences) were performed based solely on identity and represent the proportion of differences. These scale linearly and range from 0-1. The rpoB gene fragment used ranges from the Streptomyces coelicolor A3(2) gene SCO4654 coordinates 1451-1747.

Correlations were performed on subsets of each GCF to improve accuracy when looking for an exact mass. An entire GCF likely includes multiple unique structures that share a common core, and this diversity would not be reflected in our current analysis based on intact masses. The central domains or genes discussed above, e.g. ketosynthase domains or full length pepM, were clustered at a threshold of 0.9 using the program uclust³⁵. Individual domains were given unique identifiers, such that some fragmented clusters that are closely related could be used for correlations that were left out of the larger GCF analysis. These subsets of each GCF were the correlated with the binary occurrence of each individual component ID in the output from SIEVE (see below). Scoring was applied in the following fashion: GCF present, component ID present = +10; GCF absent, component ID present = −10; GCF present, component ID absent = 0; GCF absent, component ID absent = +1.

Extrapolations were performed using the program estimateS ¹⁸. Genomes were considered as sampling sites and occurrence of gene cluster families within the genomes as counts of individuals. Analyses were repeated for 100 runs. Classical measures were used for Chao1 and ACE estimates. Extrapolation was performed for a total of 750 knots out to 15000 samples.

Growth and Extraction

Seed cultures for production were grown in yeast-extract, malt-extract medium for 7 days at 30°C. 200 μl of this culture was used as inoculum for four solid media types (all ingredients given per liter, all pH should be adjusted to 7): (1) 12.5 g glycerol, 1.0 g arginine-HCl, 1.0 g NaCl, 1.0 g K₂HPO₄, 0.5 g MgSO₄-7H₂O, 0.01 g Fe₂(SO₄)₃ -6H₂O, 1 mg FeSO₄-7H₂O, 1 mg MnCl₂-4H₂O, 1 mg ZnSO₄-7H₂O, 15.0 g agar ⁴⁵; (2) 20.0 g mannitol, 20.0 g soya flour, 20.0 g agar; (3) 5 ml glycerol, 10.0 g sucrose, 5.0 g beef extract, 5.0 g casamino acids; (4) ISP medium 4 ⁴⁶ with 10 mM N-acetylglucosamine. Strains were grown on these solid media for 10 days. Plates were frozen and then manually extruded through filter paper to obtain a crude extract. Crude sample was acidified to a final concentration of 0.1% formic acid and then loaded onto an Oasis HLB extraction column (Waters, Milford, MA). Following equilibration, columns were washed with 2 mL of water, followed by 2 mL of 10% ACN. Metabolites were eluted in 1.2 ml of 80% ACN and the eluate was evaporated to dryness. All flow-through was performed by gravity-flow, and a vacuum manifold was used when gravity was not sufficient. Microcentrifuge tubes were weighed empty and after drying to obtain a sample weight.

LC-MS Analysis

Extracts were resuspended in 5% ACN with 0.2% formic acid to a final concentration of 2 μg/μL. For LC-MS analysis, 40 μg of sample was loaded onto a 150 mm × 2.1 mm i.d., 2 μm particle size Kinetex C18 RPLC column (Phenomenex, Torrance, CA). Analysis was performed using an Agilent 1150 LC system (Agilent, Santa Clara, CA) equipped with a photodiode array and placed in-line with a Q-Exactive mass spectrometer (Thermo Fisher Scientific, Waltham, MA). Chromatography was performed at a flow rate of 200 μL/min using water/0.1% formic acid (solvent A) and acetonitrile/0.1% formic acid (solvent B) with the following gradient: time 0 min., 2% B; 35 min., 60% B; 54 min., 98% B. UV spectra were acquired at a rate of 1 Hz. For every spectrum, the three most intense apices were recorded. The mass spectrometer instrument settings were as follows: capillary temperature 275 °C, sheath gas 8 (arbitrary units), spray voltage 4.2 kV. Full MS spectra were acquired at 35,000 resolution for the mass range m/z 250 to 3750 for all samples. This resulted in an average scan rate of 6 Hz. Following each full MS scan, the top 5 most intense ions were selected for a dependent MS² scan. MS² was conducted using HCD with a collisional energy of 25%.

Software, Informatics and Statistical Treatment of Data

SIEVE software was used for chromatographic alignment, component detection, removal of background, and relative quantification. Chromatographic alignment was performed for all samples using an initial tile size of 500 frames. Following alignment, feature detection was performed using an initial mass tolerance of 10 ppm and a retention time window of 3.0 minutes. Analysis parameters were optimized for low m/z (250-500), mid m/z (500-900), and high m/z (900-3,750). A minimum normalized intensity of 5×10⁷, 1×10⁷, or 5×10⁶, was selected as the threshold for defining a peak as a feature in the low-, mid- and high-m/z regimes, respectively. Deisotoping and the summation of multiple adducts observed for a single species were performed to reduce data complexity. Following complexity reduction, a final list of components was output. For each component, a reconstructed ion chromatogram was created and the integrated intensity of the peak was calculated. An integrated intensity of 2×10⁶, 1×10⁶, or 7.5×10⁵ was selected as the cut-off value for a component to be considered present in a given sample in the low-, mid- and high-m/z regimes, respectively.

Automated Dereplication

All component mass values for each of the components were searched against an accurate mass database consisting of known bacterial metabolites using a mass tolerance of only 3 parts-per-million (ppm). The database of 9,817 known natural products was prepared using Antibase (http://wwwuser.gwdg.de/~hlaatsc/antibase.htm), Dictionary of Natural Products, as well as additional bacterial natural products found in the literature. Manual analysis of fragmentation spectra was performed to confirm intact mass-based identifications. In addition, UV data was compared to literature values when available to offer additional and orthogonal confirmation of compound identity (Supplementary Table 3).

MS/MS Verification Procedure

All compounds putatively identified via accurate intact mass were confirmed using accurate mass, tandem MS (MS²) data. To ensure that low-quality spectra were not included, MS² spectra containing less than 5 peaks at >1% relative abundance were excluded from analysis. Additionally, spectra containing more than 100 peaks at >1% abundance were included only if >20% of the peaks appeared in the higher m/z half of the spectrum. A large number of very low m/z fragment ions in the absence of high molecular weight fragments was found to be characteristic of over-fragmentation, and had the potential to lead to false identifications. Structures for putatively identified compounds were fragmented in silico using the software suite Mass Frontier (Thermo Fisher Scientific, Waltham, MA). Both general fragmentation rules and fragmentation library modes were used. The fragmentation library was composed of the HighChem ESI Positive 2008 library as well as annotated spectra from in-house compound libraries. General fragmentation rules were bypassed for library reactions. For general fragmentation rules, electron sharing and charge stabilization resonance reactions were allowed. Ionization, stabilization, and cleavage were allowed on aromatic systems. For compounds with a molecular mass less than 500 Da, 8 reaction steps were allowed with an upper limit of 10,000 reactions. For compounds with a molecular mass greater than 500 Da, 12 reaction steps were allowed with an upper limit of 20,000 reactions. Confirmation required that at least 4 of the top 5 most intense peaks were consistent with theoretical fragment ions to within 5 ppm. On average, 8.4 out of the top 10 most abundant peaks were consistent with theoretical fragment ions among all confirmed species in the dataset.

1

• 1.Bérdy J. Bioactive microbial metabolites. J. Antibiot. 2005;58:1–26. doi: 10.1038/ja.2005.1. [DOI] [PubMed] [Google Scholar]
• 2.Bérdy J. Thoughts and facts about antibiotics: Where we are now and where we are heading. J. Antibiot. 2012;65:385–395. doi: 10.1038/ja.2012.27. [DOI] [PubMed] [Google Scholar]
• 3.Bentley SD, et al. Complete genome sequence of the model actinomycete Streptomyces coelicolor A3(2) Nature. 2002;417:141–147. doi: 10.1038/417141a. [DOI] [PubMed] [Google Scholar]
• 4.Lautru S, Deeth RJ, Bailey LM, Challis GL. Discovery of a new peptide natural product by Streptomyces coelicolor genome mining. Nat. Chem. Biol. 2005;1:265–269. doi: 10.1038/nchembio731. [DOI] [PubMed] [Google Scholar]
• 5.Kersten RD, et al. A mass spectrometry–guided genome mining approach for natural product peptidogenomics. Nat. Chem. Biol. 2011;7:794–802. doi: 10.1038/nchembio.684. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 6.Ziemert N, et al. The natural product domain seeker NaPDoS: a phylogeny based bioinformatic tool to classify secondary metabolite gene diversity. PLoS ONE. 2012;7:e34064. doi: 10.1371/journal.pone.0034064. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 7.Medema MH, et al. antiSMASH: rapid identification, annotation and analysis of secondary metabolite biosynthesis gene clusters in bacterial and fungal genome sequences. Nucleic Acids Res. 2011;39:W339–W346. doi: 10.1093/nar/gkr466. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 8.Conway KR, Boddy CN. ClusterMine360: a database of microbial PKS/NRPS biosynthesis. Nucleic Acids Res. 2013;41:D402–D407. doi: 10.1093/nar/gks993. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 9.Diminic J, et al. Databases of the thiotemplate modular systems (CSDB) and their in silico recombinants (r-CSDB) J. Ind. Microbiol. Biotechnol. 2013;40:653–659. doi: 10.1007/s10295-013-1252-z. [DOI] [PubMed] [Google Scholar]
• 10.Yadav G, Gokhale RS, Mohanty D. SEARCHPKS: a program for detection and analysis of polyketide synthase domains. Nucleic Acids Res. 2003;31:3654–3658. doi: 10.1093/nar/gkg607. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 11.Tae H, Kong E-B, Park K. ASMPKS: an analysis system for modular polyketide synthases. BMC Bioinformatics. 2007;8:327. doi: 10.1186/1471-2105-8-327. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 12.Ichikawa N, et al. DoBISCUIT: a database of secondary metabolite biosynthetic gene clusters. Nucleic Acids Res. 2013;41:D408–D414. doi: 10.1093/nar/gks1177. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 13.Caboche S, et al. NORINE: a database of nonribosomal peptides. Nucleic Acids Res. 2008;36:D326–D331. doi: 10.1093/nar/gkm792. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 14.Kim J, Yi G-S. PKMiner: a database for exploring type II polyketide synthases. BMC Microbiol. 2012;12:169. doi: 10.1186/1471-2180-12-169. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 15.Fischbach M, Walsh C. Assembly-Line Enzymology for Polyketide and Nonribosomal Peptide Antibiotics: Logic, Machinery, and Mechanisms. Chem. Rev. 2006;106:3468–3564. doi: 10.1021/cr0503097. [DOI] [PubMed] [Google Scholar]
• 16.Raghupathy N, Durand D. Gene cluster statistics with gene families. Mol. Biol. Evol. 2009;26:957–968. doi: 10.1093/molbev/msp002. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 17.Wang X, et al. Identification and characterization of the actinomycin G gene cluster of Streptomyces iakyrus. Mol. Biosyst. 2013;9:1286–1289. doi: 10.1039/c3mb70081j. [DOI] [PubMed] [Google Scholar]
• 18.Colwell RK, et al. Models and estimators linking individual-based and sample-based rarefaction, extrapolation and comparison of assemblages. J. Plant Ecol. 2012;5:3–21. [Google Scholar]
• 19.Doroghazi JR, Metcalf WW. Comparative genomics of actinomycetes with a focus on natural product biosynthetic genes. BMC Genomics. 2013;14:611. doi: 10.1186/1471-2164-14-611. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 20.Jensen PR, Williams PG, Oh DC, Zeigler L, Fenical W. Species-Specific Secondary Metabolite Production in Marine Actinomycetes of the Genus Salinispora. Appl. Environ. Microbiol. 2006;73:1146–1152. doi: 10.1128/AEM.01891-06. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 21.Dunbar KL, Melby JO, Mitchell DA. YcaO domains use ATP to activate amide backbones during peptide cyclodehydrations. Nat. Chem. Biol. 2012;8:569–575. doi: 10.1038/nchembio.944. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 22.Charlop-Powers Z, Owen JG, Reddy BVB, Ternei MA, Brady SF. Chemical-biogeographic survey of secondary metabolism in soil. Proc. Natl. Acad. Sci. U. S. A. 2014;111:3757–3762. doi: 10.1073/pnas.1318021111. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 23.Bunge J, Willis A, Walsh F. Estimating the Number of Species in Microbial Diversity Studies. Annual Review of Statistics and Its Application. 2014;1:427–445. [Google Scholar]
• 24.Nguyen DD, et al. MS/MS networking guided analysis of molecule and gene cluster families. Proc. Natl. Acad. Sci. U. S. A. 2013;110:E2611–E2620. doi: 10.1073/pnas.1303471110. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 25.Cote R. In: ATCC Bacteria and Bacteriophages, 19th edn. Pienta P, Tang J, Cote R, editors. American Type Culture Collection; 1996. p. 484. [Google Scholar]
• 26.Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18:821–829. doi: 10.1101/gr.074492.107. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 27.Li R, et al. De novo assembly of human genomes with massively parallel short read sequencing. Genome Res. 2010;20:265–272. doi: 10.1101/gr.097261.109. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 28.Chaisson MJ, Brinza D, Pevzner PA. De novo fragment assembly with short mate-paired reads: Does the read length matter? Genome Res. 2009;19:336–346. doi: 10.1101/gr.079053.108. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 29.Margulies M, et al. Genome sequencing in microfabricated high-density picolitre reactors. Nature. 2005;437:376–380. doi: 10.1038/nature03959. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 30.Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W. Scaffolding preassembled contigs using SSPACE. Bioinformatics. 2011;27:578–579. doi: 10.1093/bioinformatics/btq683. [DOI] [PubMed] [Google Scholar]
• 31.Peng Y, Leung HC, Yiu S-M, Chin FY. IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics. 2012;28:1420–1428. doi: 10.1093/bioinformatics/bts174. [DOI] [PubMed] [Google Scholar]
• 32.Hyatt D, et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119. doi: 10.1186/1471-2105-11-119. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 33.Li L, Stoeckert CJ, Roos DS. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 2003;13:2178–2189. doi: 10.1101/gr.1224503. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 34.Kurtz S, et al. Versatile and open software for comparing large genomes. Genome Biol. 2004;5:R12. doi: 10.1186/gb-2004-5-2-r12. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 35.Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–2461. doi: 10.1093/bioinformatics/btq461. [DOI] [PubMed] [Google Scholar]
• 36.Ester M, Kriegel H-P, Sander J, Xu X. In: KDD. Simoudis Evangelos, Han Jiawei, Fayyad Usama., editors. Vol. 96. 1996. pp. 226–231. [Google Scholar]
• 37.Shannon P, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–2504. doi: 10.1101/gr.1239303. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 38.Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–1797. doi: 10.1093/nar/gkh340. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 39.Price MN, Dehal PS, Arkin AP. FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS ONE. 2010;5:e9490. doi: 10.1371/journal.pone.0009490. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 40.Zhang Q, Yu Y, Vélasquez JE, van der Donk WA. Evolution of lanthipeptide synthetases. Proc. Natl. Acad. Sci. U. S. A. 2012;109:18361–18366. doi: 10.1073/pnas.1210393109. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 41.Yutin N, Puigbò P, Koonin EV, Wolf YI. Phylogenomics of prokaryotic ribosomal proteins. PLoS ONE. 2012;7:e36972. doi: 10.1371/journal.pone.0036972. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 42.Larkin MA, et al. Clustal W and clustal X version 2.0. Bioinformatics. 2007;23:2947–2948. doi: 10.1093/bioinformatics/btm404. [DOI] [PubMed] [Google Scholar]
• 43.Huerta-Cepas J, Dopazo J, Gabaldón T. ETE: a python Environment for Tree Exploration. BMC Bioinformatics. 2010;11:24. doi: 10.1186/1471-2105-11-24. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 44.Krzywinski M, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19:1639–1645. doi: 10.1101/gr.092759.109. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 45.El-Nakeeb MA, Lechevalier HA. Selective isolation of aerobic actinomycetes. Appl. Microbiol. 1963;11:75–77. doi: 10.1128/am.11.2.75-77.1963. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 46.Smith SE, et al. Comparative Genomic and Phylogenetic Approaches to Characterize the Role of Genetic Recombination in Mycobacterial Evolution. PLoS ONE. 2012;7:e50070. doi: 10.1371/journal.pone.0050070. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials