Streptomyces ambofaciens ATCC 23877 was obtained from the American Type Culture Collection (ATCC). The identity of this strain was verified by sequence analysis using the 16S rRNA genes specific ‘universal’ primers 5′-CAGCAGCCGCGGTAATAC-3′ and 5′-CCGTCAATTCCTTTGAGTTT-3′ (Olsen et al., 1986). The strain was stored in 1 ml cryotubes at -80°C as frozen mycelium in yeast starch (YS) medium containing 15% glycerol at a biomass concentration of approximately 0.25 g dry cell weight (DCW)/ml, or at -20°C as spores in 20% glycerol at a titre of approximately 5 × 10⁸/ml. The composition (per liter) of the media used in this study for S. ambofaciens growth and manipulation is here reported. YS broth: 2 g yeast extract, 10 g soluble starch, (18 g agar in YS agar) (pH 7.0); SMII (medium for conjugation experiments): 15 g dextrose, 10 g soybean flour, 0.5 g MgSO₄ ⋅ 7 H₂O, 5 g CaCO₃, 15 g agar; 2XYT (medium for spore resuspension): 16 g tryptone, 10 g yeast extract, 5 g NaCl.

Escherichia coli strain DH5α [F^- Φ80d lacZΔM15 endA1 recA1 hsdR17 supE44 thi-1 λ^- gyrA96 λ (lacZYA-argF) U169] was used in cloning procedures. This strain was cultivated in Luria-Bertani broth supplemented with 50 μg/ml ampicillin when required. Conjugation experiments were performed with E. coli strain GM2929/pUB307::Tn7. Lennox Broth (LeB) was used to grow, wash and resuspend this strain. The composition (per liter) of the LeB is as follows: 10 g tryptone, 5 g yeast extract, 5 g NaCl, 1 g glucose.

Streptomyces ambofaciens total genomic DNA was obtained as described by (Kieser et al., 2000). Manipulation of E. coli DNA was performed as indicated by Sambrook and Russell (2001). PCR amplification conditions were 94°C for 5 min followed by 34 cycles of 94°C for 1 min, 65°C for 1 min and 72°C for 1 min, and finally 72°C for 10 min.

Overexpression of S. ambofaciens ccr (srm4) (SAM23877_ RS26790 in S. ambofaciens ATCC 23877 genome) coding for crotonyl-CoA carboxylase/reductase was achieved by cloning the PCR-amplified coding region into an ad hoc engineered expression vector, pNGEM/OriT/P, under the control of groES promoter. To construct pNGEM/OriT/P we utilized the E. coli–Streptomyces shuttle vector pN702GEM3. pN702GEM3 (Fernandez-Abalos et al., 2003) contains the pIJ702 origin of replication (for Streptomycetes), the bifunctional Neo/Kan resistance marker from Tn5 (for both E. coli and Streptomycetes), the pUC origin of replication (for E. coli) and a polylinker derived from pGEM3Zf(+) (Promega). This vector is replicative but not conjugative in Streptomycetes. To transform it into conjugative vector pNGEM/OriT a DNA fragment containing the oriT determinant was amplified by PCR from plasmid pTYM18 (Onaka et al., 2003) using the primer pair OriTHindIII-for and OriTHindIII-rev whose sequence is reported below (HindIII restriction sites are underlined):

The amplified fragment was treated with HindIII, and the resulting 930 bp DNA fragment was cloned into pN702GEM3 using the HindIII restriction site of the polylinker.

To construct the pNGEM/OriT/P plasmid, the groES gene promoter (P_groES) region was amplified by PCR using S. ambofaciens ATCC 23877 genomic DNA as template, and the primer pair groES-for and groES-rev (XbaI and NdeI restriction sites are underlined):

The amplified DNA fragment was digested with XbaI, and the resulting 325 bp DNA fragment was inserted into the polylinker of pNGEM/OriT, generating pNGEM/OriT/P.

The crotonyl-CoA carboxylase/reductase (CRR) encoding gene srm4 (locus tag: SAM23877_ RS26790) was amplified from S. ambofaciens ATCC 23877 genomic DNA, using the primer pair CCR-for and CCR-rev (NdeI and EcoRI restriction sites are underlined):

The amplified 1475 bp DNA fragment was digested with NdeI and EcoRI and then ligated to NdeI and EcoRI –cleaved pNGEM/OriT/P, obtaining the resulting plasmid pNGEM/OriT/P/CCR.

Plasmids pNGEM/OriT/P and pNGEM/OriT/P/CCR were introduced into the S. ambofaciens strain by conjugation with E. coli GM2929/pUB307::Tn7, as described previously Kieser et al. (2000). To allow plasmid selection, conjugation medium SMII was supplemented with kanamycin (25 μg/ml).

To overproduce and purify His6-tagged S. ambofaciens CCR (Srm4) in E. coli, the entire ccr (srm4) gene was amplified by PCR using the pNGEM/OriT/P/CCR construct as template, and using the primer pair CCR-for (above reported) and pETCCR-rev (HindIII restriction site is underlined):

The amplified DNA fragment was digested with NdeI and HindIII and then ligated to NdeI and HindIII –cleaved pET21b, obtaining the resulting plasmid pET/Srm4. pET/Srm4 was then introduced into E. coli BL21(DE3) competent cells.

