Introduction

With the increasing popularity of gene expression analysis in biological research, quantitative real-time polymerase chain reaction (qRT-PCR) has become a critical and powerful tool for rapid and reliable quantification of mRNA transcriptional expression levels of target genes due to its high-throughput screening, sensitivity, simplicity, specificity and accuracy1,2. Relative quantification of target gene expression under certain stresses has been widely studied since the beginning of this century3. An accurate assay of gene expression through qRT-PCR relies on every step of sample preparation and processing, e.g., the integrity of purified RNA, the efficiency of reverse transcription, and the overall transcriptional activity of the tissues or cells analysed4; each step needs to be accurately normalized by stably expressed reference genes (RGs)5,6. Therefore, the selection of reliable RGs for normalization under given experimental conditions is a requirement for developing an accurate qPCR assay.

Housekeeping genes, such as the glyceraldehyde 3-phosphate (GAPDH), the actin gene (ACTIN), translation elongation factor EF-1 alpha (EF-1α), 18 s rRNA, 25 S rRNA and poly-ubiquitin (UBQ), have been commonly used as the normalization scalar in studies of relative quantification of plant target genes, some of which (EF-1α, GAPDH, ACTIN) have been identified as reliable RGs in certain plants under given experimental conditions7,8,9,10. However, to date, no RG has been found to exhibit perfectly stable expression in all plant species, even in the same tissue from the same plant species, but under different experimental conditions11,12,13. For instance, DcACTIN and DcUBQ have been identified as the top two stable RGs in carrot (Daucus carota L.) under abiotic stresses, but eIF-4α and GAPDH have been ranked in the top two RGs in carrots under hormone stimuli7; in tea plants (Camellia sinensis (L.) O. Kuntze), CsTIP41 was identified as the most stable RG for leaf development, but CsTBP was identified as the most stable RG for tea leaves under hormone stimuli14. Therefore, to avoid missing or overemphasizing potential biological changes of target gene expression, it is essential to identify optimum stable RGs for the proposed research object, for different tissues of the same species, for the same tissue of the same species under different biotic or abiotic stresses and their processing time.

Tea is one of the most important leaf-type woody cash crops in China, and the tender buds and leaves of this plant are the raw material for commercial tea. Since the publication of the draft genome sequence of C. sinensis var. sinensis15, the molecular mechanisms of aroma components biosynthesis, cold spells or resistance, drought resistance, barren tolerance, and other interactions of tea plants with environmental factors or with other organisms around them have been elucidated16,17,18,19,20. During the development of tea plant, it usually suffers serious damage from the infestation of insect herbivores all year round. Therefore, the chemical and molecular mechanisms under interactions between tea plants and their herbivorous pests need to be widely excavated to offer theoretical foundations for utilizing chemical signals between them to control tea pests or breeding new insect-resistant tea varieties. The RGs used previously in the studies of herbivores (Ectropis obliqua, Empoasca onukii) induced tea plant defensive responses at the gene transcriptional level, such as CsGAPDH and 18SrRNA21,22,23, were roughly selected from previously reported RGs without critical identification under given experimental conditions, which may lead to the deviation of the results to some extent and may also lead to the neglect of some important experimental phenomena. Therefore, it is important to define the RG for qRT-PCR analysis in tea plants under infestations of different pests and their related biotic stresses.

According to previous reports, CsACTIN1, Clathrin adaptor complex subunit (CsCLATHRIN1), CsEF1, CsGAPDH1, SAND family protein gene (CsSAND1), Tap42-interacting protein of 41 kDa (CsTIP41), Ubiquitin-conjugating enzyme (CsUBC1), Polypyrimidine tract-binding protein (CsPTB1), alpha-1 tubulin (CsTUA1) and TATA-box binding protein gene (CsTBP1) are frequently used as stable RGs in the process of mRNA expression analysis (Tables 1 and 2)20,24,25,26,27,28,29. In the present study, we measured mRNA transcriptional levels of the above mentioned ten RGs in different tissues of tea plants in circadian rhythms, jasmonic acid-treated tea leaves, T. aurantii infested tea leaves, E. onukii infested tea leaves, and tea leaves treated with mechanical damage plus E. obliqua regurgitant. The results were evaluated by BestKeeper, geNorm, NormFinder and the ΔCt method to identify the most stably expressed RGs firstly; secondly, RefFinder was used to integrate the results to determine the most stable RG for each treatment. Finally, to demonstrate the importance of stable RGs in the normalization process of tea plants under infestations of different pests or their related biotic stresses, CsMYC2, CsOPR3, CsPAL and CsPALc were chosen as the target genes for validation. As we all know, MYC2 was a key transcription factor of JA signaling pathway30; OPR3 is the isoenzyme relevant for JA biosynthesis22 and PAL were closely associated with the accumulation of endogenous SA31. The aim of this study was to select the most appropriate RGs for the gene expression analysis of tea plants under different biotic stresses.

