Transcriptome and translatome profiles of Streptomyces species in different growth phases

Streptomyces are efficient producers of various bioactive compounds, which are mostly synthesized by their secondary metabolite biosynthetic gene clusters (smBGCs). The smBGCs are tightly controlled by complex regulatory systems at transcriptional and translational levels to effectively utilize precursors that are supplied by primary metabolism. Thus, dynamic changes in gene expression in response to cellular status at both the transcriptional and translational levels should be elucidated to directly reflect protein levels, rapid downstream responses, and cellular energy costs. In this study, RNA-Seq and ribosome profiling were performed for five industrially important Streptomyces species at different growth phases, for the deep sequencing of total mRNA, and only those mRNA fragments that are protected by translating ribosomes, respectively. Herein, 12.0 to 763.8 million raw reads were sufficiently obtained with high quality of more than 80% for the Phred score Q30 and high reproducibility. These data provide a comprehensive understanding of the transcriptional and translational landscape across the Streptomyces species and contribute to facilitating the rational engineering of secondary metabolite production.

Bacteria can fine-tune gene expression both at the transcriptional and translational levels 16,17 . For example, Escherichia coli proteome analysis revealed that only approximately half of protein abundance is determined by transcriptional regulation, which indicates the existence of various post-transcriptional regulation 18 . In this regard, deciphering translational dynamics is important to understanding post-transcriptional regulations that are closely related to cellular protein levels 19 . Recently, ribosome profiling has been used to measure translational levels by deep sequencing of the ribosome-protected mRNA fragments (RPFs) at the position of the translating ribosome 20 . Several ribosome profiling studies in Streptomyces have been reported by our research group for S. coelicolor, S. clavuligerus, and S. lividans, which revealed translational buffering of secondary metabolism-related genes at a later growth phase and that translational abundance is more consistently maintained than transcript abundance 11,21,22 . Translational regulations are advantageous for the tight control of secondary metabolite biosynthesis, as translation requires the highest energy costs among all cellular reactions 19 . Moreover, the expression of smBGC-associated genes can be more rapidly regulated at the translational level than at the transcriptional level in response to dynamic environmental changes 23 . Given the dynamic relationship between transcription and translation, as exhibited by translational buffering 11 , integrative analysis at both levels should unravel complex regulations in Streptomyces. However, transcriptomic and translatomic data have covered only a small portion of approximately 350 reported Streptomyces genomes, which have not been systematically validated at the multi-species level.
In this study, we provide RNA-Seq and ribosome profiling data of five Streptomyces species at four different growth phases, followed by validation of the read quality. The species were S. avermitilis MA-4680, S. clavuligerus ATCC27064, S. lividans TK24, S. venezuelae ATCC15439, and S. tsukubaensis NRRL 18488, which are industrial strains that produce antifungal avermectin, β-lactamase inhibitor clavulanic acid, and immunosuppressant FK506, respectively [24][25][26] . S. lividans and S. venezuelae were characterized by their fast growth and ease of genetic manipulation, and have been employed as heterologous expression hosts [27][28][29][30] . An overview of the preparation of transcriptomic and translatomic data is illustrated in Fig. 1. A total of 12 to 83.5 million raw reads for RNA-Seq and 113 to 763.8 million raw reads for ribosome profiling were obtained. Although the RNA-Seq and ribosome profiling data of two species (S. clavuligerus 21 and S. lividans 22 ) among the five species were already reported in previous studies by our research group, this study provided a uniformly processed and mapped dataset of all five species. This facilitates the efficiency of the comparative transcriptome and translatome analysis at multi-time points between multi-species. Further, understanding the transcriptional and translational regulatory mechanisms and developing regulatory synthetic parts, such as promoters, ribosome-binding sequences, 5′ untranslated regions, and terminators 4 from the dataset allows rational genome engineering for efficient secondary metabolite production by Streptomyces 11 .