Escherichia coli BL21(DE3) harboring pET/Srm4 was grown in 1 l of LB medium supplemented with ampicillin (10 μg ml-1) to A600 = 0.6 at 37°C and 250 r.p.m. Culture was induced by the addition of 0.4 mM IPTG and incubated for 16 h at 20°C and 180 r.p.m. Cells were harvested by centrifugation at 5,500 RCF for 10 min, re-suspended in lysis buffer (20 mM Tris–HCl pH 8.0, 200 mM KCl, 25 mM imidazole, 10% glycerol and 0.1 mg ml^-1 DNase I). Cells were lysed by passing the suspension through a French Press twice at 1000 psi, and the cell debris were removed by centrifugation at 10,000 × g for 20 min at 4°C. The clarified cell lysate was bound to 1 ml of Ni-IMAC resin (Bio-Rad Laboratories) for 1 h at 4°C. The resulting slurry was transferred to a gravity-flow column, and washed twice with 1 ml of lysis buffer. The recombinant protein His6-Srm4 was then eluted with 0.2 ml of elution buffer (20 mM Tris–HCl pH 8.0, 500 mM NaCl, 500 mM imidazole and 10% glycerol), and frozen for storage at -80°C.

Spiramycin production by S. ambofaciens cultures grown in YS broth was assessed by HPLC. At different time intervals, supernatants were filtered through Phenex-RC membrane (0.45 μm; Phenomenex). Four hundred microliter filtrated samples were lyophilized and resuspended in a mixture water-acetonitrile 70:30. The concentration of spiramycins was determined by HPLC (Beckman System Gold Programmable Solvent Module125) equipped with a UV detector (232 nm). The column was a reverse phase LiChrospher^®RP-8 HPLC Column (Supelco) and the mobile phase was a mixture of 100 mM phosphate buffer (pH 2.3) containing 50 mM sodium perchlorate and acetonitrile 70:30 v/v with a flow rate of 0.8 ml/min. Purified spiramycin (Sigma–Aldrich) freshly dissolved in phosphate buffer (pH 7) was used as standard. Two-tailed paired Student’s t-test was used to assess statistical significance of antibiotic titre differences between bacterial strains during the time course. Statistical significance was declared at p-value < 0.05.

Streptomyces ambofaciens cultures (100 ml) were harvested by centrifugation at 10,000 × g at 4°C for 20 min. Supernatants were discarded and mycelial pellets were washed twice in homogenization buffer (50 mM potassium phosphate pH 7.5, 1 mM EDTA, 1 mM dithioerythritol (DTE), 10% glycerol, 0.1 mM phenylmethylsulfonyl fluoride). After resuspension in homogenization buffer, mycelial extracts were prepared by mycelial disruption at 2000 psi using a French Press and subsequent centrifugation of the homogenates at 25,000 x g for 5 min at 4°C to remove debris. The supernatant was collected and used for the enzyme activity measurement. Protein concentrations were determined by using the colorimetric Bio-Rad colorimetric assays for protein determination with BSA as a standard. Protein extracts were analyzed by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) according to standard procedures.

The crotonyl-CoA carboxylase/reductase (CCR) activity was determined as described (Erb et al., 2007) in a reaction mixture containing 50 mM phosphate buffer (pH 7.5), 1 mM EDTA, 1 mM DTE, 1 mM crotonyl-CoA, 30 mM NaHCO₃ and 0.1–0.3 mg of mycelial extract proteins, or 0.10–0.2 μg of purified CCR protein. Upon equilibration at 30°C for 5 min, 6 mM NADPH was added to initiate the reaction. CCR activity was measured spectrophotometrically by following the decrease in absorbance at 340 nm (𝜀_NADPH = 6.22 mM^-1 cm^-1), which corresponds to the oxidation of the NADPH. Specifically, the decrease in the absorbance due to CCR-independent NADPH oxidation was excluded from the calculation of CCR activity, by performing an analogous reaction in the absence of crotonyl-CoA. A unit (U) of enzyme activity is the amount of enzyme catalyzing the oxidation of 1 μmol NADPH min^-1.

Total bacterial RNA was extracted from S. ambofaciens growing in YS medium at 28°C on a rotary shaker at 200 rpm, by using the RNeasy Mini kit according to the manufacturer’s instructions (QIAGEN). RNase-free DNase I was used to eliminate traces of DNA in the samples, in accordance with the instructions of the manufacturer (Promega, Madison, WI, USA).

Semi-quantitative analyses of the ccr-specific transcripts, normalized to 16S rRNA, were carried out by Real-Time qPCR. Total RNAs (1 μg) were reverse-transcribed by using random hexamer (2.5 μM) with Superscript RT (Invitrogen). About 0.1–1% of each RT reaction was used to run real-time PCR on a SmartCycler System (Cepheid) with SYBR^®Green JumpStart Taq ReadyMix (Sigma–Aldrich) and the oligonucleotide primer pairs 16Suniv-1/16S-univ-2 (specific for 16S rRNA gene) and ccr-rtF/ccr-rtR (specific for ccr). The primer sequences (and names) were: 5′-CAGCAGCCGCGGTAATAC-3′ (16Suniv-1); 5′-CCGTCAATTCCTTTGAGTTT-3′ (16Suniv-2); 5′-CGATCGTCACCTGCGCCTCCAC-3′ (ccr-rtF); 5′-CGAAGTGGGTGCCGACGATAC-3′ (ccr-rtR). Primers were synthesized as a service by MWG-Biotech AG Oligo Production. RT-PCR product lengths were 185 bp for 16Suniv-1/16S-univ-2 and 96 bp for ccr-rtF/ccr-rtR. Real-Time qPCR samples were run in triplicate. The Real-Time qPCR conditions were: 30 s at 94°C, 30 s at 60°C, 30 s at 72°C for 35 cycles; detection of PCR products was performed at 83°C. The relative transcript levels were calculated using the 2ΔΔct method (Livak and Schmittgen, 2001).

For RNAseq experiments ribosomal RNAs were depleted starting from 1 μg of total RNA from each of the time points and biological replicates by using the RiboZero Gram positive kit (Epicentre, Illumina), strand specific RNA-seq libraries were prepared by using the ScriptSeq^TM v2 RNAseq library preparation kit (Epicentre, Illumina) starting from 50 ng of previously rRNA depleted RNA from each biological replicate and for all the time points analyzed. Then each library was sequenced on a MiSeq Illumina sequencer and 76 bp reads were produced.