Table 1 Ten housekeeping genes frequently used for qRT-PCR of tea plant.
Table 2 Sequence Information of the Candidate Reference Genes and Target Genes.

Results

Expression profiles of candidate reference genes

The expression level of RGs in all treatments is performed in terms of the cycle threshold number (Ct value). As shown in Fig. 1, the raw Ct values of all candidate RGs ranged from 13.90 (EF1) to 28.29 (TBP). EF1 (18.44), ACTIN1 (18.91), GAPDH1 (18.97) and TUA1 (19.23) were the most abundant transcripts, reaching the threshold fluorescence peak after 18 cycles. PTB1 (23.65), CLATHRIN1 (23.71), SAND1 (24.04) and TBP (24.08) were expressed at the lowest levels. The raw Ct values of the four target genes ranged from 18.72 (PALc) to 27.26 (MYC2). More details were shown in Fig. S8.

Figure 1
figure 1

Expression Profiles of Ten Candidate Reference Genes and Four Target Genes in C. sinensis. The expression level of RGs in all samples is performed in terms of the cycle threshold number (Ct value). The data are expressed as box-whisker plots; the short bar in the box refers to the Ct mean value; the box represents the 25th–75th percentiles; the median is indicated by a bar across the box; the whiskers on each box represent the distribution of the Ct values; and the dark spots refer to extreme outliers.

Diurnal expression in different tissues

Leaf

The gene expression stability of ten candidate RGs for leaves with circadian rhythm was analyzed by using geNorm, NormFinder, BestKeeper and the ΔCt method. The results showed that the gene stability ranking as analyzed by BestKeeper differed from the ranking as analyzed by the other three methods. For example, geNorm, NormFinder and the ΔCt method identified UBC1 and CLATHRIN1 as the most stable 2 of the 10 RGs in all test periods (from 0:00 am to 22:00 pm), whereas BestKeeper identified GAPDH1 and CLATHRIN1 as the most stable 2 of the 10 RGs for diurnal expression in leaves. However, all four methods identified PTB1 as the most variable RG. According to the results from RefFinder, the stability ranking of RGs from the most to the least was as follows: UBC1 > CLATHRIN1 > GAPDH1 > TBP > EF1 > SAND1 > TUA1 > ACTIN1 > TIP41 > PTB1 (Table 3). With GeNorm (Fig. 2), all pairwise variation (Vn/n + 1) was below 0.15 (the recommended cut-off), indicating that the inclusion of an additional RG was unnecessary. Based on the ranking of the RGs by RefFinder, CLATHRIN1 and UBC1 were identified as the best combination for normalizing the diurnal expression in leaves (Tables 4 and 5).

Table 3 Ranking of 10 Reference Genes Expression under Different Experimental Manipulations.
Figure 2
figure 2

Optimal Number of Reference Genes for the Normalization of C. sinensis under Different Experimental Manipulations. The pairwise variation (Vn/n + 1) was analysed by geNorm software to determine the optimal number of RGs included in the qPCR analysis. Values less than 0.15 indicate that another RG will not significantly improve normalization.

Table 4 Ranking of 10 Reference Genes Expression in Different Processing Time under Different Experimental Manipulations.
Table 5 Summary of treatments and results.

Stem

GeNorm identified SAND1 and TIP41 as the most stable RGs in all test periods (from 0:00 am to 22:00 pm) (Table 4). NormFinder and the ΔCt method identified TUA1 and CLATHRIN1 as the most stable RGs. BestKeeper identified TUA1, CLATHRIN1 and SAND1 as the top three RGs. However, all four methods identified GAPDH1 as the most unstable RG (Table 3). According to the results from RefFinder, the stability ranking of RGs from the most to the least was as follows: TUA1 > SAND1 > CLATHRIN1 > UBC1 > TIP41 > PTB1 > ACTIN1 > TBP > EF1 > GAPDH1. Based on the ranking of the RGs by RefFinder, TUA1 and SAND1 were identified as the best combination for normalizing the diurnal expression in the stem (Table 5).

Root

NormFinder and the ΔCt method identified UBC1 and SAND1 as the most stable RGs, and ACTIN1 as the least stable RG in all test period (from 0:00 am to 22:00 pm) (Table 3). GeNorm identified SAND1 as the most stable RG. BestKeeper identified TIP41 as the most stable RG. According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: SAND1 > UBC1 > TBP > TIP41 > CLATHRIN1 > PTB1 > GAPDH1 > EF1 > TUA1 > ACTIN1. The results of the geNorm analysis revealed that all V values were below 0.15 (Fig. 2). Thus, SAND1 and UBC1 were identified as the best combination for normalizing the gene diurnal expression in roots (Table 5).

JA treatment

