Identification and evaluation of reliable reference genes for quantitative real-time PCR analysis in tea plants under differential biotic stresses

The selection of reliable reference genes (RGs) for normalization under given experimental conditions is necessary to develop an accurate qRT-PCR assay. To the best of our knowledge, only a small number of RGs have been rigorously identified and used in tea plants (Camellia sinensis (L.) O. Kuntze) under abiotic stresses, but no critical RG identification has been performed for tea plants under any biotic stresses till now. In the present study, we measured the mRNA transcriptional levels of ten candidate RGs under five experimental conditions; these genes have been identified as stable RGs in tea plants. By using the ΔCt method, geNorm, NormFinder and BestKeeper, CLATHRIN1 and UBC1, TUA1 and SAND1, or SAND1 and UBC1 were identified as the best combination for normalizing diurnal gene expression in leaves, stems and roots individually; CLATHRIN1 and GAPDH1 were identified as the best combination for jasmonic acid treatment; ACTIN1 and UBC1 were identified as the best combination for Toxoptera aurantii-infested leaves; UBC1 and GAPDH1 were identified as the best combination for Empoasca onukii-infested leaves; and SAND1 and TBP1 were identified as the best combination for Ectropis obliqua regurgitant-treated leaves. Furthermore, our results suggest that if the processing time of the 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, 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. The results of this study provide a technical guidance for further study of the molecular mechanisms of tea plants under different biotic stresses.


CsACTIN1
Different organs Nitrogen stress Fe stress Sun et al. 29 ; Liu et al. 20 ; Wang et al. 24 2 CsCLATHRIN1 Different organs Leaves with Cold and short photoperiod treatments Shoots after auxin antagonist auxinole treatments Hao et al. 28 3

CsEF1
Diurnal expression in leaves Hao et al. 28 4 CsGAPDH1 Different maturity of leaves Leaves with Cold and drought treatments Nitrogen stress Drought, cold, Al, and NaCl stresses Sun et al. 29 ; Ma et al. 25 ; Liu et al. 20 5

CsTIP41
In various tea leaf developmental stages Wu et al. 26 7

CsUBC1
Shoots with cold and short photoperiod treatments Mn stress Hao et al. 28 ; Wang et al. 24 8

CsTUA1
Physical damages Ma et al. 25 10 CsTBP In various tea leaf developmental stages Leaves with hormone treatments Mn stress Post-harvest leaves Posharvest Wu et al. 26 ; Wang et al. 24 ; Zhou et al. 27 Table 1. Ten housekeeping genes frequently used for qRT-PCR of tea plant.
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: (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).
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: 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:   (Table 3). According to the results of RefFinder, the stability ranking of RGs from the most to the least was as follows: 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: 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:  (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: 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).

T. aurantii infestation in the time interval from 6 h to 24 h
Recommended comprehensive ranking  www.nature.com/scientificreports www.nature.com/scientificreports/ 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).
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).

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 organs 26,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 treatment 28 . 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 organs 33 . 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 stress 24 . 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 NormFinder 26 , 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 regulators 37 . 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 chrysanthemum 39 , and EF1A 2a, EF1A 1a1 and EF1A 2b were also identified as the best RG in JA-treated leaves Scientific RepoRtS | (2020) 10:2429 | https://doi.org/10.1038/s41598-020-59168-z www.nature.com/scientificreports www.nature.com/scientificreports/ of soybean 40 . 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 conditions [41][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 peach 44 . 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 together 45 , but the UBC2 genes were not the proper RGs in soybean (Glycine max) and watermelon (Citrullus lanatus) exposed to cadmium or under abiotic stress 46,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 normalization 13,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. (2020) 10:2429 | https://doi.org/10.1038/s41598-020-59168-z www.nature.com/scientificreports www.nature.com/scientificreports/ 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 insects 51,52 , and JA and SA responsive genes could be expressed upon herbivore attack or hormone stimuli 51,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.