Bowtie 2 (v2.2.6) (Langmead and Salzberg, 2012) was used to align raw reads to S. ambofaciens ATCC 23877 genome (GCF_001267885.1). Multi-mapping reads alignment was forced using non-deterministic option to obtain a single mapping locus for the reads located on the terminal arm region (nt 1-202694/8101246-8303940), which is known to be duplicated. High quality reads were selected imposing the following criteria: mapping quality for uniquely mapping reads MAPQ greater than 30 and alignment score greater than -15; for multi-mapping reads the alignment score was set equal or greater than -15. These criteria define high quality reads only those showing a maximum of three variations (mismatches or short indels): the sum of penalty scores assigned to mismatches and short indels should not be greater than -15.

rRNA depletion, strand specificity and gene coverage were evaluated using BEDTools (v2.20.1) (Quinlan and Hall, 2010) and SAMtools (v0.1.19) (Li et al., 2009) to verify the library preparation and sequencing performances (see Supplementary Data Sheet 7 and Table S1). To obtain the raw read counts of protein coding genes on unique chromosomal regions and on the pSAM1 plasmid, strand specific reads overlapping for at least 50% of their length to the genes present in S. ambofaciens ATCC 23877 RefSeq annotation were considered in each sample. For protein coding genes located in duplicated regions it was impossible to determine the exact chromosomal origin of the reads. To overcome this problem we reported only one of the two homologs in the final list of genes and the relative raw count is the sum of all the reads aligning to both loci (strand specificity was not considered). The R (R Development Core Team, 2012) package DESeq2 (v1.6.1) (Love et al., 2014) was then used to normalize the counts basing on the median-of-ratio method to estimate the scaling factors.

Time course analysis was performed with maSigPro (Conesa et al., 2006; Nueda et al., 2014). To find the genes whose expression was significantly altered during the time-course, we considered only those expressed in at least five over eight samples. The regression parameters were left to default values (α = 0.05; Q = 0.05) in both the two steps of the analysis, in particular a forward step regression model was adopted in the second step. Significant genes lists were obtained by imposing R-square >0.7, in order to select genes correctly fitting the model as suggested by Conesa et al. (2006). Finally hierarchical clustering was applied to the expression values of the significant genes to individuate 9 time-dependent patterns of gene expression.

The enrichments of functional categories were calculated using Fisher test and BH correction for multiple testing, the categories with an adjusted p-value ≤ 0.05 were considered significantly enriched (see next paragraph for functional categories definition).

Raw data are publicly available at Sequence Reads Archive under accession number BioProject PRJNA342588.

The Refseq S. ambofaciens ATCC 23877 genome annotation was manually refined by using RAST (Overbeek et al., 2014) and antiSMASH (Weber et al., 2015): functional categories were defined according to RAST classification, and the genes encoding secondary metabolites, annotated by antiSMASH, were manually added to the category of secondary metabolism. This version of the genome annotation is reported in Supplementary Data Sheet 1. A draft metabolic reconstruction of S. ambofaciens ATCC 23877 was initially assembled using RAST annotation system with default parameters and then downloaded from the ModelSEED database (Henry et al., 2010). Afterward, this reconstruction was extensively curated by manual inspection following the main steps listed in Thiele and Palsson (2010). Several additional information sources were interrogated, including KEGG (Kanehisa and Goto, 2000), BRENDA (Scheer et al., 2011) and MetaCyc (Caspi et al., 2006). In addition, models of more/less related organisms (S. coelicolor and M. tuberculosis) were used as reference in a comparative genomics workflow for identifying potentially missing reactions. BLAST (Altschul et al., 1990) searches (adopting the Bidirectional Best Hit criterion) on the S. ambofaciens ATCC 23877 genome were carried out to confirm/reject the inclusion of further genes/reactions from related organisms to the S. ambofaciens model. The list of S. ambofaciens cellular transporters was refined probing the Transporter Classification Data Base (TCDB, Saier et al., 2014). Biosynthetic pathways leading to specific antibiotics known to be produced by S. ambofaciens were not initially included in the draft reconstruction. For this reason, extensive bibliographical data mining was carried out to include spiramycin (Richardson et al., 1990; Karray et al., 2007; Nguyen et al., 2010), stambomycin, antimycin, and congocidine into the reconstruction.

Overall, at the end of the manual refinement process, 112 reactions, 80 metabolites and 81 genes were added to the initial reconstruction.

To date, no detailed information on the biomass composition of S. ambofaciens ATCC23877 is available in scientific literature. Accordingly, the biomass forming reaction from the S. coelicolor A3(2) metabolic model (Borodina et al., 2005) was used to include a reaction accounting for biomass assembly in S. ambofaciens ATCC23877 reconstruction.

The FBA method was employed to simulate flux distribution in different conditions. Briefly, FBA is a constraint-based method relying on the representation of the biochemical system under investigation in the form of a stoichiometric matrix S (m × n), where m is the number of metabolites and n the number of reactions. FBA is based on the assumption of the cellular pseudo-steady state, according to which the net sum of all the production and consumption rates of each internal metabolite within a cell is considered to be zero. Under this assumption, the system can be described by the set of linear equations:

in which X_i is the concentration of metabolite i, S_ij is the stoichiometric coefficient of the i^th metabolite in the j^th reaction, v_j is the flux of the j^th reaction, N the entire set of metabolites and M the entire set of reactions.

Upper and lower bounds of flux through each reaction act as further constraints and are expressed as:

where l_b and u_b are the lower and upper limits for reaction j, respectively. Finally, FBA exploits linear programming to determine a feasible steady state flux vector that optimizes a given objective function (e.g., biomass production).

The reconstructed model was analyzed using COBRAToolbox-2.0 (Hyduke et al., 2011; Schellenberger et al., 2011) in MATLAB^®R2009b (Mathworks Inc.). Gurobi 5.6¹ and GLPK 4.32² solvers were used for computational simulations presented herein.