GeNorm, NormFinder and the ΔCt method identified CLATHRIN1, GAPDH1 and UBC1 as the top three stable RGs in all test periods (from 0.5 h to 48 h) (Table 3). BestKeeper identified SAND1, PTB1 and TIP41 as the top three stable RGs. All four methods identified TUA1 as the most unstable RG (Table 3). According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: CLATHRIN1 > GAPDH1 > UBC1 > TIP41 > PTB1 > SAND1 > TBP > ACTIN1 > EF1 > TUA1. The results of the geNorm analysis revealed that all V values were below 0.15 (Fig. 2). Thus, CLATHRIN1 and GAPDH1 were identified as the best combination for normalizing JA-treated leaves. With further analysis, RefFinder identified CLATHRIN1 and UBC1 as the best combination for JA treatment in the time interval from 0.5 h to 1.5 h, GAPDH1 and TIP41 as the best combination in the time interval from 3 h to 6 h, and CLATHRIN1 and GAPDH1 as the best combination in the time interval from 12 h to 48 h (Tables 4 and 5).

T. aurantii infestation

NormFinder and ΔCt identified ACTIN1 and UBC as the most stable 2 of the 10 RGs in all test periods (from 6 h to 48 h) (Table 4). BestKeeper ranked ACTIN1 and EF1 as the top two stable RGs. GeNorm ranked ACTIN1 and TBP as the top two RGs. According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: ACTIN1 > UBC1 > GAPDH1 > CLATHRIN1 > TBP > EF1 > PTB1 > SAND1 > TIP41 > TUA1 (Table 3). The results of the geNorm analysis revealed that almost all V values were below 0.15 (Fig. 2). Thus, ACTIN1 and UBC1 were identified as the best combination for normalizing T. aurantii-infested leaves. With further analysis, RefFinder identified ACTIN1 and UBC1 as the best combination in the time interval from 6 h to 24 h, ACTIN1 and EF1 as the best combination at 48 h (Tables 4 and 5).

E. onukii infestation

The GeNorm, NormFinder and ΔCt methods identified GAPDH1 and UBC1 as the most stable 2 of the 10 RGs, while PTB1 was the least stable RG in all test periods (from 12 h to 144 h) (Table 3). BestKeeper identified EF1, GAPDH1 and CLATHRIN1 as the top three stable RGs. According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: UBC1 > GAPDH1 > EF1 > TIP41 > SAND1 > CLATHRIN1 > TBP > TUA1 > ACTIN > PTB1. The results of the geNorm analysis revealed that all V values were below 0.15 (Fig. 2). Thus, UBC1 and GAPDH1 were identified as the best combination for normalizing E. onukii-infested leaves. With further analysis, RefFinder identified GAPDH1 and UBC1 as the best combination in the time interval from 12 h to 72 h, PTB1 and TBP as the best combination at 96 h, TIP41 and EF1 as the best combination in the time interval from 120 h to 144 h (Tables 4 and 5).

Mechanical damage and E. obliqua regurgitant treatment

GeNorm, NormFinder and the ΔCt method identified SAND1 and TBP1 as the most stable 2 of the 10 RGs, while TUA1 was the least stable RG in all test periods (from 1.5 h to 48 h) (Table 3). BestKeeper identified ACTIN1, CLATHRIN1 and TBP as the top three stable RGs. According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: SAND1 > TBP > CLATHRIN1 > PTB1 > ACTIN1 > TIP41 > UBC1 > EF1 > GAPDH1 > TUA1. The results of geNorm revealed that all V values were below 0.15 (Fig. 2). Thus, SAND1 and TBP1 were identified as the best combination for normalizing regurgitant-treated leaves. With further analysis, RefFinder identified TIP41 and TBP as the best combination in the time interval from 1.5 h to 3 h, TBP and CLATHRIN1 as the best combination at 6 h, and SAND1 and TBP as the best combination in the time interval from 12 h to 48 h (Tables 4 and 5).

Validation of proposed RGs