Methods
Strains and cell growth. Streptomyces strains were inoculated from their 20% glycerol stock of spores into 50 mL of R5− liquid medium with 8 g of glass beads (3 ± 0.3 mm diameter) in a 250 mL baffled flask, grown at 30 °C, and pre-cultured at 250 rpm. The R5− liquid medium consists of 103 g L −1 sucrose, 0.25 g L −1 K 2 SO 4 , 10.12 g L −1 MgCl 2 •6H 2 O, 10 g L −1 glucose, 0.1 g L −1 casamino acids, 5 g L −1 yeast extract, 5.73 g L −1 TES (pH 7.2), 0.08 mg L −1 ZnCl 2 , 0.4 mg L −1 FeCl 3 •6H 2 O, 0.02 mg L −1 CuCl 2 •2H 2 O, 0.02 mg L −1 MnCl 2 •4H 2 O, 0.02 mg L −1 Na 2 B 4 O 7 •10H 2 O, 0.02 mg L −1 (NH 4 ) 6 Mo 7 O 24 •4H 2 O, and 0.28 g L −1 NaOH. The grown mycelium was inoculated to fresh R5− medium with an initial optical density of 0.05 at 600 nm for the main culture as biological duplicates and grown under the previously mentioned conditions. The cells were sampled at four different time points based on the growth profile of each strain, as follows: early-exponential (E), transition (T), late-exponential (L), and stationary (S) phases. The E, T, L, and S time points were 13, 17, 19.5, and 33.5 h for S. avermitilis, 26, 80, 105.5, and 125 h for S. clavuligerus, 9.5, 14, 16, and 20 h for S. lividans, 12.5, 24.5, 30.5, and 48.5 h for S. venezuelae, and 15, 18.5, 28, and 48 h for S. tsukubaensis after inoculation, respectively ( Fig. 1c). At the sampling time points for the ribosome profiling samples, thiostrepton (Sigma-Aldrich, St. Louis, MO, USA) was added to the cultures to a final concentration of 20 μM to compartment the translating ribosomes on the mRNA, which is a highly sensitive drug for Streptomyces compared to chloramphenicol or other drugs 31,32 . The cultures were then incubated for 5 min at 30 °C, and subsequently harvested for the construction of ribosome profiling libraries.