We calculated the growth dependencies between biomass and antibiotics biosynthesis using Pareto optimality as described in Heinken et al. (2013). We built Pareto optimum curves by (i) determining the minimal and maximal flux through biomass assembly reaction, (ii) fixing the flux through this reaction at different steps spanning minimal to maximal flux and (iii) maximizing the flux through each antibiotic synthesis reaction for each step. The procedure was repeated with biomass and antibiotics assembly reactions exchanged.

The FSEOF algorithm for flux scanning based on enforced objective flux was implemented for the selection of gene targets to be amplified as described in Choi et al. (2010). Briefly, initial fluxes were calculated by constraint-based flux analysis using the objective function of maximizing biomass formation and setting (i) upper and lower boundaries to exchange reactions as to simulate growth in complex medium, (see below) and (ii) the spiramycin production rate to the experimentally measured value of 0.0017 mmol^∗h^-1g^-1. Next, the objective function was set to spiramycin production and the theoretical maximum product formation rate was calculated. Afterward, FSEOF was performed under the objective function of maximizing cell growth while the spiramycin production rate is gradually enforced from the flux value used in the first simulation to the theoretical maximum production rate obtained with the second simulation (i.e., performing ten different model optimizations). Candidate target genes were then collected identifying those fluxes that increased upon the application of the enforced objective flux without changing the reaction’s direction. More specifically, we selected those reactions (and corresponding coding genes) showing an increased flux in at least one of the 10 simulations, without changing the reaction’s direction. Please refer to Choi et al. (2010) for a detailed and formal explanation of the FSEOF algorithm.

Mycelial growth and spiramycin production were measured during S. ambofaciens ATCC 23877 growth in shake flask batch cultivation using complex YS medium (Figure 1 and Supplementary Data Sheet 2). Biomass data displayed the typical two-stage growth curve of actinomycetes consisting of an initial stage of rapid growth (RG1) (12–40 h), a transition phase of stalled growth (T) (40–60 h), a secondary growth phase (RG2) (60–72 h) and a stationary (S) phase (72–120 h) (Figures 1A,B). The analysis of growth kinetics also showed that during RG1 specific growth rates (μ_max) were not constant over time but exhibited three distinct peaks at 12, 18, and 24 h, respectively (Figure 1B). Spiramycin production started at 18 h and reached the highest titre in the medium at 48 h (Figure 1A). Then its concentration decreased progressively till 78 h and afterward it remained nearly constant till 117 h (Figure 1A). The evaluation of specific production rates (q_p) led us to distinguish three distinct phases of spiramycin production (Figure 1B and Supplementary Data Sheet 2): a brief pulse of antibiotic production culminating at 18 h (phase I) during RG1, a robust and prolonged phase of production (between 27 and 51 h) (phase II) starting during RG1 and lasting over almost the entire T phase, and a third phase of production (between 81 and 108 h) (phase III) starting after RG2 and lasting over almost the entire S phase. Notably, production phases II and III seem to be preceded by growth rate peaks and are characterized by either growth rate decline (phase II) or no growth (phase III) demonstrating a strict interconnection between growth pulses and activation of spiramycin production (Supplementary Data Sheet 2).

To gain a better insight into the physiological and genetic mechanisms governing mycelial growth and antibiotic production, after RG1, global transcriptional analysis was carried out by RNA-Seq technology. For RNAseq experiments, samples were collected from shake flask batch cultures in YS medium at four different time points (T1 = 48 h, T2 = 72 h, T3 = 96 h, and T4 = 120 h) for RNA extraction. Two biological replicates were analyzed for each time point.

The maSigPro bioconductor package (Conesa et al., 2006) was adopted to analyze the RNA-Seq time series dataset starting from the normalized gene counts obtained after DESeq2 (Love et al., 2014) standard normalization (see Material and Methods). maSigPro was initially developed for the analysis of single and multi-series time course microarray experiments, and was recently adapted to deal also with Next Generation-Sequencing (NGS) series of data (Nueda et al., 2014). This package follows a two-step regression strategy to find genes with significant temporal expression changes and/or significant differences between experimental groups. Following the first step of regression the software individuated 3551 significant genes in our dataset and then, in a second step, it refined this number to a total of 1186 variable genes, 598 upregulated and 588 downregulated during the time course. Among them we found 30 upregulated genes coding for biosynthetic enzymes for secondary metabolites, most of which involved in butyrolactone, stambomycin, and alpomycin/kinamycin production; while among the downregulated genes 37 participate in spiramycin production. Moreover 34 genes were positioned on the duplicated arm region where the cluster of genes coding for alpomycin and kinamycin production is located, while 32 genes were expressed by the extrachromosomal plasmid pSAM1 (Supplementary Data Sheet 3).

The 1186 variable genes were finally subdivided into nine clusters according to their gene expression profiles (Supplementary Data Sheets 3, 7 and Figure S1) and median profile was then inferred for each cluster (Figure 2). Clusters 1, 4, 6, and 8 consist of genes whose expression levels increase during the time course, although with distinct patterns and basal expression levels. In particular the expression levels of genes in cluster 8 has high basal expression at T1 (48 h), rapidly increased from T1 (48 h) to T2 (72 h) and then remained roughly constant throughout T3 (96 h) and T4 (120 h). Expression levels of genes in cluster 4 also exhibited high basal expression at T1 (48 h), progressively increased from T1 to T3 and remained high at T4. These two clusters are significantly enriched in functional categories “cell division and cell cycle” (cluster 4), “cell wall and capsule” (clusters 4 and 8), “DNA metabolism” (clusters 4 and 8), and “sulfur metabolism” (cluster 8) (Supplementary Data Sheet 3). The enrichment in these categories is consistent with growth kinetics indicating a secondary growth phase (RG2) between 60 and 72 h (Figure 1).