CsMYC2 was chosen as the target gene to validate the rationality of the recommended RGs used in diurnal expression analysis (Fig. 3A–C). The expression level of CsMYC2 in leaves at 14:00 pm was significantly higher than that in the time period from 0:00 am to 12:00 am (NF 9–10, F = 14.098, P = 0.000; P = 0.000; P = 0.000; P = 0.000; P = 0.000; P = 0.000) and that at 16:00 pm, 20:00 pm and 22:00 pm (NF 9–10, F = 14.098, P = 0.000; P = 0.000; P = 0.000) when normalized with the two unstable RGs, TIP41 and PTB1 (NF 9–10); these expression level trends were quite similar to that normalized with the combination of UBC1 and CLATHRIN1 (NF 1–2, F = 10.169, P = 0.000; P = 0.000; P = 0.003; P = 0.003; P = 0.005; P = 0.000), except for 10:00 am (NF 1–2, F = 10.169, P = 0.138) (Fig. 3A); the expression level of CsMYC2 in leaves at 4:00 am was significantly higher than that at 0:00 am and 2:00 am when normalized with the combination of UBC1 and CLATHRIN1 (NF 1–2, F = 10.169, P = 0.000; P = 0.002), but no significant differences were detected when normalized with the combination of TIP41 and PTB1 (NF 9–10, F = 14.098, P = 0.141; P = 0.485) (Fig. 3A). The expression level of CsMYC2 in stem at 10:00 am was significantly higher than that at the time period from 0:00 am to 6:00 am and from 12:00 am to 22:00 pm when normalized either with the combination of TUA1 and SAND1 (NF 1–2, F = 3.743, P = 0.000; P = 0.003; P = 0.019; P = 0.000; P = 0.003; P = 0.008; P = 0.002; P = 0.030; P = 0.001) or with the combination of EF1 and GAPDH1 (NF 9–10, F = 6.969, P = 0.000; P = 0.001; P = 0.005; P = 0.000; P = 0.000; P = 0.005; P = 0.000; P = 0.005; P = 0.006), except for 16:00 pm (NF 1–2, F = 3.734, P = 0.383; NF 9–10, F = 6.969, P = 0.000); however, the expression level of CsMYC2 in stem at 16:00 pm was significantly higher than that at 12:00 am and 18:00 pm when normalized with the combination of TUA1 and SAND1 (NF 1–2, F = 3.734, P = 0.030; P = 0.023), and no significant differences were detected when normalized with the combination of EF1 and GAPDH1 (NF 9–10, F = 6.969, P = 0.145; P = 0.256) (Fig. 3B). The expression level of CsMYC2 at 16:00 pm in root was significantly higher than that at4:00 am, 12:00 am, 14:00 pm, 20:00 pm and 22:00 pm when normalized with the most stable combination of SAND1 and UBC1 (NF 1–2, F = 3.610, P = 0.013; P = 0.000; P = 0.000; P = 0.002; P = 0.003;), but the expression level of CsMYC2 at 16:00 pm has no significant differences with that at all the time points (NF 9–10, F = 3.972, P = 0.521; P = 0.080; P = 0.464; P = 0.179; P = 0.604; P = 0.173; P = 0.360; P = 0.789; P = 0.525; P = 0.200), except for 10:00 am(NF 9–10, F = 3.972, P = 0.001), when normalized with the most unstable combination of TUA1 and ACTIN1 (NF 9–10) (Fig. 3C).

Figure 3
figure 3

Validation of the gene stability measure. Expression profiles of CsMYC2, CsOPR3, CsPAL and CsPALc under different experimental conditions using different RGs. (A) Diurnal expression profile of CsMYC2 in leaves, NF (1–2) were UBC1 and CLATHRIN1, NF (9–10) were TIP41 and PTB1; (B) Diurnal expression profile of CsMYC2 in stems, NF (1–2) were TUA1 and SAND1, NF (9–10) were EF1 and GAPDH1; (C) Diurnal expression profile of CsMYC2 in roots, NF(1–2) were SAND1 and UBC1, NF (9–10) were TUA1 and ACTIN1; (D) Expression profile of CsOPR3 at 3 h normalized with the best combination (GAPDH1 and TIP41) at 3 h, the best combination (CLATHRIN1 and UBC1) at 0.5–1.5 h, and the best combination (CLATHRIN1 and GAPDH1) at 12–48 h RGs under JA treatment; (E) Expression profile of CsPAL at 48 h normalized with the best combination (ACTIN1 and EF1) at 48 h, and the best combination (ACTIN1 and UBC1) at 6–24 h RGs under T. aurantii infestation; (F) Expression profile of CsPALc at 96 h normalized with the best combination (PTB1 and TBP) at 96 h, the best combination (GAPDH1 and UBC1) at 12–72 h, and the best combination (TIP41 and EF1) at 120–144 h under E. onukii infestation; (G) Expression profile of CsOPR3 at 6 h normalized with the best combination (TBP1 and CLATHRIN1) at 6 h, the best combination (TIP41 and TBP) at 1.5–3 h, and the best combination (SAND1 and TBP) at 12–48 h RGs under E. obliqua infestation; (H) Expression profile of CsOPR3 normalized with the stable and unstable RGs at 3 h under JA treatment. NF1 was GAPDH1, NF (1–2) were GAPDH1 and TIP41, NF10 was ACTIN1, NF (9–10) were TUA1 and ACTIN1; (I) Expression profiles of CsPAL normalized with the stable and unstable RGs at 6 h under T. aurantii infestation. NF1 was ACTIN1, NF (1–2) were ACTIN1 and UBC1, NF10 was TUA1, NF (9–10) were PTB1 and TUA1; (J) Expression profile of CsPALc normalized with the stable and unstable RGs at 96 h under E. onukii infestation. NF1 was PTB1, NF (1–2) were PTB1 and TBP, NF10 was TUA1, NF (9–10) were TIP41 and TUA1; (K) Expression profile of CsOPR3 normalized with the stable and unstable RGs at 6 h under E. obliqua infestation. NF1 was TBP, NF (1–2) were TBP and CLATHRIN1, NF10 was TUA1, NF (9–10) were EF1 and TUA1; Data are means ± SE. One-way ANOVA (Tukey’s test) was used to analyze significant difference among treatments (A~C,F,G,J,K); different letters indicate significant differences among treatments (lowercase letters, P < 0.05; uppercase letters, P < 0.01). Two samples were compared by using Student’s t-test (D, E, H, I); **P < 0.01.