RNA-Seq library preparation and high-throughput sequencing.
The overview of the library construction of RNA-Seq is illustrated in Fig. 1a. The harvested cells were washed with polysome buffer (20 mM Tris-HCl, pH 7.5; 140 mM NaCl and 5 mM MgCl 2 ), and then resuspended with 500 μL lysis buffer (0.3 M sodium acetate, pH 5.2; 10 mM ethylenediaminetetraacetic acid and 1% Triton X-100). The resuspended cells were frozen with liquid nitrogen and grounded using a mortar and pestle. The ground mycelium was thawed and centrifuged at 4 °C for 10 min at 16,000 × g. The supernatant was collected and stored at −80 °C. Following the preparation of lysates from four growth phases as biological duplicates, the lysates were mixed with a solution of phenol:chloroform:isoamyl alcohol (25:24:1, v/v), and the mixtures were separated by centrifugation. DNA in the extracted RNA samples were removed by treatment with 2 μL DNase I (NEB, Ipswich, MA, USA), 5 μL 10 × DNase I buffer, and 1 μL SUPERase-In RNase Inhibitor (Thermo Scientific, Waltham, MA, USA). Lastly, the DNase I-treated RNA samples were purified using phenol:chloroform:isoamyl alcohol (25:24:1, v/v) and ethanol precipitation. To eliminate rRNAs in the recovered RNA samples, the Ribo-Zero rRNA Removal Kit for Bacteria (Epicentre, Madison, WI, USA) was used according to the manufacturer's instructions. The quality of rRNA-depleted RNA samples was checked using 2% agarose gel electrophoresis. The suitable RNA samples were then used to construct RNA sequencing libraries using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA). The size distributions of the final libraries were checked using the Agilent 2200 TapeStation System (Agilent, Santa Clara, CA, USA). The constructed libraries were sequenced on the HiSeq. 2500 platform using either a 100-bp (S. lividans, S. avermitilis, S. clavuligerus, and S. venezuelae) or 50-bp (S. tsukubaensis) single-end read recipe (Fig. 1a).
Ribosome profiling library preparation and high-throughput sequencing. An overview on the library construction of ribosome profiling is illustrated in Fig. 1a 21 . The mycelium that was treated by thiostrepton was collected by centrifugation at 4 °C for 10 min at 3,000 × g, and the cell pellet was washed with 2 mL of polysome buffer that was composed of 20 mM Tris-HCl (pH 7.5), 140 mM NaCl, and 5 mM MgCl 2 with 20 μM thiostrepton. The washed pellet was re-suspended in 1 mL of lysis buffer composed of 950 μL of polysome buffer www.nature.com/scientificdata www.nature.com/scientificdata/ and 50 μL of 20% Triton X-100 with 20 μM thiostrepton. The resuspended cells were dripped into a mortar filled with liquid nitrogen and then grounded with a pestle. The cell debris was removed by centrifugation at 4 °C for 5 min at 3,000 × g. The supernatant was further clarified and collected by centrifugation at 4 °C for 10 min at 16,000 × g. To digest RNA in the lysate (containing 50 μg RNA), the S. avermitilis and S. tsukubaensis samples were treated with 750 U of RNase I (Invitrogen, Waltharn, MA, USA) at 37 °C for 45 min, and the remaining strains were treated with 400 U of Micrococcal Nuclease (MNase) (NEB), 20 μl of 10× MNase buffer, and 2 μl of 100× Bovine Serum Albumin (BSA) (NEB) at 37 °C for 2 h. The samples were then loaded onto Illustra MicroSpin S-400 HR Columns (GE Healthcare Life Sciences, Marlborough, MA, USA) that were previously washed three times with 500 μL of washing buffer, which was composed of 50 mM Tris-HCl (pH 8.0), 250 mM NaCl, 50 mM MgCl 2 , 25 mM EGTA, and 1% Triton X-100. The column was centrifuged at 4 °C for 2 min at 400 × g, and the flow-through was further purified by a phenol-chloroform-isoamyl alcohol extraction and ethanol precipitation. rRNA was depleted with the Ribo-Zero rRNA Removal Kit (Epicentre) according to the manufacturer's instructions. The ribosome-protected RNA fragments (RPF) of between 26 and 34 bp were separated by electrophoresis for 65 min at 200 V using 15% polyacrylamide TBE-urea gel (Invitrogen), and eluted in 400 μL of RNA gel extraction buffer, which was composed of 300 mM sodium acetate pH 5.5, 1 mM EDTA, and 0.25% (w/v) SDS. The samples were frozen for 30 min at −80 °C and then incubated at 37 °C for 4 h. The eluted RNAs were isolated by ethanol precipitation and purified once again with the RNeasy MinElute Column (Qiagen, Hilden, Germany) using the manufacturer's protocol. The enriched RPFs were then denatured for 90 s at 80 °C and incubated for 1 h at 37 °C with 5 mL of 10× T4 Polynucleotide Kinase (PNK) buffer (NEB), 20 U of SUPERase-In RNase Inhibitor, and 10 U of T4 PNK (NEB) to dephosphorylate the 3′ end. The dephosphorylated RNAs were purified using the RNeasy MinElute Column (Qiagen). The sequencing library was constructed from the end-repaired RPFs using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB) according to the manufacturer's instructions. The final library of approximately 150-160 bp was size-selected by gel electrophoresis for 90 min at 100 V using a 2% agarose gel that was dyed with SYBR Gold Nucleic Acid Gel Stain (Bio-Rad, Hercules, CA, USA). The concentration of the final library was measured using a Qubit 2.0 Fluorometer (Invitrogen) and the Qubit dsDNA HS Kit. The size distribution was assessed using the Agilent 2200 TapeStation System (Agilent). The constructed library was sequenced on the Illumina HiSeq. 2500 platform using the 50-bp single-end read recipe (Fig. 1b).