Genes in clusters 1 and 6 showed low expression levels at T1 (48 h) and gradually reached their highest expression level during T4, although with different kinetics. Cluster 1 is enriched in the functional category “secondary metabolism” (Supplementary Data Sheet 3). In particular, it contains five type I polyketide synthase (PKS)-encoding genes involved in stambomycin biosynthesis, and several regulatory genes (belonging to TetR and SARP families) of the alpomycin/kinamycin gene cluster. Genes involved in butyrolactone biosynthesis also belong to this cluster (Supplementary Data Sheet 3). Cluster 6 is enriched in the functional category “fatty acid, lipid, isoprenoids,” but it also contains four type I PKS-encoding genes of the stambomycin gene cluster (Supplementary Data Sheet 3). This result indicates that all nine PKSs-encoding genes involved in the biosynthesis of the large (51-membered) glycosylated macrolide stambomycin (Laureti et al., 2011) achieve the highest expression levels in the late stationary phase (S). This finding is of interest for management of spiramycin production. In fact, the biosynthetic pathways leading to spiramycin and stambomycin compete with each other for common metabolic precursors (i.e., malonyl-CoA and methylmalonyl-CoA) during the assembly of the macrolide scaffold.

It is also noteworthy that the up-regulated clusters 1, 4, 6, and 8 contain more than 95% of the genes modulated on pSAM1 indicating a global gene expression increase of this extrachromosomal element during late growth.

Such a competition is possibly alleviated by differential temporal expression patterns of the genes involved in spiramycin and stambomycin biosynthesis (Supplementary Data Sheet 3) as suggested by maSigPro analysis that assigned 4 spiramycin biosynthetic genes to cluster 2, 16 genes to cluster 3, 5 genes to cluster 7 and 2 genes to cluster 9. Indeed, although with different kinetics, the spiramycin genes belonging to the clusters 2, 3, 5, 7, and 9 exhibited a general trend of down-regulation during the time course while the stambomycin genes increased their expression level during the time course (Figure 2 and Supplementary Data Sheet 7 and Figure S1). Consistently with growth kinetics (Figure 1), clusters 2 and 3 contains many genes coding for ribosomal proteins and translation factors, and genes involved in central carbon (and energy) metabolism. Indeed, all these genes are expected to be strongly down-regulated beyond the RG1. Cluster 2 is enriched in the functional categories “membrane transport,” “virulence, disease and defense,” while cluster 3 is enriched in genes involved in “secondary metabolism” (essentially, spiramycin biosynthesis). Cluster 7 exhibits an enrichment in the category “fatty acid, lipid, isoprenoids,” while clusters 5 and 9 exhibit no significant enrichment in any considered functional categories.

A S. ambofaciens ATCC 23877 draft reconstruction was obtained from ModelSeed as described in “Material and Methods.” As it is usually expected for automatically reconstructed metabolic networks, this draft model failed in predicting S. ambofaciens growth, both on known defined media and on an arbitrarily complex one. Thus, the model underwent a manual refinement process which led to the inclusion of 112 reactions, 80 metabolites and 81 genes, the re-formulation of the biomass assembly reaction (following the one available for S. coelicolor) and the inclusion of gap-filling reactions (i.e., not gene-encoded) to reconcile S. ambofaciens growth on defined minimum media (see below). The obtained reconstruction was named iMF1244 according to the standard nomenclature and its main features are reported in Table 1. The SBML-formatted version of iMF1244 is available as Supplementary Data Sheet 4.

The reliability of the reconstructed S. ambofaciens model was inferred by comparing the predicted growth rates on a set of different minimal media against a set of experimentally determined ones available from previous studies (Lounès et al., 1995, 1996a,b; Colombie et al., 2005). Accordingly, in silico minimal growth media were defined using exchange reactions present in the model, and biomass optimization was selected as the model objective function (O.F.). More in detail, lower bounds of exchange reactions accounting for all the salts present in the experimental medium were set to -10 mmol g^-1 × h^-1, in order to mimic non-limiting conditions. The carbon source of each available minimum growth medium was then chosen as the unique carbon source of this in silico medium and its uptake flux was set to the corresponding experimentally determined value. The compositions of all the simulated in silico media are reported in Supplementary Data Sheet 7 and Table S2. The growth rates that were inferred by the in silico simulation were then compared to those experimentally determined, revealing an overall good agreement (Pearson’s product-moment correlation 0.97, p-value 0.02, Figure 3). This indicates that the model should be able to reproduce a realistic picture of the S. ambofaciens ATCC 23877 metabolic landscape.

Afterward, we tested the response of the model in a simulated complex medium. The assimilation of yeast extract (YE) amino acids and its effect on the production of spiramycin in culture of S. ambofaciens has been extensively analyzed by Benslimane et al. (1995). In this work the authors showed that in a YE medium a specific growth rate (μ_max) of 0.15 h^-1 was soon established after setting up the culture. The specific growth rate then decreased sharply to reach zero value at 24 h and rose up again when a second growth phase with a μ_max of 0.012 h^-1 was observed. During this second growth phase, spiramycin biosynthesis was observed, with a specific production rate (q_p) of 1.46 mg^-1h^-1g^-1 (corresponding to 0.0017 mmol^∗h^-1g^-1). The growth on YE of S. ambofaciens was characterized by the consumption of all the amino acids available in the medium (Benslimane et al., 1995). Accordingly, to account for such behavior, we added to the model all the 20 amino acids exchange and transport reactions. Since no detailed information is available concerning the single amino acids uptake rates in the second growth phase on YE, each uptake flux was arbitrary set to the value of 0.015 mmol^∗h^-1g^-1. In such conditions, setting the growth rate exactly to the measured value of 0.012 h^-1, the model predicts a spiramycin production rate of 0.0016 mmol^∗h^-1g^-1, which resembles the one experimentally determined (0.0017 mmol^∗h^-1g^-1).