CsOPR3 was chosen as the target gene to validate the rationality of the recommended RGs used in exogenous application of JA (Fig. 3D,H). When the best combination of the time interval from 3 h to 6 h, GAPDH1 and TIP41 (NF 1–2, F = 1.426, P = 0.028) was used for normalization, the expression level of CsOPR3 in JA-treated leaves was significantly higher than that in the control at 3 h, but no significant difference was found when normalized with the best combination of the time interval from 0.5 h to 1.5 h, CLATHRIN1 and UBC1 (NF 1–2, F = 0.163, P = 0.091) or 12 h to 48 h, CLATHRIN1 and GAPDH1 (NF 1–2, F = 0.599, P = 0.126) (Fig. 3D). When the most appropriate RG–GAPDH1 (NF 1, F = 0.023, P = 0.037) or the best combination of GAPDH1 and TIP41 (NF 1–2, F = 1.426, P = 0.028) of the time interval from 3 h to 6 h was used for normalization, the expression level of CsOPR3 in JA-treated leaves at 3 h was significantly higher than that in the control, but no significant difference was found when normalized with the combination of the two unstable RGs, TUA1 and ACTIN1 (NF 9–10, F = 0.138, P = 0.204), or with the most unstable RG (NF 10, F = 3.888, P = 0.259) (Fig. 3H).

CsPAL was chosen as the target gene to validate the rationality of the recommended RGs used in T. aurantii infestation (Fig. 3E,I). When the best combination at 48 h, ACTIN1 and EF1 (NF 1–2, F = 2.458, P = 0.047), was used for normalization, the expression level of CsPAL in treated leaves at 48 h was significantly higher than that in control, but no significant difference was found when normalized with the most stable combination of the time interval from 6 h to 24 h, ACTIN1 and UBC1 (NF 1–2, F = 2.921, P = 0.063) (Fig. 3E). When the most appropriate RG–ACTIN (NF 1, F = 0.116, P = 0.041) or the best combination of ACTIN1 and UBC1 (NF 1–2, F = 0.245, P = 0.030) of the time interval from 6 h to 24 h was used for normalization, the expression level of CsPAL in treated leaves at 6 h was significantly higher than that in control, but no significant difference was found when normalized with the most unstable combination of PTB1 and TUA1 (NF 9–10, F = 0.820, P = 0.141) or with the most unstable RG (NF 10, F = 2.355, P = 0.120) (Fig. 3I).

CsPALc was chosen as the target gene to validate the rationality of the recommended RGs used in E. onukii infestation (Fig. 3F,J). When the best combination of PTB1 and TBP at 96 h was used for normalization, the expression level of CsPALc at 96 h in pre-pregnant female-infested leaves was significantly higher than that of pregnant female-infested leaves (NF 1–2, F = 13.471, P = 0.002) and control leaves (F = 13.471, P = 0.008), but a relatively slight difference between pre-pregnant female-infested leaves and pregnant female-infested leaves was found when normalized with the combination of the two stable RGs in 12–72 h, GAPDH1 and UBC1 (NF 1–2, F = 4.838, P = 0.040) or in 120–144 h, TIP41 and EF1 (NF 1–2, F = 5.934, P = 0.018) (Fig. 3F). When the most appropriate RG–PTB1, or the most stable combination of PTB1 and TBP at 96 h was used for normalization, the expression level of CsPALc at 96 h in pre-pregnant female-infested leaves was significantly higher than that of pregnant female-infested leaves (NF 1, F = 10.566, P = 0.005; NF 1–2, F = 13.471, P = 0.002) and control leaves (NF 1, F = 10.566, P = 0.017; NF 1–2, F = 13.471, P = 0.008), but a relatively slight difference between pregnant female-infested leaves and pre-pregnant female-infested leaves was found when normalized with the most unstable combination, TIP41 and TUA1 (NF 9–10, F = 4.938, P = 0.037), and no significant difference was found when normalized with the most unstable RG (NF 10, F = 4.769, P = 0.072) (Fig. 3J).

