rpoB and efp are stable candidate reference genes for quantitative real-time PCR analysis in Saccharopolyspora spinosa

Abstract

Spinosad (spinosyn A and spinosyn D), the secondary metabolite produced by Saccharopolyspora spinosa, is a potent insecticide with low effects on the environment and mammals. Strategies such as metabolic engineering, mutagenesis and fermentation process optimization have been employed for its production enhancement. Quantitative real-time polymerase chain reaction(qRT-PCR) is one of the preferred methods for the evaluation of gene transcript levels, whose accuracy and sensitivity depend on the normalization using optimal reference genes. However, no single reference gene is universally appropriate for all strains under various conditions. In this study, the transcriptional stability of 35 candidate reference genes, including 23 traditionally used reference genes and 12 novel reference genes identified in S. spinosa homologous species, was analysed in S. spinosa ATCC 49460 and spinosad high-yield strain S. spinosa S3-3 under three fermentation phases. The transcriptional stability of these genes was assessed by three statistical algorithms, geNorm, NormFinder and BestKeeper. The overall rankings suggested that rpoB and efp were the most stable candidate reference genes, and they may be the most promising reference genes for further study on the measurement of expression levels of target genes involved in the biosynthetic process of spinosad.

Keywords:

• Saccharopolyspora spinosa
• spinosad
• quantitative real-time PCR
• reference genes
• fermentation periods

Introduction

Saccharopolyspora spinosa, a kind of soil actinomycete, has attracted great attention because of its secondary metabolites, spinosad. Spinosad (spinosyn A and spinosyn D) is a member of polyketide-derived macrolide insecticide with high effectiveness and broad-spectrum insecticidal activity [1]. Taking the advantages of low mammalian and avian toxicity and low environmental impact, spinosad was awarded ‘The Presidential Green Chemistry Challenge (U.S.A.)’ three times and was permitted to be used on organic food by the European Union, leading to a high market demand [2–4].

In order to increase the fermentation yield of spinosad, researchers had made great efforts through chemical mutagenesis or fermentation optimization of S. spinosa [3, 5]. Recently, the whole genome sequence of S. spinosa has been published, which enables further studies on the construction of high-yield spinosad strains by metabolic engineering [6]. The whole genome sequence of S. spinosa laid the foundation to search for new target genes in S. spinosa to increase spinosad yield through molecular biology approach. Thereafter, transcriptional level research played a key role in revealing the regulation network and catalytic actions during cell processes.

Several methods have been used to investigate gene transcription levels, such as Northern blotting, semi-quantitative reverse transcription polymerase chain reaction (RT-PCR), in situ hybridization and quantitative real-time PCR (qRT-PCR). qRT-PCR is the most widely used technique to investigate gene expression due to its high sensitivity, specificity, accuracy, reproducibility, convenience and high throughout [7–9]. However, multiple factors can affect the accuracy of the qRT-PCR data, including RNA quality and quantity, cDNA synthesis methods and its concentration, the efficiency of amplification and so on [10]. For normalization and precise quantification of the expression of the target genes, the selection and validation of stable references genes are crucial steps before proceeding to gene expression analysis.

A number of housekeeping genes such as 16S ribosomal RNA (16S rRNA), ribosome protein L13 (rbL13), RNA polymerase factor, DNA polymerase factor and elongation factor have been used as reference genes in different expression profiling studies in bacteria [11–14]. However, it is reported that the expression of housekeeping genes may vary among various species and different experimental conditions. To identify the key steps in spinosad biosynthesis of S. spinosa, the most commonly used reference gene was DNA polymerase sigma factor subunit A (sigA) [15–17]. Also, 16S rRNA was selected as the reference gene for revealing the metabolism difference between Saccharopolyspora pogona and S. spinosa, as well as seeking the link between primary metabolism and spinosad biosynthesis [18, 19]. And ribosome protein L13 (rbL13) was also selected to determine the transcript levels of rex, ctyA and ctyb, which were related to the redox control in S. spinosa [20]. In the study of heterologous expression of spinosyn biosynthetic gene cluster, DNA gyrase subunit A (gyrA) was chosen as the reference gene [21]. Additionally, the polyketide synthase (spnA), which was involved in spinosad biosynthesis, was also considered to be a reference gene. It was reported that 16S rRNA and rbL13 were considered to be the most suitable reference genes in different S. spinosa strains and under different culture conditions [22]. Overall, it is obvious that the normalization of reference genes in S. spinosa has not been standardized.

In this study, the most stable reference genes for gene transcriptional analysis in S. spinosa were screened with two different strains and three fermentation phases. The study will contribute to more complex studies, for instance, searching for the limiting enzyme during spinosad biosynthesis and new target genes for S. spinosa yield improvement, as well as exploring other macrolide secondary metabolites.