In a typical scenario, the metabolic pathways for the production of secondary products (such as antibiotics) should compete for the pool of available metabolic intermediates with all the other cellular biosynthetic processes (e.g., biomass formation and/or homeostasis maintenance). This is typically accounted for by the so-called Pareto front between the two objectives (Marler and Arora, 2004). Pareto front refers to the set of resources allocation that is Pareto efficient, i.e., the situation in which it is impossible to make any one objective better off without making at least one individually worse. In other words, a pair of objective functions (such as a biomass and antibiotic assembly reactions) can be balanced to find the set of optimal solutions for which one objective can only be improved at the expense of the other. Recently, it has been shown that cells’ metabolism operates close to the Pareto-optimal (Schuetz et al., 2012) and that the calculation of Pareto fronts represents a technique that is more effective than weighting the objectives (Gevorgyan et al., 2011), especially when genetic algorithms are used to approximate the front (Boghigian et al., 2010; Angione and Lió, 2015).

Accordingly, we computed (see Materials and Methods for details) Pareto optimality for each of the four antibiotics that were included in the reconstruction. As shown in Supplementary Data Sheet 7 and Figure S2 the maximum theoretical yields were observed in the condition of μ = 0 and were 0.0052, 0.0105, 0.00012 and 0.00086 mmol^∗h^-1g^-1 for spiramycin, antimycin, stambomycin and congocidine, respectively. Afterward, the model (correctly) predicts a decrease in the production rates of each antibiotic as a consequence of an increase of the growth rate. These results, together with the consistence of growth rates prediction, support the iMF1244 as being a reliable reconstruction of the metabolism of this S. ambofaciens ATCC 23877.

Here we used our reconstruction to provide a global description of S. ambofaciens metabolism in the analyzed time points, by integrating gene expression data with constraint-based metabolic modeling. In general, many methods to achieve this task have been developed to date (Machado and Herrgard, 2014). However, a recent work (Kim et al., 2016) showed that, among them, iMAT (Shlomi et al., 2008) led to more accurate predictions in the analysis of Streptomyces coelicolor metabolic model. Briefly, iMAT uses gene expression values to divide reactions into two groups: highly and lowly expressed. The algorithm next seeks to find the flux distribution that maximizes the consistency with this classification. Accordingly, we first computed the first and third quartiles for RPKM values in each time point to define down- and up-regulated genes, respectively (obtained values are reported in Supplementary Data Sheet 7 and Table S3). Next, we performed four different FBA optimizations (one for each of time points analyzed) using iMAT to constrain fluxes across reactions according to gene expression values. For each optimization the flux through biomass assembly and spiramycin production reactions was set to the values experimentally determined (see Figure 1 and Supplementary Data Sheet 7 and Table S4) and the resulting fluxes distribution analyzed. To account for the cellular maintenance cost, an arbitrary cellular growth rate of 1e^-10 h^-1 was set in those cases in which no growth had been experimentally observed. Results obtained, however, were not affected when this arbitrary cellular growth rate was replaced by a simulated ATP maintenance cost. In this case a distribution of overlapping fluxes was observed. Figure 4A shows the flux across each reaction for every time step and the number of flux-carrying reactions.

Our model predicts that, on average, 392 reactions are active across all the time points (standard deviation, SD = 63.4). Accordingly, around 1000 reactions are predicted to be inactive in each time point (73% of the model reactions, on average). This is in agreement with the fact that the bacterium is growing in a complex medium and thus can find in the surrounding environment most of its cellular building blocks. Interestingly, the second time point is the one displaying the highest number of flux-carrying reactions (475), in agreement with the fact that, unlike the other time points, S. ambofaciens is actively producing biomass in this moment of the growth phase. After this time point, S. ambofaciens experiences a decreasing trend in the number of flux-carrying reactions (Figure 4) reflecting with both the reduced growth rate and lower spiramycin productivity observed in vivo in last phases of the growth.

From a functional viewpoint, our modeling predicts that the set of active pathways is largely shared by the four time points (Figure 4B). Indeed, a similar scheme of reactions is apparently carrying flux across all the time points, including the main central routes of S. ambofaciens metabolism (e.g., purine and pyrimidine metabolism, glycolysis, TCA cycle). This suggests the absence of a massive metabolic reprogramming in the four time points analyzed, i.e., the final part of S. ambofaciens growth. The most evident exception to this general trend is represented by the second time point, the one also showing the highest number of active reactions. Our functional annotation of reactions predicted to be carrying flux suggests that this difference is due to the activation of four main pathways i.e., amino-sugar metabolism, steroid biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, riboflavin metabolism (Figure 4B). Thus, according to our simulations, these pathways (or part thereof) play a role when the cell is trying to achieve two objective functions simultaneously, i.e., biomass and antibiotic synthesis. The activation of riboflavin metabolism and a consequent increase in FAD cofactor availability, for example, is consistent with an increased cellular demand of acetyl-CoA (a substrate for the synthesis of many spiramycin intermediates) that, in turn, can be obtained through the FAD-dependent degradation of fatty acids. Similarly, the use of amino sugar metabolism might be propaedeutic to the recycling of nucleotide hexoses (e.g., glucose-1P), valuable precursors for the sugar moieties of the spiramycin molecule.

All these data show how changes in gene expression may redirect the cellular metabolic fluxes within the cell, leading to activate specific routes that, in turn, may represent key nodes in the utilization of intermediates that are common to multiple cellular objectives (biomass formation and spiramycin synthesis in this case). The hypothetical scenario depicted here indicates both FAD production and nucleotide sugars recycling as key factors in allowing the cell to achieve the imposed double optimization.