Data processing of ribosome profiling reads. The libraries of seven samples of S. avermitilis-ex-
cept for the E1 sample-and six samples of S. venezuelae-except for the T2 and S2 samples-were prepared and sequenced twice to increase the output and merged for further data processing ( Table 2). The sequencing results were de-multiplexed and processed by CLC Genomics Workbench (CLC Bio). A total of 113,065,267 to 763,831,282 raw reads were generated for each replicate and were exported in the FASTQ format for the data upload. The reads were then mapped to the PhiX control sequences (NCBI Genbank accession number: NC_001422) to eliminate the PhiX control reads with the following parameters: mismatch cost: 2; insertion cost: 3; deletion cost: 3; length fraction: 0.9; similarity fraction: 0.9; and non-specific matches were randomly mapped. A total of 112,376,633 to 661,109,040 reads were unmapped. As these reads were sequenced from the 5′ end of the enriched RPF to 50 bp downstream, which is longer than the size-selected RPF (26 to 34 bp), the 5′ end sequences of the 3′ adapter sequence of the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB) were also included. To remove the adapter sequences from the reads prior to mapping, the sequences were trimmed by the following parameters: action: remove adapter; strand: minus; mismatch cost: 2; gap cost: 3; internal match minimum score: 3; and end match minimum score: 3. Ultimately, the removed adapter sequence was 5′−ATACGAGATNNNNNNCGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTT−3′, in which the NNNNNN sequences were CACTGT, ATTGGC, TACAAG, and TTTCA for index 5, 6, 12, and 19, respectively. The reads were additionally trimmed based on their overall quality (score: 0.05, maximum ambiguous nucleotides: 2) and length (>15 bp). The trimming steps yielded 90.68 to 98.47% of the PhiX control unmapped reads. To confirm the data quality and reproducibility of the reads to analyze the translational abundance of the genes, the reads were mapped to their genome sequence. A total of 84,947,464 to 590,644,871 reads with an average read length of 25.8 to 33.7 bp were mapped with random mapping of non-specific matches, while 1,833,155 to 103,819,037 reads with an average read length of 25.8 to 33.1 bp were mapped with ignored mapping of non-specific matches (mismatch cost: 2; insertion cost: 3; deletion cost: 3; length fraction: 0.9; similarity cost: 0.9). The overall statistics of the data processing are summarized in Table 2. The mapped information was exported in a BAM file format, and the number of mapped reads at each genomic position was counted as the read count. Normalized read value and principal component analysis (PCA) plots were generated using the DESeq. 2 package in R 34 .

Data Records
Raw read FASTQ files, trimmed read FASTQ files, mapped read BAM files, and the gene expression text files of all samples were uploaded to the public databases (Tables 1 and 2). Raw read FASTQ files of RNA-Seq and ribosome profiling of three species (S. avermitilis, S. clavuligerus, S. tsukubaensis) were deposited at the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) 35-37 . Raw read FASTQ files of RNA-Seq and ribosome profiling of S. lividans were deposited at the European Nucleotide Archive (ENA) 38 . Raw read FASTQ files of RNA-Seq of S. venezuelae were deposited at the ENA 39 . Raw read FASTQ files of ribosome profiling of S. venezuelae were deposited at the NCBI SRA 40-47 . Trimmed read FASTQ files and mapped read BAM files of the raw read FASTQ files in the NCBI SRA (RNA-Seq and ribosome profiling data of S. avermitilis, S. clavuligerus, S. tsukubaensis, and ribosome profiling data of S. venezuelae) were deposited at the ENA with a new accession 48 . Trimmed read FASTQ files and mapped read BAM files of the raw read FASTQ files in the ENA (RNA-Seq and ribosome profiling data of S. lividans, and RNA-Seq data of S. venezuelae) were also deposited at the ENA with the same accession as each corresponding raw read FASTQ file 38,39 . The gene expression profile as raw read counts of www.nature.com/scientificdata www.nature.com/scientificdata/ RNA-Seq and ribosome profiling data of S. avermitilis 49 , S. clavuligerus 50 , S. tsukubaensis 51 , and ribosome profiling data of S. venezuelae 52 are available in a text file format in the Gene Expression Omnibus (GEO) database. Also, the gene expression profiles of all datasets (RNA-Seq and ribosome profiling of the five species), including raw read counts, DESeq. 2 normalized values, fold change values between growth phases, and p-values for the fold changes, are available in a text file format in the Figshare 53 . The raw read FASTQ data of S. clavuligerus in NCBI SRA 36 was published in the previous study 21 . Also, the raw read FASTQ data of S. lividans in ENA 38 was published in the previous study 22 . Note that the ribosome profiling data of Streptomyces griseus was uploaded under the same accession with those of S. venezuelae, but they are not described in this study.  www.nature.com/scientificdata www.nature.com/scientificdata/ remained, which indicated high-sequencing quality. The remaining reads were used as input to generate sequencing QC reports in the CLC Genomics Workbench to validate the quality of the reads. At first, the overall read lengths were extremely long, corresponding to the sequencing read recipe (Fig. 2a, Table 1). For the four species that were sequenced with the 100-bp read recipe, the percentage of read lengths that were over 100 bp was more than 97.9%, and for S. tsukubaensis, which was sequenced with the 50-bp read recipe, the percentage of read lengths that were over 50 bp was more than 93.8%. Further, more than 98.6% (S. avermitilis), 98.9% (S. clavuligerus), 98.9% (S. lividans), 98.9% (S. venezuelae), and 96.4% (S. tsukubaensis) of the total reads exhibited an average Phred score of greater than 30, which indicates 99.9% base call accuracy (Fig. 2b). In addition, the quality of each base of the obtained reads was examined. The overall base positions of the sequencing reads were highly covered, and even the lowest average values of the coverage were 97.5% (S. avermitilis), 97.6% (S. clavuligerus), 97.8% (S. lividans), 98.3% (S. venezuelae), and 95.4% (S. tsukubaensis) at the last position, respectively (Fig. 2c). Moreover, the median values of the Phred scores per base position of the reads were consistently high across reads, with 40 scores in four species and 38 scores in S. tsukubaensis (Fig. 2d). From these quality validation results, we validated the quality of all obtained RNA sequencing reads for subsequent analysis.