Materials and methods

Strains and cultivation

The wild-type strain, S. spinosa ATCC 49460 was purchased from China General Microbiological Culture Collection Center (CGMCC). The high-yield mutant, S. spinosa S3-3 derived from strain ATCC 49460, was obtained by nitrosoguanidine mutagenesis, UV irradiation combined with chloramphenicol screening successively, and shows stable inheritance.

For strain activation, 0.5 mL of spore suspensions (approx. 10⁷–10⁸ spores/mL) was inoculated to seed medium and incubated for 54 h at 29 °C, 240 rpm on a rotary shaker. Each litre of seed medium contained the following: 4 g yeast extract, 4 g tryptone, 4 g casein hydrolysate (acid), 10 g glucose, 1.36 g K₂HPO₄, 0.5 g MgSO₄, with a pH of 7.2 ± 0.1 before autoclaving. Next, 10% (v/v) of the seed culture was inoculated into 30 mL of seed medium to prepare fermented seed liquid. Then, 0.3 mL of seed culture (10%, v/v) was inoculated into 30 mL of fermentation medium in 250-mL Erlenmeyer flasks and cultured for 7 d at 29 °C, 240 rpm. The composition of the fermentation medium (per litre) was 10 g peptonized milk, 20 g cottonseed protein, 20 g dextrin, 1 g NaCl, 1.0 g MgSO₄, 0.5 g KH₂PO₄, 3 g CaCO₃, 10 g yeast extract, 60 g glucose, 20 g soybean oil. The pH value was adjusted to 7.2 with 1 mol/L of NaOH before autoclaving. All experiments had three biological replicates.

Analysis of spinosad production

The yield of spinosad was assessed using high-performance liquid chromatography (HPLC) as described by Huang et al. [15]. For determination of spinosad productivity, fermentation broth was extracted with 3× volume of methanol overnight at 4 °C, protected from light. Then, the leaching liquor was centrifuged at 14,000 rpm for 15 min before being analysed by HPLC. Samples (10 μL) were injected into an Agilent C18 reverse-phase column (5 μm, 150 × 4.6 mm) operated at 40 °C using an isocratic mobile phase composed of 45% (v/v) methanol, 45% (v/v) acetonitrile, and 10% 6.5 mmol/L ammonium acetate with a flow rate of 1.0 mL/min using a UV detector at 244 nm.

RNA isolation and cDNA synthesis