We next exploited the iMF1244 metabolic reconstruction for the genome-wide identification of gene amplification targets potentially allowing spiramycin overproduction in S. ambofaciens. To this purpose we used the FSEOF approach described in Choi et al. (2010), as reported in Section “Materials and Methods.” From the set of identified potential targets, we excluded from further consideration those reactions that (i) were added to the model for gap-filling, (ii) were involved in the import (or export) of nutrients into the system for modeling purposes (i.e., exchange reactions) and (iii) were either encoded by multiple genes (enzymatic complexes) or, at least, two alternative genes (paralogs) in the genome. In all these cases, no sufficiently detailed clues on the gene(s) to be amplified for spiramycin overproduction could be derived. Following these selection criteria, our approach led to the identification of 65 potential reactions (corresponding to 54 genes) whose increased flux may lead to an increased flow through spiramycin synthesis pathway. The complete list of the reactions identified is reported in Supplementary Data Sheet 5.

A large fraction of this set of reactions is involved in the synthesis of precursors of the macrolactone backbone of spiramycin. These include, for example, glycerone phosphate (produced by sn-glycerol-3-phosphate: NAD+ 2-oxidoreductase) involved in the biosynthesis of methoxymalonyl-ACP, crotonoyl-CoA (the product of butanoyl-CoA: oxygen 2-oxidoreductase) precursor of ethylmalonyl-CoA and acetyl-CoA, the substrate of 3-oxopropanoate: NAD+ oxidoreductase for the generation of malonyl-CoA. Reactions leading to the production of intermediates necessary for the synthesis of the deoxyhexose sugars embedded in the spiramycin molecule (e.g., alpha-D-glucose 1-phosphate 1,6-phosphomutase responsible for the conversion of D-glucose 1-phosphate into D-glucose 6-phosphate) were also identified.

Among the proposed targets, we selected the reaction carried out by crotonyl-CoA carboxylase/reductase and presumably involved in the synthesis of ethylmalonyl-CoA. The spiramycin macrolactone backbone is synthesized by a PKS that, apart from the most common polyketide precursors, malonyl-CoA and methylmalonyl-CoA, was predicted to incorporate the less common precursor, ethylmalonyl-CoA, and the functionalized extender unit, methoxymalonyl-CoA (Kuhstoss et al., 1996; Karray et al., 2007). In particular, the ethylmalonyl-CoA is used as building block to assemble only a few known polyketides, including the polyether antibiotic monensin A, whose biosynthetic pathway was the first established pathway to require ethylmalonyl-CoA (Liu and Reynolds, 1999).

Although there is more than one pathway to ethylmalonyl-CoA, the pathway involving the crotonyl-CoA carboxylase/reductase (CCR) (Erb et al., 2007) (commonly referred to as “the ethylmalonyl-CoA pathway”) appears to be the dominant source of this precursor for polyketide biosynthesis (Wilson and Moore, 2012). This pathway was originally described in purple non-sulfur bacterium Rhodobacter sphaeroides and then discovered also in methylotrophic bacteria such as Methylobacterium extorquens and in Streptomyces spp. In the ethylmalonyl-CoA pathway a C4-compound, acetoacetyl-CoA, derived from two molecules of acetyl-CoA, is converted to the C5-compound mesaconyl-CoA (2-methylfumaryl-CoA) that is then transformed to (2R,3S)-β-methylmalyl-CoA by hydration and is finally cleaved to glyoxylate and propionyl-CoA. The key enzyme of the pathway is CCR, which simultaneously carboxylates and reduces the C4-compound crotonyl-CoA, forming the C5-compound (2S)-ethylmalonyl-CoA (herein indicated as “ethylmalonyl-CoA”). As illustrated in Supplementary Data Sheet 7 and Figure S3, the ethylmalonyl-CoA and its direct metabolic precursor crotonyl-CoA lies at the crossroad between key metabolic pathways devoted to energy metabolism and energy storage, lipid and amino acid metabolism and secondary metabolism.

In S. ambofaciens the ethylmalonyl-CoA is used as a building block to synthesize the polyketide backbone of spiramycin and antimycin, and the full set of genes necessary for a functioning ethylmalonyl-CoA pathway is clustered together in the S. ambofaciens chromosome not very distant from the spiramycin biosynthetic cluster. The cluster comprises the ORFs: SAM23877_RS28800 coding for methylsuccinyl-CoA dehydrogenase (msd gene); SAM23877_RS28805 coding for mesaconyl-CoA hydratase (mcd gene); SAM23877_RS28810 coding for L-malyl-CoA/beta-methylmalyl-CoA lyase (mclA gene); SAM23877_RS28815 encoding ethylmalonyl-CoA mutase (ecm gene); SAM23877_RS28820 encoding the CCR (ccr gene); SAM23877_RS28825 coding for a TetR family transcriptional regulator; SAM23877_RS28830 coding for 3-hydroxybutyryl-CoA dehydrogenase (Supplementary Data Sheet 7 and Figure S3).

As well as in other polyketides gene clusters, a secondary copy of the CCR-encoding gene is located in both spiramycin (SAM23877_RS26790, srm4) and antimycin (SAM23877_RS01780) biosynthetic gene clusters. RNAseq data showed that under tested experimental conditions all three CCR-encoding genes were poorly expressed, the gene RS26790 has an expression level higher than the gene RS28820 which is more expressed than the gene RS01780 (RS26790 > RS28820 > RS01780), among them only the gene SAM23877_RS26790 (srm4) shows an expression level above the baseline threshold (RPKM > 10) in all the time points analyzed, thus allowing us to consider it effectively expressed along the growth curve (Table 2). Although the CCR enzymatic activity of the ccr paralog gene products has not been proven so far, their role is likely to synchronize ethylmalonyl-CoA synthesis from acetoacetyl-CoA with polyketide assembly, and/or to increase intracellular ethylmalonyl-CoA levels in order to ensure adequate levels of the precursor during polyketide biosynthesis. Therefore, the carbon flow through the ethylmalonyl-CoA pathway was consistently associated with spiramycin productivity by our metabolic model.