RNA-Seq read quality validation.
Assessment of transcriptome data. The qualified reads were mapped to each reference genome with a uniquely mapped percentage that ranges from 76.91% to 95.09% (S. avermitilis), 52.19% to 90.88% (S. clavuligerus), 60.15% to 87.08% (S. lividans), 67.35% to 86.23% (S. venezuelae), and 76.33% to 96.26% (S. tsukubaensis) ( Table 1). The number of uniquely mapped reads at each gene was counted and normalized using the DESeq. 2 package in R 33 to reduce variation between samples. Using the normalized values, principal component analysis (PCA) was performed, which validated the high reproducibility of the sequencing data (Fig. 2e). The distribution of log 2 (DESeq normalized value + 1) broadly ranged from 0 to 20 in the different growth phase samples (Fig. 2f).
Ribosome profiling read quality validation. A total of 40 ribosome profiling reads were obtained from five Streptomyces species at four time points as duplicates. Unlike the RNA-Seq data, the trimmed reads were considered as raw sequences of the enriched RPF sequences, as additional PhiX control and adapter sequences that were involved in the ribosome profiling steps must be removed (Fig. 1a,b). Since the RPF fragments were selected by size, ranging from 26 to 34 bp, the 3′ end of the total 50 bp sequencing read contained non-RPF sequences, such as the 3′ adapter sequences, which do not represent the quality of the RPF reads. Thus, the QC reports on the trimmed reads were exported from CLC Genomics Workbench (Qiagen) to assess the quality of the RPF reads. The read length distribution exhibited a broad range from 20 to 40 bp, with one or two enriched peaks (Fig. 3a). The enriched peak sizes were comparable to the monosome-protected sizes, and they varied for different species, while they were more conserved for different growth phase samples of the same species. The differences in RNA degradation efficiency of RNase I or MNase across species may be the primary reason for the observed size differences 54 . Further, the read quality that was measured by the average Phred scores was generally high in all samples of the five species; the quality of more than 94% of the reads was higher than Q20, and more than 80% were higher than Q30 (Fig. 3b). Both per-sequence and per-base analyses of the read quality were observed. As the read size ranged mostly between 25 and 35 bp after adapter trimming, the base number coverage at each position of the 50 bp read dramatically decreased at the 3′ end (Fig. 3c). In terms of species, most of them exhibited the highest decline at 28 to 30 bp, which was consistent with the read length distribution (Fig. 3a). Given the base coverage, the median Phred score per base was demonstrated to be from 1 to 35 bp (Fig. 3d). The overall median of the quality score was approximately Q38, while the median score at the 5′ end was slightly lower than that of the middle section, and the score at the 3′ end of select species showed dramatic reductions. The low quality at the 3′ end may be due to some portions of identical long reads, which were somehow enriched during the size selection step of library construction, which stimulates wrong base calling. For the S. clavuligerus E2 sample, enriched peaks of less than 20 bp in length for approximately 10% of the total reads were unexpectedly found, along with decreased coverage at 15 bp, but these peaks did not seem to affect the overall read quality (Figs. 3a,c,d). Overall, most of the reads were shorter than 35 bp, and the read quality of all samples was high and suitable for downstream analyses.
Assessment of translatome data. To examine the additional quality of the reads for the translational abundance of each gene, the trimmed reads were mapped to their corresponding genome. Based on the mapping parameter, some reads would be non-specifically aligned to more than one genomic position due to highly repetitive genomic regions, including rRNA genes. Approximately, 43.1 to 590.6 million reads (75.3 to 96.8% of the trimmed reads) were mapped when the non-specifically matched reads were randomly assigned to one of the mapped positions, while 1.8 to 103.9 million reads (1.3 to 54.9% of the reads) were uniquely mapped when the non-specifically matched reads were excluded (Table 2). These results suggest that the non-specifically matched reads were generally more than half of the total mapped reads. Further, the proportion of these reads varied in different samples even within the same species, which is because the rRNA was enriched during the monosome recovery step, and the efficiency of rRNA removal differed across samples 55 E1  E2  T1  T2  L1  L2  S1  S2   E1  E2  T1  T2  L1  L2  S1  S2   E1  E2  T1  T2  L1  L2  S1  S2   E  T  L  S   E1  E2  T1  T2  L1 L2 S1 S2 Fig. 2 Read quality analysis of RNA-Seq samples of five Streptomyces species at four growth phases. The replicate of each growth phase is represented as "1" or "2" after the growth phase. www.nature.com/scientificdata www.nature.com/scientificdata/ obtained are, based on several bacterial transcriptome studies, considered sufficient for analysis of the whole translational profile and differential expression levels of genes, as 1 to 5 million reads are suggested for high statistical significance [56][57][58][59] . Among the uniquely mapped reads, some reads were mapped within RNA genes, rather   0  10  20  30  40  50  60  70   E1  E2  T1  T2  L1  L2  S1  S2   E1  E2  T1  T2  L1  L2  S1  S2   E1  E2  T1  T2  L1  L2  S1  S2   E  T  L  S   E1  E2  T1  T2  L1 L2 S1 S2 Fig. 3 Read quality analysis of ribosome profiling samples of five Streptomyces species at the four growth phases. The replicate of each growth phase is represented as "1" or "2" after the growth phase. www.nature.com/scientificdata www.nature.com/scientificdata/ than protein-coding genes, which mostly corresponded to tRNA genes. These reads may be the fragments of tRNA and rRNA that were bound to the ribosome and then enriched during monosome recovery 60 . Therefore, further validation was performed using only the mapped reads of the protein-coding genes. A total of 1.1 to 19.5 million reads were mapped to protein-coding genes that were 4.3 to 81.0% of the uniquely mapped genes, which indicates a high ratio of tRNA gene-mapped reads (Table 2). To validate the mapped read quality, the reproducibility of the mapped read number among biological replicates was investigated by PCA. All replicates were found to exhibit high reproducibility (Fig. 3e). The mapped read quality for quantitative analysis, such as the differential translational abundance of genes during growth, was examined by the distribution of the normalized values at four different growth phases, as described in the "Methods" section. The overall log 2 value (DESeq normalized value + 1) broadly ranged from 0 to 20, which was considered significant to analyze the translational abundance in different growth phases (Fig. 3f). In conclusion, the mapped reads were confirmed to exhibit high quality in terms of sequencing depth, reproducibility, and translational abundance.

Code availability
Versions and parameters of all the bioinformatic tools that were used in this work are described in the "Methods" section.