Total RNA was extracted using TRNzol reagent (Tiangen Biotech Co., LTD. cat. no. DP424) from the frozen powder of S. spinosa mycelia of different fermentation time [23]. The possibility of contamination of genomic DNA was eliminated with RNase-free DNaseI (Takara Bio, Dalian, China; cat. no. 2270 A), and then, the digested RNA was used as a template for the synthesis of cDNA with Reverse Transcriptase (Promega Co., cat. no. M170A) with random hexamer primer (Takara Bio, Dalian, China; cat. no. 3801) according to the manufacturer’s protocol (https://www.takarabiomed.com.cn/DownLoad/3801.pdf). During the treatment of RNA, RNase Inhibitor (Takara Bio, Dalian, China; cat. no. 2313 A) was added to prevent the degradation of RNA. The integrity and size distribution of the RNA were verified by 2% agarose gel electrophoresis, and the concentrations of RNA and cDNA were determined spectrophotometrically.

Candidate gene selection and primer designing

It was reported that the reliable reference genes of the model organism Streptomyces coelicolor were identified by transcript abundance analysis [24]. Based on the research, we compared the amino acid sequences of selected reference genes of S. coelicolor through BLAST in the NCBI database, and then, 12 candidate genes were found. Besides, according to previous studies, another 23 commonly used housekeeping genes such as 16S rRNA, rbL13, sigA, GAPDH, tufA, hrdB and so on were included in the study [12–14, 22, 25]. All primers for candidate genes were designed by software online (https://sg.idtdna.com). The definition and primer sequences of genes are listed in Table 1.

Table 1. Description of candidate reference genes, primer sequences and amplicon size in the study.

Gene	Gene description	Primer sequence (5’–3’)	Amplicon size (bp)	qRT-PCR efficiency (%)	R^2
rpoC	DNA-directed RNA polymerase subunit beta	F-CCGTGGTTGGAGATGAAGTAC	109	101.79	0.998
		R-GTATGCGTGGTCTGGTGTC
sigA	Major sigma factor	F-CGCAGACCACCAAGCCAGAC	125	98.19	1.000
		R-CTCTTCCTGTGCGGTGAGCA
GAPDH	Type I glyceraldehyde-3-phosphate dehydrogenase	F-GTAGGTATCAACGGCTTCGG	144	96.27	0.999
		R-GAGGATGCTGTCGTACTTGAG
16S rRNA	16S ribosomal RNA	F-ACCCTTGCCCTATGTTGCC	73	107.89	0.994
		R-CCTTCCTCCGAGTTGACCC
tufA	Elongation factor Tu	F-GAAGTTCGAGAGGGACAAGC	102	102.58	0.999
		R-TCGTGCAGAACCTTGGTG
Ts	Elongation factor Ts	F-GAGTTCATCGAGCTGGCC	146	97.02	0.993
		R-AGCTTCTCGCCGATCTTG
greAB	Elongation factor GreAB	F-GTCTACTCGCCCGATTCAC	114	109.54	0.992
		R-GGCACGCATCAAGGTCAC
greA	Transcription elongation factor GreA	F-GAAGATCAACGAGGCCCG	107	101.20	0.999
		R-CCTGAAGCTGGCGGATG
lepA	Elongation factor 4	F-AACCTGAACTTGATCTCCACC	114	92.24	0.997
		R-GTCCTTGATCTTGCCCTCG
efp	Elongation factor P	F-TGAACTTCCAGCACGTCAAG	118	96.72	0.999
		R-GTTTCCACCTTCACCCCG
fusA	Elongation factor G	F-GGTGTTCGACGGTAAGGAAG	134	98.84	0.999
		R-GAACCGTGAAGTAGAAGTCCG
hprT	Hypoxanthine phosphoribosyltransferase	F-ACATCGCTTCTGTGCTCATC	150	108.96	0.993
		R-CGAAGTCGGTCATGAACATG
rbL13	Ribosomal protein L13	F-GGCGTAGACCTTGAGCTTC	79	102.35	0.998
		R-GCTCGAAAAGGCGATCAAG
rpoB	DNA-directed RNA polymerase subunit beta	F-TTTCTCCGATCCACGCTTC	135	97.05	0.999
		R-GTCTGGCTCTTGATCTCGC
ubQ	Ubiquitin	F-TGTCCCAGGAAAAGGTTCAG	128	106.28	0.999
		R-GTCGATTTCGTCCAGGATGG
hrdB	DNA-directed RNA polymerase subunit beta	F-GTATGCGTGGTCTGGTGTC	109	103.91	0.990
		R-CCGTGGTTGGAGATGAAGTAC
thyA	FAD-dependent thymidylate synthase	F-GTGTCGTTCTACCTGTCCG	80	107.81	1.000
		R-AGCTGTGAGTAGGAGAGGTG
gyrB	DNA topoisomerase (ATP-hydrolyzing) subunit B	F-TCAACAGCTACGCCAAGG	143	103.85	1.000
		R-CAGCTTGGTCTTGGTCTGG
hrdA	DNA-directed RNA polymerase subunit alpha	F-GTCGTACCCTGCTCTCCT	136	99.95	1.000
		R-ACCAGCTCCTTGATGTTCAG
hrd-omega	DNA-directed RNA polymerase subunit omega	F-GCAGGTCAGCTCCAAGTAC	73	104.73	1.000
		R-CGCGTAGTAGTCGTTGATCTG
hrdB2	RNA polymerase sigma factor	F-AGGTGATCAACAAGCTC	101	101.56	0.998
		R-CGGGGTGATGTCCATTTC
hprT2	Hypoxanthine phosphoribosyltransferase	F-ATCCTCAAGGACCTCGACC	117	98.69	0.999
		R-AGGATTGCGCGAACTCAGG
ubQ2	ubiquitin-like protein Pup	F-GATGCCATGTTGGACGAAATC	75	104.88	0.999
		R-TTTCTGGATGTACGAGCGC
RS0111025	TIGR03619 family F420-dependent LLM class oxidoreductase	F-GGCAAGTTCTGGACCATCAG	72	134.56	0.983
		R-CAGGTAGATCGGCAGGTG
RS0122145	RNA polymerase sigma-70 factor	F-CCACACCATGATCCTCGATG	72	135.79	0.996
		R-AGCTTGTCGGGATTGATCAC
ruvA	Holliday junction branch migration protein RuvA	F-CTGTGGACGTCGCTGATC	188	106.95	0.993
		R-TTCGTCAGGTTGCCATCG
mshA	D-inositol-3-phosphate glycosyltransferase	F-TCAAGAACGCCAACCTCG	198	26034.28	0.961
		R-GGGTGAAGCATTCGAGGTC
RS0103485	Thioredoxin	F-AAGTTCAGCGAGATGTACCAC	188	111.26	0.987
		R-ATCTTGGCGTTGACCTGAG
RS0138655	HAMP domain-containing histidine kinase	F-ATCTGGACAAGCGGATGC	100	131.54	0.999
		R-TTCGTGATCTGGTGCTTGAG
phoX	PhoX family phosphatase	F-CGGATCACCCTGAACACC	182	112.90	0.996
		R-TTGCCGAAGTAGAAGTGGATG
RS0100110	Acid phosphatase	F-GGCCGTAGTGTTCGCGA	171	106.44	0.999
		R-TGCGTCCACACCGTCCAGT
RS0117990	ABC transporter ATP-binding protein	F-AGCGTGTTCAAGGAGATGAC	181	332.97	0.973
		R-CGAACTCCACCTGCTTCTG
RS0103685	SAM-dependent methyltransferase	F-GATGTAGTCGATCGTGTCCTG	136	116.86	0.974
		R-AGCACTGGAAACCGAACG
RS0140840	5-formyltetrahydrofolate cyclo-ligase	F-ACTGGAATGGGCGGAGTT	180	154.91	0.987
		R-TCGTAGTAACCGGCACCG
RS0101290	Hypothetical protein	F-GTACAGCTTCGATTTGGCG R-ATTTGGACTCGGCACGAC	196	96.82	0.988

For qRT-PCR, the primer specificity was determined by the PCR products and melting curve analysis. The amplification efficiency (E) for each primer and the R² value were evaluated based on the standard curve generated from a 10-fold serial dilution of pooled cDNA (10°, 10⁻¹, 10⁻² and 10⁻³). The amplification efficiency was calculated using the given formula [E=(10^(–1/slope)–1 × 100] [26, 27].

SYBR green qRT-PCR

For analysis of the transcriptional levels of target genes, 0.07 μg of cDNA was used for qRT-PCR on the CFX96 touch Real-Time PCR System (Bio-Rad). Each 15 μL reaction system consisted of 1.0 μL cDNA (70 ng/μL), 0.75 μL each primer (10 μmol/L), 7.5 μL 2 × SYBR® Select Master Mix (Applied Biosystems; cat. no. 44729) and 5.0 μL RNase-free ddH₂O. The PCR procedures were as follows: 2 min at 50 °C for activating uracil-DNA glycosylase, 2 min at 95 °C for pre-denaturation, followed by 40 cycles of 15 s at 95 °C for denaturation and 30 s at 60 °C for annealing and extension. To reduce the experimental error, each sample had four parallels. The CFX Manager Software V3.1 was used for visualizing and analysing the raw data.

Data analysis

To analyse the expression stability of candidate reference genes, the cycle threshold (Ct) values of the candidate genes were analysed with three statistical tools: geNorm, NormFinder and BestKeeper [28–30].

Results

Fermentation and sampling

The whole fermentation period of S. spinosa was 7 days. According to the fermentation process, the S. spinosa strain showed different physiological statuses, and the samples of two different strains were collected at three fermentation stages: initial stage of spinosad production (fermentation for 3 days), logarithmic period of spinosad production (fermentation for 5 days) and stationary phase of spinosad production (fermentation for 7 days). As shown in Figure 1, at these three time points, the spinosad production of wild-type strain and high-yield mutant was comparable, while the yield of spinosad production of S. spinosa S3-3 was 6–16 times higher than that of S. spinosa ATCC 49460.

Figure 1. Spinosad production by two S. spinosa strains during different fermentation phases.

Primer amplification specificity and efficiency evaluation

To assess the primer specificity of the 35 candidate genes in our study, the amplified products were checked by 2% agarose gel electrophoresis. Figure 2 shows that the amplicon sizes were consistent with the expected sizes. The amplification of a single PCR fragment with no dimer formation and a single sharp peak in the melting curves (data not shown) ensured the specificity of the primers for the selected reference genes. In addition, the primer efficiency of all selected genes was higher than 90% with R² values of the standard curve ranging from 0.961 to 1.000 (the detailed description is shown in Table 1). It was reported that the acceptable amplification efficiency varied from 90% to 110% [26, 31]. Based on the principle, some of the genes, including RS0111025 (134.56%), RS0122145 (135.79%), mshA (26034.28%), RS0138655 (131.54%), RS0117990 (332.97%), RS0140840 (154.91%), were excluded.

Figure 2. Amplification product of candidate reference genes. PCR products in 2% agarose gel stained with ethidium bromide (EB). (A) PCR analysis of the 17 genes (rpoC, sigA, GAPDH, 16S rRNA, tufA, Ts, greAB, greA, lepA, efp, fusA, hprT, rbL13, rpoB, ubQ, hrdB and thyA); (B) PCR analysis of the 6 genes (gyrB, hrdA, hrd-omega, hrdB2, hprT2 and ubQ2); (C) PCR analysis of the 12 genes (RS0111025, RS0122145, ruvA, mshA, RS0103485, RS0138655, phoX, RS0100110, RS0117990, RS0103685, RS0140840 and RS0101290). M: D15000 + 2000 DNA marker (Tiangen Biotech Co., LTD. cat. no. MD116).

Candidate reference gene transcription-level profiling

The expression rates of all the selected reference genes were measured based on Ct values from the qRT-PCR for different samples. The mean Ct values gave insight into the approximate gene expression levels. As shown in Figure 3A, the mean Ct values for the remaining 29 candidate reference genes ranged from 9.82 for 16S rRNA to 36.65 for ruvA for S. spinosa ATCC 49460. Similarly, for S. spinosa S3-3, the mean Ct values ranged from 9.21 for 16S rRNA to 38.05 for ruvA (Figure 3B), which revealed a relatively wide range of expression differences between genes. Among these selected reference genes, 16S rRNA showed a maximum expression variation, whereas RS0140840, ruvA and RS0117990 were the least abundantly expressed genes with higher Ct values (>35). According to the MIQE guideline that the Ct value of the reliable reference gene should be less than 31 [32], the candidate genes ubQ2, ruvA, RS0101290, RS0103685, phoX and Ts were not suitable to be selected as the internal reference genes.

Figure 3. Threshold cycle (Ct) values of 29 candidate reference genes. The box plots of S. spinosa ATCC 49460 (A) and S. spinosa S3-3 (B) for Ct values from the qRT-PCR analysis. Note: For each candidate gene, the line inside the box is the median and the little square inside the box is the mean. The top and bottom lines of the box are 75th and 15th quartiles, respectively. The top and bottom whiskers represent the 5th and 95th percentiles.

geNorm analysis

The software geNorm was used to calculate the gene expression stability value (M value) to rank the candidate genes. Based on the geNorm algorithm, the gene with the least M value was considered as the most stable reference gene, whereas with the highest M value was regarded as the least stable gene [30]. For S. spinosa ATCC 49460, the M value of the selected reference genes excluding hrdB2 was lower than 1.5. The M value of the genes ranged from 0.815 to 2.374. The gene hprT2 (M = 0.815) was the most stable reference gene followed by hrdA (M = 0.873), RS0103485 (M = 0.877), rpoB (M = 0.883) and fusA (M = 0.9), while hrdB2 was the least stable candidate reference gene with the highest M value (2.374) (Figure 4(A); Table 2).

Figure 4. Expression stability analysis of each candidate reference gene in S. spinosa ATCC 49460 (A) and S. spinosa S3-3 (B) using geNorm.

Note: The lower stability value (M value) indicates higher stability of the candidate reference genes. The direction of the arrow indicates the most and least stable reference genes, respectively.

Table 2. geNorm analysis for 23 candidate reference genes in S. spinosa wild-type and mutant strains.

	ATCC 49460		S3-3			ATCC 49460		S3-3			ATCC 49460		S3-3
Gene	M value	Rank	M value	Rank	Gene	M value	Rank	M value	Rank	Gene	M value	Rank	M value	Rank
rpoC	1.177	19	0.712	4	efp	0.956	8	0.698	2	gyrB	1.161	18	1.129	22
sigA	0.932	6	1.052	19	fusA	0.9	5	0.777	7	hrdA	0.873	2	0.768	6
GAPDH	1.138	17	0.834	12	hprT	1.022	12	0.842	13	hrd-omega	1.116	15	0.965	17
16S rRNA	0.976	9	0.785	9	rbL13	1.136	16	0.81	10	hrdB2	2.374	23	2.33	23
tufA	1.013	11	0.86	15	rpoB	0.883	4	0.668	1	hprT2	0.815	1	0.86	14
greAB	1.286	21	0.715	5	ubQ	0.941	7	0.82	11	RS0103485	0.877	3	0.935	16
greA	1.113	14	0.976	18	hrdB	1.331	22	0.701	3	RS0100110	1.096	16	1.059	22
lepA	1.269	20	1.078	21	thyA	0.995	10	0.784	8

As for S. spinosa S3-3, except for hrdB2, the other 22 candidate reference genes showed M value less than 1.5. The M value of genes varied from 0.668 for rpoB to 2.33 for hrdB2. The gene rpoB (M = 0.668) was the most stable reference gene followed by efp (M = 0.698). Similarly with S. spinosa ATCC 49460, hrdB2 was the least stable candidate reference gene with the highest M value (2.33) (Figure 4(B); Table 2).

Moreover, geNorm also determines an optimum number of reference genes based on the pairwise variation (V_n/V_n+1) between sequential normalization factor. A value of V_n/V_n+1 less than 0.15 indicates that the value of n being the optimal number of candidate reference genes is enough to be used as an internal reference gene [30]. In this study, the values of V₂/V₃ in two strains were lower than 0.15, suggesting that two reference genes were sufficient for normalization of genes (Figure 5).

Figure 5. geNorm analysis to determine optimal number of reference genes required for gene normalization. Pairwise variation (V) analysis for S. spinosa ATCC 49460 (A) and S. spinosa S3-3 (B).

Note: V1 to V23 represent the variation in candidate reference genes ranked based on their stability. V1 is the variation for the most stable, whereas V23 is the variation for the least stable gene.

NormFinder analysis

The stability of the selected reference genes was further estimated by the NormFinder based on inter-group and intra-group variations, and the genes with lower values indicated higher stability [28]. Based on the stability value (S value) of all genes, hrdB2 was the least stable reference gene with S value higher than 1.5 either in S. spinosa wild-type strain or in mutant. hprT2, hrdA, efp, fusA and rpoB showed higher stability with lower values. For S. spinosa ATCC 49460, hprT2 was identified as the most stable reference gene with an S value of 0.101 followed by hrdA (S = 0.257) (Figure 6(A); Table 3). While for S. spinosa S3-3, geNorm analysis and NormFinder algorithm revealed almost similar results for selecting candidate reference gene. rpoB was ranked as the most stable reference gene with an S value of 0.094. efp was the third stable candidate reference gene with an S value of 0.15 (Figure 6(B); Table 3).

Figure 6. Expression stability analysis of each candidate reference gene in S. spinosa ATCC 49460 (A) and S. spinosa S3-3 (B) using NormFinder. The lower stability value (S value) indicates higher stability of the candidate reference genes. The direction of the arrow indicates the most and least stable reference genes, respectively.

Table 3. NormFinder analysis for 23 candidate reference genes in S. spinosa wild-type and mutant strains.

	ATCC 49460		S3-3			ATCC 49460		S3-3			ATCC 49460		S3-3
Gene	S value	Rank	S value	Rank	Gene	S value	Rank	S value	Rank	Gene	S value	Rank	S value	Rank
rpoC	0.606	16	0.157	4	efp	0.386	7	0.15	3	gyrB	0.643	19	0.73	22
sigA	0.408	8	0.595	21	fusA	0.327	5	0.33	8	hrdA	0.257	2	0.328	7
GAPDH	0.598	15	0.379	12	hprT	0.465	11	0.335	9	hrd-omega	0.609	18	0.549	18
16S rRNA	0.381	6	0.3	6	rbL13	0.609	17	0.393	13	hrdB2	1.586	23	1.585	23
tufA	0.46	10	0.414	15	rpoB	0.226	3	0.094	1	hprT2	0.101	1	0.365	11
greAB	0.735	21	0.192	5	ubQ	0.41	9	0.408	14	RS0103485	0.308	4	0.456	16
greA	0.543	14	0.499	17	hrdB	0.752	22	0.13	2	RS0100110	0.503	13	0.562	19
lepA	0.672	20	0.575	20	thyA	0.465	12	0.365	10

BestKeeper analysis

Another program, the BestKeeper could identify potential reference genes based on two significant criteria: standard deviation (SD) of the average Ct values and coefficient of correlation (r) value [29]. Higher r value means more stable candidate reference genes, and lower SD value represents higher stability. Here, the r value was considered as the prime element for evaluation. In the wild-type strain, greAB (r = 0.91) was ranked as the most stable reference gene, followed by tufA and rpoB with r values of 0.89 and 0.88, respectively. And hrdB2 was the least stable candidate reference gene with the minimum r value (Table 4), which was identical with the results from the geNorm and NormFinder analyses. For the mutant strain, rbL13 (r = 0.95) was the most stable gene, followed by rpoB (r = 0.94), thyA (r = 0.94) and greAB (r = 0.94). RS0103485 (r = 0.35), sigA (r = 0.48) and hrdB2 (r = 0.56) were less stably expressed genes with lower r values. Considering that an SD value less than 1.0 suggests consistence and stability of gene expression levels, the gene greAB was not eligible for a reference gene with SD value higher than 1.0.

Table 4. The 23 candidate reference genes stability and ranking analysed by BestKeeper.

	ATCC 49460			S3-3				ATCC 49460			S3-3				ATCC 49460			S3-3
Gene	r value	SD	Rank	r value	SD	Rank	Gene	r value	SD	Rank	r value	SD	Rank	Gene	r value	SD	Rank	r value	SD	Rank
rpoC	0.52	0.87	20	0.92	0.64	6	efp	0.80	0.59	7	0.90	0.55	7	gyrB	0.60	0.99	18	0.56	0.87	20
sigA	0.72	0.59	12	0.48	0.71	22	fusA	0.85	0.82	6	0.90	0.67	8	hrdA	0.67	0.75	14	0.80	0.67	14
GAPDH	0.86	1.05	5	0.87	0.70	10	hprT	0.63	0.94	17	0.83	0.73	13	hrd-omega	0.79	0.80	8	0.74	0.58	16
16S rRNA	0.64	0.27	16	0.86	0.42	11	rbL13	0.75	0.71	11	0.95	0.72	1	hrdB2	0.15	1.68	23	0.56	2.28	21
tufA	0.89	0.96	2	0.84	0.75	12	rpoB	0.88	0.57	3	0.94	0.53	2	hprT2	0.66	0.49	15	0.70	0.49	18
greAB	0.91	1.16	1	0.94	0.66	3	ubQ	0.76	0.48	9	0.93	0.74	5	RS0103485	0.69	0.66	13	0.35	0.51	23
greA	0.76	0.79	10	0.77	0.82	15	hrdB	0.44	1.04	21	0.88	0.59	9	RS0100110	0.36	0.84	22	0.68	1.03	19
lepA	0.59	0.98	19	0.72	0.99	17	thyA	0.88	0.78	4	0.94	0.66	3

Comprehensive analysis of expression stability

For a comprehensive assessment of the most stable reference genes, the above results from geNorm, NormFinder and BestKeeper were taken into consideration. As shown in Table 5, rpoB was the most stable reference gene followed by fusA and hprT2 in S. spinosa ATCC 49460. For S. spinosa S3-3, rpoB was the most stable reference gene followed by efp, greAB and rpoC. While fusA ranked the seventh in S. spinosa S3-3, comparing to efp ranked the sixth, we suggested efp as a preferable one. Overall, the gene rpoB ranked as the most reliable reference gene, and efp was the second most stable candidate reference gene in both S. spinosa wild-type strain and its mutant S3-3. hrdB2 and gyrB performed as the least stable candidate reference genes.

Table 5. Overall ranking of the candidate reference genes by the results from geNorm, NormFinder and BestKeeper.

Method	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23
	Ranking order of selected reference genes in S. spinosa ATCC 49460 at different fermentation stages (Better-Good-Average)
geNorm	hprT2	RS0103485	hrdA	rpoB	efp	sigA	fusA/ubQ		thyA	16S rRNA	hrd-omega	tufA	hprT	rbL13	gyrB	RS0100110	rpoC	GAPDH	greA	lepA	greAB	hrdB	hrdB2
NormFinder	hprT2	hrdA	RS0103485	efp	thyA	rpoB	ubQ	sigA	16S rRNA	fusA	hrd-omega	rbL13	gyrB	tufA	hprT	lepA	RS0100110	rpoC	GAPDH	greA	greAB	hrdB	hrdB2
BestKeeper	greAB	tufA	rpoB	thyA	GAPDH	fusA	efp	hrd-omega	ubQ	greA	rbL13	sigA	RS0103485	hrdA	hprT2	16S rRNA	hprT	gyrB	lepA	rpoC	hrdB	RS0100110	hrdB2
Recommended comprehensive ranking	rpoB	fusA	hprT2	hrdA	RS0103485	efp	tufA	ubQ	thyA/sigA		16S rRNA	GAPDH	greA	hprT	hrd-omega	greAB	rbL13	RS0100110	rpoC/gyrB		lepA	hrdB	hrdB2
	Ranking order of selected reference genes in S. spinosa S3-3 at different fermentation stages (Better-Good-Average)
geNorm	rpoB	efp	hrdB	rpoC	greAB	thyA	hrdA	rbL13	16S rRNA	fusA	ubQ	hprT	GAPDH	hprT2	tufA	hrd-omega	RS0103485	sigA	greA	gyrB	lepA	RS0100110	hrdB2
NormFinder	rpoB	efp	rpoC	hrdB	greAB	thyA	16S rRNA	hrdA	rbL13	ubQ	fusA	hprT	hprT2	GAPDH	hrd-omega	tufA	RS0103485	sigA	lepA	greA	gyrB	RS0100110	hrdB2
BestKeeper	rbL13	rpoB	greAB/thyA		ubQ	rpoC	efp	fusA	hrdB	GAPDH	16S rRNA	tufA	hprT	hrdA	greA	hrd-omega	lepA	hprT2	RS0100110	gyrB	hrdB2	sigA	RS0103485
Recommended comprehensive ranking	rpoB	efp	greAB	rpoC/hrdB		thyA	fusA	rbL13	16S rRNA	hrdA	ubQ	GAPDH	hprT	tufA	hprT2	greA	hrd-omega	RS0103485	lepA/RS0100110		sigA	gyrB	hrdB2

Discussion

qRT-PCR is an essential technique for analysing transcriptomes and investigating the regulatory molecular mechanisms [33]. But an inappropriate reference gene could result in invalid results. In general, a reference gene is not universal under different conditions. Therefore, it is extremely necessary to select suitable references to guarantee the accuracy of the transcriptional analysis of the target genes. In previous studies, stable reference genes have been identified in multiple species of cyanobacteria [34], bacteria [22, 35, 36], yeast [37], mycetes [38] and plants [39, 40] under different conditions.

This study is a systematic assessment of candidate reference genes in S. spinosa wild-type strain and high-yield mutant during different fermentation periods. Twenty-three traditionally used reference genes (encoding ribosomal RNA, ribosomal protein, elongation factor, RNA polymerase subunit, ubiquitin, RNA polymerase sigma factor, glyceraldehyde-3-phosphate dehydrogenase, etc.) and 12 novel genes from its homologous species Streptomyces coelicolor were selected as the candidate genes. The specificity of synthesized primers was determined by agarose gel electrophoresis and melting curve analysis, and amplification efficiency of primers was calculated using a series of 10-fold dilutions of pooled cDNA.

Three different statistical algorithms named geNorm, NormFinder and BestKeeper were regularly used to analyse the datasets to identify reference genes [28–30]. In this study, the geNorm and NormFinder programs displayed the same results of the most stable and least suitable reference genes, but differed from the order analysed by BestKeeper in two different strains. In the S. spinosa ATCC 49460, hprT2 and hrdB2 were identified as the reference gene with the highest and the lowest stabilities by geNorm and NormFinder softwares, respectively, while greAB and hrdB2 were the most and the least stable genes obtained by BestKeeper. As for S. spinosa S3-3, rpoB and hrdB2 were identified as the most and the least stable genes by geNorm and NormFinder, whereas in BestKeeper, rbL13 was ranked as the highest stable gene and RS0103485 as the lowest stable one. Similar results that the ranking orders obtained from geNorm and NormFinder varied from those obtained by BestKeeper have been described in previous reports [33, 41]. Based on the comprehensive evaluation of the results from three methods of analysis, the genes rpoB and efp were the most reliable reference genes in S. spinosa strains at different fermentation stages. And hrdB2 and gyrB exhibited relatively lower stability in S. spinosa strains.

rpoB, encoding DNA-directed RNA polymerase subunit beta, was identified as stable reference in Bifidobacterium animalis under different culture conditions [42]. Also, rpoB had demonstrated its transcriptional stability in Actinobacillus suis, Staphylococcus aureus, Xanthomonas, Listeria monocytogenes and other bacteria [43–46]. In addition, elongation factor evaluated among multiple species showed a quite stable expression level under different conditions [11, 22, 27]. Although hrdB was widely used as a reference gene for detection of transcription levels of target genes in Streptomyces [47], it was more stable in S. spinosa S3-3 than in S. spinosa ATCC 49460. Zhang et al. [22] found that 16S rRNA and rbL13 were the most suitable reference genes in different S. spinosa strains and under different culture conditions. On the contrary, our study showed rpoB and efp were ranked as top candidate reference genes in different strains, under the condition of different fermentation periods, yet 16S rRNA and rbL13 were less stable.

Conclusions

In summary, we have successfully evaluated and validated stable reference genes in different S. spinosa strains at different periods of fermentation. The analysis revealed the reliable reference genes in S. spinosa wild-type strain and high-producing spinosad mutant. The results showed that rpoB and efp were the best candidate reference genes for the normalization of gene expression in the process of spinosad synthesis. These two genes may be the most promising reference genes for future studies on gene function, spinosad biosynthesis process, as well as molecular mechanisms in S. spinosa.

Data availability

The data that support the findings reported in this study are available from the corresponding author upon reasonable request.

Disclosure statement

The authors declare that they have no conflict of interest. Xiaomeng Liu, Yunpeng Zhang, Tie Yin, Qi Li, Qiulong Zou and Xiaolin Zhang are affiliated with COFCO, and conceived the experiments, designed the experiments, performed the experiments, analysed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft; Kexue Huang is affiliated with Qilu Pharmaceutical, and authored or reviewed drafts of the paper.

Additional information

Funding

This work was supported by the National Nature Science Foundation of China [Grant No. 31670055].