To validate this prediction, we planned to over-express the ccr paralog srm4 of the spiramycin gene cluster in S. ambofaciens. However, we first verified that srm4 actually codes for an enzyme with CCR activity. To this purpose, the His6-tagged S. ambofaciens CCR (Srm4) protein was overproduced in E. coli BL21(DE3) strain and purified (Supplementary Data Sheet 7 and Figure S4). Then its CCR enzymatic activity was determined. The results demonstrated an enzymatic activity of 449 ± 11 U mg-1. To over-express ccr (srm4), the gene was placed under the control of the PgroES promoter (Luo et al., 2015) in pN702GEM3, a multicopy E. coli-Streptomyces shuttle plasmid vector (Supplementary Data Sheet 7 and Figure S5). The pNGEM/OriT/P/CCR or the control plasmid pNGEM/OriT/P were then introduced into S. ambofaciens ATCC 23877 by intergeneric conjugation with E. coli, and mycelial growth (Figure 5A), spiramycin production (Figure 5B) ccr (srm4) mRNA levels (Figure 5C), and CCR enzymatic activity in crude cell extracts (Figure 5D) were analyzed in the trans-conjugant strains grown in YS medium. With respect to growth kinetics, compared to control strain (harboring pNGEM/OriT/P) the ccr-over-expressing strain exhibited major differences during the RG1 with about 14.6% lower final biomass values but similar growth rates (0.0697 vs. 0.0622 between 24 and 48 h) (Figure 5A and Supplementary Data Sheet 6b).

Strain manipulation did not alter significantly either soluble protein content (expressed as mg protein/g DCW) or electrophoretic pattern as determined, respectively, by Bio-Rad colorimetric assays and SDS-PAGE (Supplementary Data Sheet 7 and Figure S6). In each time point the ccr (srm4) mRNA expression levels of the strain over-expressing ccr (pNGEM/OriT/P/CCR), were compared at the same time both with the expression levels of the wild type and of the strain with control plasmid (pNGEM/OriT/P). As shown in Figure 5C the strain over-expressing ccr (srm4) exhibited an average maximum increase in the mRNA levels of 2–3-fold when compared respectively with wild type or with the control plasmid strain, thus indicating a moderate gene overexpression at each time point (Figure 5C). It is noteworthy that during the time course, in all the three strains ccr (srm4) mRNA levels reached the highest values of expression at 48 h. So in the light of these results we can infer that the percentage of the overexpressed Srm4 protein in the total soluble S. ambofaciens proteome is therefore almost negligible.

In the ccr-over-expressing strain spiramycin productivity was enhanced during the RG1 and T phases. Compared to control strain spiramycin specific production rates (q_p) were increased by a factor of 3.63 at 24 h (due to early triggering of antibiotic biosynthesis) and 1.52 at 48 h (Supplementary Data Sheet 6b), but antibiotic titres were only modestly increased by a factor of 1.503 at 48 h (i.e., an average increase of 50.3% with a standard deviation of 36.5%, p-value = 0.01078), when spiramycin titres reached the highest values during the time course. Antibiotic titres were instead significantly increased by a factor of 4.336 at 72 h (i.e., an average increase of 333.6% with a standard deviation of 235.6%, p-value = 0.00024) (Figure 5B and Supplementary Data Sheet 6a). In contrast, no statistically significant differences were observed during the late phase of antibiotic production (during RG2 and S phases).

The achieved increase was rather limited when compared to that obtained with similar metabolic engineering strategy that resulted in about 2- and 1.8-fold increase in actinorhodin production by S. coelicolorA3(2) strains overexpressing, respectively, the ribulose 5-phosphate 3-epimerase and the NADP-dependent malic enzyme (Kim et al., 2016). However, the increase was comparable to those obtained in other streptomycetes by genetic modification of single targets (Huang et al., 2012, 2013). For instance, in Streptomyces tsukubaensis 35, 40, or 50% improvement in FK506 (a 23-membered polyketide macrolide also known as tacrolimus) production could be achieved, respectively, by overexpression of zwf2 (coding for glucose-6-phosphate dehydrogenase), accA2 (coding for acetyl-CoA carboxylase) or dahp (coding for 3-deoxy-D-arabino-heptulosonate-7-phosphate synthase). However, the combined effect of gdhA-deletion and dahp-, accA2-, zwf2-overexpression enhanced the concentration of FK506 up to 1.47-fold in fed-batch fermentations (Huang et al., 2013). In another streptomycete, such as for example Streptomyces roseosporus, co-overexpression of zwf2, dptI, and dptJ genes resulted in a 34.4% higher daptomycin concentration compared with the parental strain (Huang et al., 2012). These comparisons indicate that the achieved increase may not be irrelevant, and that spiramycin production may be further enhanced by the combined effect of ccr-overexpression, and genetic manipulation of other identified targets and/or improvement of fermentation conditions. Preliminary data indicate very long-chain acyl-CoA dehydrogenases, and several enzymes involved in glyoxylate cycle and glycine, serine and threonine metabolism as further good targets for strain improvement.

The results of biochemical assay with crude mycelial extracts demonstrated that the increase in spiramycin productivity paralleled with an increase in CCR activity (Figure 5D). In particular, at 48 h of growth, the enzyme activity in the ccr-over-expressing strain was about two-fold higher than in control strain (harboring pNGEM/OriT/P). It should be noted also that the order of magnitude of CCR activity that was determined in the protein extracts from S. ambofaciens is similar to those assessed in protein extracts from Rhodobacter sphaeroides (0.7 U mg^-1), Methylobacterium extorquens (0.8 U mg^-1) and Streptomyces coelicolor A3(2) (0.4 U mg^-1) (Erb et al., 2007). Altogether these results validate the metabolic model providing evidence that spiramycin production may be increased by genetic manipulation of the ethylmalonyl-CoA metabolic node.