CsOPR3 was chosen as the target gene to validate the rationality of the recommended RGs used in E. obliqua regurgitant treatment (Fig. 3G,K). When the best combination at 6 h, TBP and CLATHRIN1 was used for normalization, the expression level of CsOPR3 at 6 h in wounding leaves was significantly higher than that of regurgitant-treated leaves (NF 1–2, F = 32.921, P = 0.015) and intact leaves ((NF 1–2, F = 32.921, P = 0.000), but no significant difference between regurgitant-treated leaves and wounding leaves was found when normalized with the combination of the most two stable RGs in 1.5–3 h, TIP41 and TBP ((NF 1–2, F = 23.023, P = 0.051) or in 12–48 h, SAND1 and TBP (NF 1–2, F = 14.784, P = 0.176) (Fig. 3G). When the most appropriate RG–TBP (NF 1), or the most stable combination of TBP and CLATHRIN1 (NF 1–2) at 6 h was used for normalization, the expression level of CsOPR3 at 6 h in wounding leaves was significantly higher than that of regurgitant-treated leaves (NF 1, F = 26.647, P = 0.023; NF 1–2, F = 32.921, P = 0.015) and intact leaves (NF 1, F = 26.647, P = 0.001; NF 1–2, F = 32.921, P = 0.000), but no significant difference between regurgitant-treated leaves and wounding leaves was found when normalized with the most unstable combination, EF1 and TUA1 (NF 9–10, F = 7.557, P = 0.277) or with the most unstable RG (NF 10, F = 10.295, P = 0.117) (Fig. 3K).

Discussion

Normalizing results with one or more appropriate internal RGs is a simple and popular method for controlling error in qRT-PCR assays. To date, a few housekeeping genes have been rigorously identified and used as RGs in tea plants under abiotic stresses, such as cold, barrenness, drought, photoperiod and exogenous application of plant hormones (auxin, ABA, GA, IAA, MeJA and SA)25,26,28,32,33,34, leaf developmental stages and even different organs26,35. These results demonstrate that identifying appropriate RGs for target gene expression analysis under different experimental conditions is an essential prerequisite for developing a qPCR assay of tea plants. To the best of our knowledge, the present study is the first to define the proper RGs for qRT-PCR analysis in tea plants under infestations of different herbivorous pests and their related biotic stresses.

In the present study, ten candidate RGs were selected from those already identified as stably expressed RGs with high efficiency in tea molecular studies (Table 1). Previously, CsACINT1 was identified as one of the most unstable RGs under different experimental manipulations, such as different organs, cold or photoperiod treatment of leaves and shoots, diurnal expression in leaves, auxinole and lanolin treatment28. In the current study, our results showed that CsACINT1 was ranked as one of the five most unstable RGs for diurnal variation of different organs, JA-treated leaves, infestation of E. onukii, and mechanical damage plus E. obliqua regurgitant; however, this gene was determined as the best RG in T. aurantii infested leaves (Table 4). Similarly, CsACINT1 was found to be the most stably expressed RG in tea plants under Fe stress and in different organs33. CsUBC1 was identified as the most stable RG in almost all treatments, except for E. obliqua regurgitant treatment, while CsUBC1 was identified as the suitable RG when tea plants were under Mn stress24. CsTUA1 was ranked as the most unstable RG for tea plants across most of our experimental conditions, except for diurnal expression in stems (Table 4), while previous results revealed that CsTUA1 was the most stable RG for damage stresses of tea shoots. CsTBP was identified as one of the top two appropriate RGs for qRT-PCR analysis in hormonal stimuli tea leaf samples by GeNorm and NormFinder26, which includes ABA, GA, IAA, MeJA and SA. However, among the 10 RGs tested in this study, CsTBP was recommended as the seventh stable RG in JA stimuli samples, and CsGAPDH1 and CsCLATHRIN1 were recommended as the best RG combination for JA treatment (Table 4). The main reason for the difference is probably because different proposed RGs were adopted to rank the order. The results described above indicate, unsurprisingly, that no RG has been found to exhibit perfectly stable transcript accumulation in tea plants across different experimental conditions, even the already identified stable RGs.

The stability of the same RG varies with different plant species under diverse experimental conditions. TIP41-like protein (TIP41) was appraised as the best RG in different stages during development of bamboo (Phyllostachys edulis), reproductive stages of rapeseed (Brassica napus)36, and cucumber (Cucumis sativus) subjected to abiotic stresses and growth regulators37. Our results verified that TIP41 was the second most stable RG in JA-treated leaves in the time interval from 3 h to 6 h and the most stable RG in tea leaves infested by E. onukii in the time interval from 120 h to 144 h (Table 5). EF1 has been proven to be an appropriate RG for normalization of flower buds at different stages of female flower bud differentiation in the English walnut (Juglans regia)38, and EF1 was the second stable RG in tea leaves infested by E. onukii in the time interval from 120 h to 144 h or infested by T. aurantii at 48 h as well (Table 5). Similarly, EF1-a gene was found to perform well for aphid-infested chrysanthemum39, and EF1A 2a, EF1A 1a1 and EF1A 2b were also identified as the best RG in JA-treated leaves of soybean40. GAPDH, ACTIN and UBC are the commonly used RGs for qRT-PCR analysis in varied plant, whose function is maintaining cell survival irrespective of physiological conditions41,42,43. In this study, we found that ACTIN, UBC and GAPDH were the top three appropriate RGs for the whole samples of T. aurantii-infested leaves (Table 4), but GAPDH and ACTIN were less stable in peach44. CsUBC1 was also identified as an appropriate RG in almost all treatments, except for E. obliqua regurgitant treatment. HbUBC2a and HbUBC4 were identified as the most stable RGs in Brazilian rubber trees (Hevea brasiliensis) when all samples were analysed together45, but the UBC2 genes were not the proper RGs in soybean (Glycine max) and watermelon (Citrullus lanatus) exposed to cadmium or under abiotic stress46,47. Consequently, our results emphasize that the selection of reliable RGs for normalization under any given experimental design is a requirement for developing a proper qPCR assay.

Multiple RGs have been suggested for normalizing target gene expression, which will reduce the probability of biased normalization13,48. In the current study, our results demonstrated using multiple RGs simultaneously in qRT-PCR analysis would increase the sensitivity of gene expression in E. onukii infested leaves (Fig. 3J) or E. obliqua regurgitant treatment (Fig. 3K). Furthermore, our results suggest that if the processing time of treatment was long, the best RGs for normalization should be recommended according to the stability of the proposed RGs in different time intervals when intragroup differences were compared (Table 5; Fig. 3D–G), which would strongly increase the accuracy and sensitivity of target gene expression in tea plants under biotic stresses. However, when the differences of intergroup were compared, the RGs for normalization should keep consistent across different time points.

In summary, we screened a series of RGs to study the gene expression profile of different organs of tea plants with circadian rhythm, JA-treated tea leaves, tea leaves attacked by T. aurantii or E. onukii, and tea leaves treated with mechanical damage plus E. obliqua regurgitant. Our results provide a technical guidance for further study of the molecular mechanisms of tea plants under different biotic stresses.

Methods

Insects

The tea aphid (Toxoptera aurantii), the tea leafhooper (Empoasca onukii) and the tea looper (Ectropis obliqua) were caught from the experimental tea garden of the Tea Research Institute of the Chinese Academy of Agricultural Sciences (TRI, CAAS, N 30°10′, E 120°5′), Hangzhou, China. The insects were reared on the potted tea shoots in the controlled climate room at 26 ± 2 °C, 70 ± 5% rh, and a photoperiod of 14:10 h (L:D). Newly hatched larvae/nymphs were fed on tender tea shoots that were enclosed in net cages (75 × 75 × 75 cm) and kept in the room. After one generation, mixed age nymphs of T. aurantii were used for plant treatment. The 4th-instar E. onukii nymphs were collected individually and maintained in separate plastic tubes (1.5 cm wide × 9 cm high) with fresh tea stems, and then the newly molted adults were separated by sex according to morphological characteristics. One newly molted adult female and two males were kept in a plastic container (12 cm high × 7 cm diameter) with fresh tea shoots for 5 days to obtain a fully mated female. One-day-old virgin female adults were used as feeding adults, and 6-day-old fully mated females were used as pregnant females. Our biological bioassay results showed that the pre-oviposition period is 5 d, and 6-day-old fully mated females have similar food consumption to that of 1-day-old virgin females (unpublished data). Forth-instar larvae of E. obliqua were used for collecting regurgitants.

Regurgitant collection

As the method proposed by Yang et al.49, regurgitant was absorbed from E. obliqua oral cavity with a P200 Pipetteman (Gilson, Middleton, WI, USA). The collected regurgitant was homogenized at first. The homogeneous regurgitant was centrifuged for 5 min (10,000 × g), then the supernatant was collected and stored at −80 °C until use.

Tea plants and treatments

Longjing 43 tea plants (three-year-old) were used for experiments, which were planted individually in a plastic pot (14 cm diameter × 15 cm high), incubated in the greenhouse programmed at12-h photophase, 26 ± 2 °C, and 70–80% relative humidity. All materials were incubated under such conditions unless otherwise stated. Plants were fertilized with fertilizer once a month and irrigated once every other day. Day before processing, tea leaves were washed under the running water. Leaves in the same position but in different branches of the same tea plant were selected for each time points. Treatments were prepared as follows.

Different tissues in circadian rhythm

The second leaves (numbered sequentially from the most apically unfolded leaf down the stem), stems (tender internodes between the first and the second) and fibrous roots of tea plants were harvested every 2 h of a day in the autumn of 2018. Four replications were carried out.

Exogeneous application of JA

JA (Sigma Chemical Co., St. Louis, MO, USA) was dissolved in a small amount of ethylalcohol and made up to a concentration of 0.15 mg/mL in 50 mM sodium phosphate buffer (titrated with 1 M citric acid until pH 8). Treatments were individually sprayed with 8 mL of JA solution. Tea plants were individually sprayed with 8 mL of the buffer were used as control. Plants were treated at 10:00 am in the climate chamber. The second leaves were harvested at 0.5, 1.5, 3, 6, 12, 24 and 48 h after the start of treatment. Each treatment was replicated five times.

T. aurantii infestation

Fifty aphids were inoculated on the tender bud and the 1st leaves. A fine-mesh sleeve was used to cover the 2nd leaf to prevent aphid infestation and honeydew pollution. The second leaves that covered with mesh sleeves only were used as controls. The 2nd leaves were harvested at 6, 12, 24, 48 h after the start of treatment. Each treatment was replicated five times.

E. onukii infestation

The 2nd tender leaf was covered with a mesh sleeve into which 4 one-day-old virgin adult females or 4 six-day-old fully mated adult females that had been starved for 2 h were introduced at 9:00 pm. Plants with only their 2nd leaves covered with mesh sleeves were used as controls. Seventy-two hours after the start of treatments, E. onukii adults were carefully removed. Then, the 2nd leaves were harvested at 12, 24, 48, 72, 96, 120 and 144 h after the start of removal. Each treatment was replicated six times.

Mechanical damage plus E. obliqua regurgitant treatment

A fabric pattern wheel was used to damage tea leaves following the method described previously (2004)50. Each leaf was rolled 6 times, and 15 μL regurgitant was immediately painted to the puncture wounds. Deionized water in equal amounts was painted to the wounds for wounding treatment. The intact 2nd leaf was used as control. The treated and control 2nd leaves were harvested at 1.5, 3, 6, 12, 24 and 48 h after the start of treatment. Each treatment was replicated five times.

All treatments are briefly summarized below (Table 5).

Total RNA isolation, cDNA synthesis and qPCR analysis

The TRIzol™ kit (TIANGEN, Beijing, China) was used to isolate plant total RNA according to the protocol. The ratios of A260/280 and A260/230 of isolated RNA were examined by a spectrophotometer (Nanodrop ND 1000, Wilmington, DE, USA), and their ratios ranging from 2.0 to 2.2 and 2.0 to 2.3 individually suggested a high purity. One µg of total RNA was used to synthesize the first-strand cDNA by using a PrimerScript® RT Reagent Kit (Takara, Dalian, China) according to the protocol. A five gradient dilutions of cDNA was used as a template for each treatment to create the standard curves. After reverse transcription, the synthesized cDNA was stored at −20 °C until use.

Ten candidate RGs, including CsACTIN1, CsCLATHRIN1, CsEF1, CsGAPDH1, CsSAND1, CsTIP41, CsUBC1, CsPTB1, CsTUA1 and CsTBP, were chosen from previous reports for their high stability under different stresses of tea plant (Table 2). The qPCR reactions were carried out on a LightCycle® 480 Real-Time PCR System (Roche Diagnostics, Mannheim, Germany) with a 10-μl reaction system, which contains 0.5 μl forward and reverse primers (10 μM), 5 μl FastStart Essential DNA Green Master and 25 ng first-strand complementary DNA. The programs for all genes included a preliminary step at 95 °C for 10 min, 45 cycles of denaturation amplification at 95 °C for 15 s, at 60 °C for 15 s and at 70 °C for 12 s. Finally, a melting curve analysis from 60 °C to 95 °C was carried out to confirm the specificity of the PCR products. The standard curve method was used to calculate the gene relative expression level. Each sample was analyzed in triplicate.

Validation of selected reference genes

JA and SA signaling pathways play key roles in plant defense against herbivorous insects51,52, and JA and SA responsive genes could be expressed upon herbivore attack or hormone stimuli51,53. A key transcription factor of JA signaling–CsMYC2, a key enzyme in the biosynthesis of JA–CsOPR3, two enzyme involved in the biosynthesis of SA–CsPAL and CsPALc were selected as target genes to validate the rationality of diurnal expression in different tissues, JA treatment and E. onukii infestation, T. aurantii infestation or E. obliqua regurgitant treatment individually. RefFinder is a comprehensive tool, which was used to determine the geometric mean of genes. Based on the geometric mean of the genes, two different normalization factors (NFs) were the lowest and the highest mean values, and a single RG was the lowest or the highest mean value. Raw Ct values were transferred to relative quantities by the ΔΔCt method.

Data analysis

BestKeeper, geNorm, NormFinder, the ΔCt method and RefFinder were used to evaluate the stability of the candidate RGs. All the above methods can recommend the most stable RGs. While NormFinder, geNorm and the ΔCt method rely on transforming Ct values of (1 + E) ± ΔCt, original Ct values were used in RefFinder and BestKeeper. GeNorm software was used to identify the optimum number of RGs through the cut-off value. The Vn/n + 1 value means the pair-wise variation between two sequential NFs and the optimal number of RGs required for a perfect normalization. One-way ANOVA (Tukey’s test) was used to compare the differences among more than two treatments. The difference between two samples was analyzed by Student’s t-test.