Selection of optimal reference genes for qRT-PCR analysis of shoot development and graviresponse in prostrate and erect chrysanthemums.

The prostrate cultivars of ground-cover chrysanthemum have been used in landscape gardening due to their small stature, large crown width and strong branching ability. qRT-PCR is a rapid and powerful tool for gene expression analysis, while its accuracy highly depends on the stability of reference genes. The paucity of authentic reference genes presents a major hurdle in understanding the genetic regulators of prostrate architecture. Therefore, in order to reveal the regulatory mechanism of prostrate growth of chrysanthemum stems, here, stable reference genes were selected for expression analysis of key genes involved in shoot development and graviresponse. Based on transcriptome data, eleven reference genes with relatively stable expression were identified as the candidate reference genes. After the comprehensive analysis of the stability of these reference genes with four programs (geNorm, NormFinder, BestKeeper and RefFinder), we found that TIP41 was the most stable reference gene in all of the samples. SAND was determined as a superior reference gene in different genotypes and during the process of shoot development. The optimal reference gene for gravitropic response was PP2A-1. In addition, the expression patterns of LA1 and PIN1 further verified the reliability of the screened reference genes. These results can provide more accurate and reliable qRT-PCR normalization for future studies on the expression patterns of genes regulating plant architecture of chrysanthemums.

Introduction

Quantitative real-time reverse-transcription polymerase chain reaction (qRT-PCR), a powerful tool for analyzing gene expression profiles, is commonly used in many research fields because of its high sensitivity, accuracy, specificity, throughput and low cost [1,2]. However, the accuracy of qRT-PCR analysis is highly dependent on an appropriate choice of reference genes [3]. The use of improper reference genes can lead to conflicting expression data [4,5]. The expression level of optimal reference genes for specific experimental system must be stable both in time and space, and should be unaffected by any treatment or genetic manipulation [6]. Whereas there are no universal reference genes that are stably expressed in all of the tissues and under all conditions [7]. EF1α (elongation factor 1α) was stably expressed during aphid infestation but with unstable expression during waterlogging stress, and PP2A (protein phosphatase 2A) was stably expressed under heat and waterlogging stress but the least stable reference gene for aphid infested plants [8]. MTP (metalloprotease) and ACT (actin) were the most stable in diploid and tetraploid Chrysanthemum nankingense, while PSAA (photosynthesis-related plastid gene representing photosystem I) and EF1α were the most stable in tetraploid and hexaploid C. zawadskii [6]. SAND (SAND family protein) was the most stable reference gene in floral developmental process in C. lavandulifolium [9]. Therefore, it is essential to earnestly evaluate and validate the stability of reference genes for every specific experimental design before expression analysis of target genes by qRT-PCR experiments.

Chrysanthemum has long been cultivated worldwide as a cut flower, in the garden and as a potted flower owing to its rich germplasm among ornamental plants [10-13]. Ground-cover chrysanthemum, a cultivar group of Chrysanthemum morifolium, is widely used in landscape gardening due to its large canopy, strong branching ability, large numbers of capitula, strong resistance and wide adaptability [14,15]. Cultivars of ground-cover chrysanthemum with vertical architecture are common, while the creeping or prostrate type is rare. If the prostrate architecture of C. yantaiense (abbreviated as YT), an outstanding trait with application value in landscape gardening, is introduced into ground-cover chrysanthemums, the ground-covered ability will be increased and greening costs will be reduced [16]. Most mutants with prostrate habit in model plants have been identified to be related to the loss of gravitropism or reduced gravitropism of the above-ground part. During graviresponse, UBQ (Ubiquitin) and ACT were primarily used as reference genes. UBQ was used for data normalization in qPCR analysis after gravity stimulation in Arabidopsis inflorescence stems [17] and rice shoot base [18,19]. ACT was used as an internal control to analyze relative expression level under graviresponse in peanut gynophores [20]. In the case of chrysanthemum, however, previous researchers have mainly focused on anatomical physiological characteristics of the creeping stem [21,22]. Few studies have reported expression patterns of gravitropic-related genes in chrysanthemum. When analyzing the expression pattern of gravitropic-related genes between upper and lower side of the creeping stem after gravistimulation, ACT was used as the housekeeping gene in chrysanthemum cultivar ‘Yuhuajinhua’ [23]. On account of no systematic evaluation of reference genes during shoot development and gravitropic response in chrysanthemum, it is imperative to identify stable reference genes before analyzing expression patterns of key genes involved in shoot development and graviresponse in prostrate and erect chrysanthemum.

In the present study, eleven reference genes, including four conventional reference genes (ACT, EF1α, GAPDH, and UBQ) and seven new ones (PP2A-1, PP2A-2, SAND, TIP41, PGK, MTP, and SKIP16), were selected based on model plants and chrysanthemum transcriptome libraries. The stability of these reference genes in shoot development and gravitropic response was then evaluated by geNorm, NormFinder, BestKeeper and RefFinder. So as to validate the effectiveness of the screened reference genes, the expression patterns of LA1 and PIN1 were analyzed. This study provides more accurate and reliable qRT-PCR normalization for future researches on expression patterns of genes regulating plant architecture of chrysanthemums.

Materials and methods

Plant materials

YT displays a prostrate growth habit with a height of less than 20 cm and the stem GSA (gravitropic setpoint angle) value of 90°–100°. While, ‘Fanhuasijin’ (CNA20090874, abbreviated to FH), a cultivar of ground-cover chrysanthemum, shows a vertical architecture with high stature of 60–70 cm and the stem GSA of 170°–180°. F1 and BC1 populations segregating for prostrate and erect architecture were constructed by crossing YT and FH. Sixteen F1 progenies with prostrate growth habit were selected to construct the ‘prostrate’ bulk and 16 F1 progenies with erect architecture were chosen to generate the ‘erect’ bulk. Uniformly-sized cuttings from each selected strain were planted into a 1:1 mixture of vermiculite and perlite. Rooted cuttings were then transplanted into individual pots (10*8.5*9.5 cm) with a 1:1 mixture of peat and perlite. To achieve the great consistency and obtain the gene expression profile of the early stage, the terminal buds were removed uniformly two weeks after transplanting. Based on the morphologic changes of the prostrate bulk from preliminary experiment, the consecutively developmental process of stem was divided into three stages: Stage I (1 week after cutting the terminal bud, GSA of stem from prostrate bulk was 150°–180°), Stage II (3 weeks after cutting the terminal bud, GSA of stem from prostrate bulk was 105°–150°) and Stage III (5 weeks after cutting the terminal bud, GSA of stem from prostrate bulk was 90°–105°). Thus, the stems of prostrate and erect bulk at three stages were collected as set 1 in this study (Fig 1A). Set 2 consisted of 16 prostrate and erect strains at Stage III (Fig 1B). Set 3 was composed of eight new cultivars of ground-cover chrysanthemum from BC1 population, including four prostate cultivars (‘Fukanbowu’ ‘Fukannongyun’ ‘Fukanfendai’ and ‘Fukanhongxiu’) and four vertical cultivars (‘Beilinqiuyun’ ‘Fukanchenlu’ ‘Fukanchiyan’ and ‘Fukanxiaoyue’) at Stage III (Fig 1C). The cutting seedlings were maintained in a greenhouse under long-day photoperiod (16 h light/8 h dark), the day/night temperature of 25/18°C, and the relative humidity of 70%.

Fig 1

Different developmental stages, genotypes and gravitropic responses of prostrate and erect chrysanthemums.

(A) Set 1: shoots of prostrate and erect bulk at three developmental stages. (B) Set 2: sixteen prostrate and erect strains. (C) Set 3: eight new cultivars of ground-cover chrysanthemum. (D) Set 4: gravitropic responses of stems of YT at different time points (0, 3, 6, 9, 12, 24, and 48 h). (E) Set 5: gravitropic responses of stems of FH at different time points (0, 3, 6, 9, 12, 24, and 48 h). I, Stage I; II, Stage II; III, Stage III; P, prostrate; E, erect; BW, ‘Fukanbowu’, CNA20170990; NY, ‘Fukannongyun’, CNA20170985; FD, ‘Fukanfendai’, CNA20170984; HX, ‘Fukanhongxiu’, CNA20170986; QY, ‘Beilinqiuyun’, CNA20170987; CL, ‘Fukanchenlu’, CNA20170988; CY, ‘Fukanchiyan’, CNA20170989; XY, ‘Fukanxiaoyue’, CNA20170983.

Cutting-seedlings of YT and FH at the 10–12 leaf stage were subjected to gravistimulation by gently rotating the pot 90° in a growth chamber in the dark with the temperature of 22±1°C and the relative humidity of 60% [18]. Stems of FH displayed obviously asymmetric growth at 3 h after horizontal treatment and regained vertical growth in less than 12 h. However, shoots of YT began to bend at 6 h and it took more than 24 h for YT to completely recover upright growth. Segments from the curved internode were separately collected at 0, 1, 3, 6, 9, 12, 24, and 48 h. Then, these stem segments were rapidly divided into the upper and lower half by cutting along the mid longitudinal axis of the stem [23]. Dissected segments from YT and FH at eight time points composed set 4 (Fig 1D) and set 5 (Fig 1E), respectively. All of the samples were collected in triplicates, immediately frozen in liquid nitrogen and then stored at −80°C.

RNA extraction, DNase I digestion and first-strand cDNA synthesis

Total RNA from all samples was isolated using E.Z.N.A Plant RNA Kit (OMEGA BIO-TEK, Norcross, USA) according to the manufacturer’s instruction. Purity and concentration of total RNA were detected by a UV-Vis Spectrophotometer Q5000 (Quawell, San Jose, USA). The RNA samples showing an A260/A280 ratio of 1.9−2.2 and an A260/A230 ratio greater than 2.0 were used for subsequent experiments. The integrity of RNA was further assessed by 1% (w/v) agarose gel electrophoresis. RNA samples with a 28S/18S ratio of 1.5–2.0 and without smears were used for further analysis. 1.0 μg total RNA was reverse-transcribed by the FastQuant RT Kit (with gDNase) (TIANGEN, Beijing, China) according to operating manual. All of the synthesized cDNAs were stored at −20°C until PCR experiments.

Selection of candidate reference genes and primer designing

The candidate reference genes were derived from model plants and related species. Those nucleotide sequences were used as query sequences to search the Chrysanthemum morifolium transcriptome libraries (NCBI SRA accession: SRP173747) using the TBLASTX program. Several sequences with high similarity (e-value < 10−9) were obtained. From these sequences, a number of sequences stably expressed (|log2Ratio| < 0.3) were screened as candidates. The primers of those candidates were designed using Primer Premier 5.0 and shown in S2 Table. PCR products of reference genes were purified and then cloned into pMD18-T vector (TaKaRa, Japan). The positive clones were sequenced by Sangon Biotech Co. Ltd. (Sangon, China). Specific primers of candidate reference genes for qRT-PCR were designed using PrimerQuest® tool (https://sg.idtdna.com/Primerquest/Home/Index) with primer Tm of 55−65°C, primer length of 19−25 bp and short PCR products of 100−250 bp.

Test of specificity and amplification efficiency of primers

The performance of the primers was evaluated by 1% (w/v) agarose gel electrophoresis and qRT-PCR. Each amplicon was purified using TIANgel Midi Purification Kit (TIANGEN, China) and cloned into pMD18-T vector (TaKaRa, Japan). The positive clones were sequenced by Sangon Biotech Co. Ltd. (Sangon, China). The specificity of each primer pair was accepted when it conform to the following criteria: (1) Only a specific product was generated for cDNA template, (2) no product was generated for genomic DNA template, and (3) the melting curve of qRT-PCR showed a single peak. Amplification efficiency (E) and correlation coefficient (R2) of each primer pair were determined by the slope of a standard curve generated from serial dilutions (×1, ×5, ×25, and ×125) of pooled cDNA samples as the template [2].

qRT-PCR assays

qRT-PCR experiments were conducted using the SYBR Premix Ex Taq II Kit (TaKaRa, Japan) based on a PikoReal Real-time PCR system (Thermo Fisher Scientific, USA). Each reaction was prepared in a total volume of 20.0 μl containing 2.0 μl of diluted cDNA (~20 ng), 0.8 μl of each primer (10 μM, Sangon Biotech, Shanghai), 10.0 μl of 2 × SYBR Premix Ex Taq II, and 6.4 μl of double-distilled water. The amplification in all reactions was performed according to following program: 95°C for 30 s, 40 cycles of 95°C for 5 s and 62°C for 30 s, 72°C for 30 s. Melting curve was recorded at the end of the qPCR by heating from 60°C to 95°C with 0.2°C increment every 1 s. Negative controls with no-RT RNA were conducted to ensure there was no genomic contamination. No template controls were also conducted with double-distilled water as template. Each reaction was carried out with three technical replicates.

Stability analysis of reference genes

The expression levels of the candidate reference genes in all of the samples were determined by quantification cycle (Cq) values. Box plots of raw Cq values were generated using SPSS v22.0 software (IBM, Chicago, IL, USA) to exhibit variation of each reference gene in all tested samples.

The expression stability of 11 candidate reference genes was assessed using four computational programs, including geNorm [24], NormFinder [25], BestKeeper [26] and RefFinder (http://150.216.56.64/referencegene.php). For geNorm and NormFinder analysis, Cq values were converted into relative quantities by the 2-ΔCq method, in which ΔCq = each corresponding Cq value—minimum Cq value [27]. The geNorm program calculates stability measure (M) of gene expression based on the average pairwise variations (V) of a particular gene against all other reference genes. The lowest M value represents the most stable expression, and M value < 0.5 is taken as an indicator of stable expression [28]. NormFinder program (V0953, Aarhus, Denmark) calculates intra- and inter-group variations, and then combines the two results into a stability value of each candidate gene. Genes with lowest stability value are the most stable. BestKeeper (Version1.0, Munich, Germany) ranks the expression stability of reference genes according to the coefficient of variance (CV) and the standard deviation (SD) of the Cq values, and the reference gene with lowest CV and SD is identified as the most stable one. RefFinder is an online tool that integrates four methods (geNorm, NormFinder, BestKeeper, and ΔCT) to compare and rank the stability of candidate reference genes comprehensively.

Validation of selected reference genes

LAZY1 (LA1) is the principal determinant of branch angle and mediates plant architecture [29,30,18,31-34]. PIN1, one of the PIN-FORMED (PIN) family members, encodes an auxin efflux carrier involved in the auxin redistribution in gravitropic response [35,36]. The primers of LA1 and PIN1 for qRT-PCR are presented in Table 1. The transcription of these genes was quantified using both stable and unstable reference genes to validate the reliability of the results. The expression levels of these target genes were calculated using the 2-ΔΔCq method [27].

Table 1

Primer sequences and amplicon characteristics of target genes.

Gene symbol	GenBank ID	Primer sequence (5’–3’)Forward/Reverse	PCR efficiency (%)	Regression coefficient (R2)	Amplicon length(bp)	Melting TM(°C)
LA1	MK381414	GCAACATTCCACAGGCTACA/	90.1	0.995	95	82.7
		CAGCTCCAACACCAGGTAATC
PIN1	MK381416	GAAGTGGTGTGAGTCCAAGAA/	106.3	0.996	95	79.8
		CCCTCCTTGACTACCATTAACC

Results

Identification of candidate reference genes from transcriptome libraries

Based on previous studies on selection of reference genes in Arabidopsis [3,37] and chrysanthemums [6,8,9,38,39], fifteen reference genes were ubiquitously expressed in six different transcriptome libraries of C. morifolium and passed the BLAST test. Eleven of these genes, including PGK (phosphoglycerate kinase), MTP (metalloprotease), PP2A-1 (protein phosphatase 2A-1), PP2A-2 (protein phosphatase 2A-2), ACT (actin), EF1α (elongation factor 1α), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), TIP41 (TIP41-like family protein), UBQ (ubiquitin extension protein), SAND (SAND family protein) and SKIP16 (SKP1/ASK-interacting protein 16), showed the most stable expression with |log2Ratio| < 0.3 (S1 Table). Therefore, they were selected as candidates for normalizing the gene expression during shoot development and gravitropic response. Moreover, all candidate reference genes were cloned and sequenced. The corresponding NCBI accession numbers are given in Table 2.

Table 2

Primer sequences and amplicon characteristics of candidate reference genes.

Gene symbol	GenBank ID	Primer sequence (5′–3′) for qPCRForward/Reverse	PCR efficiency (%)	Regression coefficient (R2)	Amplicon length(bp)	Melting TM(°C)
PGK	MK381403	CGTTGGTTATTCTTGTATGTGGC/	102.8	0.995	175	83.2
		CTGAAGTCTCGTGCCCATATAG
MTP	MK381404	GATTAAAGCCAACAGTCTTGCG/	94.3	0.999	158	84.8
		ACGTTCCAAGTATCTCAATCCTG
PP2A-1	MK381405	TTGGCGGATATGGTGATTAGG/	96.6	0.997	109	79.9
		GTTGTGTTGCTTCAAGAACCTC
PP2A-2	MK381406	ATCAGAACAGGAGGTCAGG/	102.9	0.997	171	82.8
		TAATTTGTATCGGGGCACTT
ACT	MK381407	AGCCGTTCTTTCCCTGTATG/	106.9	0.997	186	83.8
		GAATACCCACGCTCTGTAAGG
EF1α	MK381413	CAATTGCTAAACCATCTGCCG/	103.8	0.994	230	82.0
		AGGCTTGAACTGTGAACGAG
GAPDH	MK381408	CCCACCTTCTCAAATACGACTC/	98.4	0.999	175	84.4
		ACTCCTGTCCCTTCTATCACC
TIP41	MK381409	TGAGTTGGCTGATAATGGAGTC/	98.0	0.999	172	82.1
		TCGGACTATAACAGGCTTTGC
UBQ	MK381410	ATCTTGTGTTGAGGTTGAGGG/	105.4	0.999	246	83.3
		GGAGACGAAGGACAAGATGAAG
SAND	MK381411	TTACCTGTTGACCCATCTGC/	106.6	0.994	135	82.3
		CATAAAGTGCCAAAGTCCAGC
SKIP16	MK381412	AGCTGTCTGAACCTGTTGATC/	98.9	0.997	139	82.2
		CTTCTGATTCTGTCCCAAACG

Performance of primers for each reference gene

The primer sequences and amplicon length of eleven reference genes are presented in Table 2. The performance of the primers was evaluated by agarose gel electrophoresis and qRT-PCR. The agarose gel electrophoresis results showed that all eleven primer pairs amplified a single product of expected size (Fig 2). The sequences of amplicons were almost identical (similarity of 98%–100%) to that of transcriptome data of C. morifolium. Melting curve analysis was performed by qRT-PCR after 40 cycles of amplification, and the presence of a single peak in melting curve further confirmed the specificity of all tested primer pairs (S1 Fig). The amplification efficiencies of these primers ranged from 94.3% to 106.9% and the linear standard curve for each gene from a five-fold serial dilution of cDNA showed R2 > 0.99 (Table 2 and S2 Fig). These results reflected the high specificity and quality of the qRT-PCR, thus these eleven pairs of primer of candidate genes were used in the further assay.

Fig 2

Agarose gel electrophoresis of amplicons of eleven candidate reference genes.

Marker (M) from top to bottom: 500, 250 and 100 bp. The primer pairs in each separation are ordered from left to right as PGK, MTP, PP2A-1, PP2A-2, ACT, EF1α, GAPDH, TIP41, UBQ, SAND, SKIP16 (two technical replicates).

Expression profiles of candidate reference genes

Quantification cycle (Cq) values of candidate reference genes are shown in boxplots (Fig 3). The Cq values ranged from 17.77 to 28.88 across all samples. Among the eleven candidate reference genes, ACT showed the lowest expression level with a mean Cq value of 27.13, while UBQ displayed the highest expression level with a mean Cq value of 19.25 (Fig 3F). The raw Cq values of candidate reference genes fluctuated in varying degrees in each set and across all samples, indicating that none of these genes had a constant expression level (Fig 3A–3F). Preliminary analysis of raw Cq values with boxplots could not provide enough information on expression stability. Thus, four computational programs were used to further evaluate the stability of reference genes for chrysanthemum in different shoot development stages and gravitropic response.

Fig 3

The Cq values of the eleven candidate reference genes depicted with boxplots.

(A) Set 1: stems of ‘prostrate’ and ‘erect’ bulk at three stages. (B) Set 2: stems of sixteen prostrate and erect strains at Stage III. (C) Set 3: stems of eight new cultivars of ground-cover chrysanthemum. (D) Set 4: gravitropic responses of stems of YT at different time points (0, 1, 3, 6, 9, 12, 24, and 48 h). (E) Set 5: gravitropic responses of stems of FH at different time points (0, 1, 3, 6, 9, 12, 24, and 48 h). (F) The Cq values of the eleven candidate reference genes across all samples. The upper and lower edges of the box represent the upper and lower quartiles, and the middle black line is the median. The upper and lower whiskers depict the smallest and largest unbooked values. Circles and asterisks indicate mild and extreme outliers, respectively.

The stability of candidate reference genes

We analyzed the stability of candidate reference genes among each set individually and the overall stability among total samples. Four computational programs, including geNorm, NormFinder, BestKeeper, and RefFinder, were used to evaluate the stability of candidate reference genes.

geNorm analysis

The geNorm program was used to rank the expression stability of candidate reference genes by calculating the average expression stability (M) [24]. A lower M value indicates a more stable gene expression, and the M value should be lower than 0.5 [28]. The ranking orders of each set based on the M value are depicted in Fig 4A. Most genes were stable with M value below 0.5, except for EF1α (M = 0.549) in set 2. In set 1, SAND and GAPDH were the most stable genes with the lowest M value of 0.067, while SKIP16 was the least stable one with the highest M value of 0.164. The most stably expressed genes in set 2 were SAND and PGK with an M value of 0.207, whereas EF1α performed poorly. In set 3, SAND and PGK were the most highly ranked with an M value of 0.247, while the least stable gene was UBQ. TIP41 and PP2A-1 performed best both in set 4 and set 5, whereas ACT and GAPDH had the largest M value, indicating that these genes were the least stably expressed. In conclusion, TIP41 and SAND showed the highest stability among all samples, and GAPDH was the least stable reference gene.

Fig 4

Average expression stability (M) and pairwise variation (V) of the candidate reference genes calculated by geNorm.

(A) Average expression stability values M after stepwise exclusion of the least stable genes across all sets. (B) Determination of optimal number of reference genes required for accurate normalization.

geNorm also calculates the pairwise variation (V) to determine the optimal number of reference genes for reliable normalization. If the Vn/n+1 value is lower than the threshold of 0.15, the minimum number of the most suitable reference genes is n. The V2/3 values for set 1 to set 5 were far lower than 0.15, indicating that one stable reference gene is sufficient to obtain accurate results. The V2/3 value (0.138) of the total set demonstrated that two reference genes were suitable for normalization of all the samples (Fig 4B).

NormFinder analysis

The ranking orders of candidate reference genes in all sets determined by stability value calculated using NormFinder are shown in Table 3, with lower value indicating higher stability. The ranking orders determined by this method were similar to the results generated by geNorm. TIP41 was the most stable genes among all tested samples both in geNorm and NormFinder analysis. For set 1 and set 2, the top two most stable and the least stable reference genes generated by NormFinder were highly consistent with those determined by geNorm. In set 3 and set 4, the three most stable genes and the two least stable gene generated by NormFinder were almost the same as those generated by geNorm. NormFinder determined PGK and EF1α as the two most stable genes in sample set 5, whereas PGK and EF1α were ranked sixth and fifth by geNorm, respectively.

Table 3

Ranking order of candidate reference genes and stability values calculated by NormFinder.

Rankingorder	Set 1		Set 2		Set 3		Set 4		Set 5		Total
	Gene	Stability	Gene	Stability	Gene	Stability	Gene	Stability	Gene	Stability	Gene	Stability
1	SAND	0.050	PGK	0.073	TIP41	0.073	PGK	0.034	PGK	0.041	TIP41	0.141
2	GAPDH	0.056	SAND	0.137	SAND	0.141	PP2A-1	0.088	EF1α	0.041	PP2A-2	0.185
3	PP2A-2	0.057	PP2A-2	0.168	PGK	0.159	TIP41	0.090	TIP41	0.051	PGK	0.233
4	PGK	0.062	PP2A-1	0.205	SKIP16	0.160	EF1α	0.105	PP2A-1	0.069	MTP	0.292
5	PP2A-1	0.086	TIP41	0.246	EF1α	0.196	SAND	0.124	MTP	0.079	SAND	0.304
6	ACT	0.090	GAPDH	0.248	GAPDH	0.221	MTP	0.139	PP2A-2	0.094	EF1α	0.366
7	EF1α	0.092	MTP	0.264	PP2A-2	0.227	PP2A-2	0.147	UBQ	0.126	PP2A-1	0.399
8	TIP41	0.092	SKIP16	0.331	PP2A-1	0.249	SKIP16	0.167	SAND	0.132	SKIP16	0.463
9	UBQ	0.098	UBQ	0.346	ACT	0.261	UBQ	0.211	SKIP16	0.176	UBQ	0.540
10	MTP	0.098	ACT	0.354	MTP	0.270	ACT	0.449	GAPDH	0.546	ACT	0.669
11	SKIP16	0.102	EF1α	0.484	UBQ	0.325	GAPDH	0.480	ACT	0.734	GAPDH	1.351

BestKeeper analysis

BestKeeper ranks the candidate reference genes based on the standard deviation (SD) and the coefficient of variation (CV) of their Cq values [26]. Lower SD and CV values represent higher stability. If the SD values of reference genes are larger than 1.0, these reference genes are considered to be unsuitable for gene expression normalization. Our results showed that almost all the reference genes had SD values smaller than 1.0 except GAPDH in all samples (Table 4). The ranking orders were different from the results calculated by GeNorm and NormFinder. For set 1, SKIP16 was the most stable gene by BestKeeper, while in geNorm and NormFinder analysis, SKIP16 was the least stable gene. In set 3, UBQ ranked second, whereas it ranked at the bottom both in geNorm and NormFinder analysis. For set 2, set 4, set 5 and total, the top three most stable and the least stable reference genes determined by BestKeeper were similar to those generated by geNorm and NormFinder (Fig 4A, Tables 3 and 4).

Table 4

Ranking order of candidate reference genes and their SD calculated by BestKeeper.

Rankingorder	Set 1		Set 2		Set 3		Set 4		Set 5		Total
	Gene	CV±SD	Gene	CV±SD	Gene	CV±SD	Gene	CV±SD	Gene	CV±SD	Gene	CV±SD
1	SKIP16	0.79±0.20	SAND	0.53±0.12	TIP41	0.96±0.23	PP2A-1	0.50±0.12	PP2A-1	0.60±0.14	TIP41	1.01±0.23
2	ACT	0.75±0.21	PP2A-1	0.64±0.15	UBQ	1.37±0.28	SAND	0.53±0.13	TIP41	0.63±0.15	SAND	0.97±0.23
3	PP2A-1	0.87±0.21	PGK	0.69±0.17	GAPDH	1.47±0.30	TIP41	0.57±0.13	UBQ	1.07±0.20	MTP	1.58±0.37
4	PGK	1.01±0.26	TIP41	1.21±0.27	PP2A-1	1.20±0.31	PGK	0.62±0.16	SAND	0.84±0.20	PP2A-2	1.56±0.39
5	MTP	1.14±0.28	PP2A-2	1.28±0.31	SAND	1.37±0.35	MTP	0.75±0.17	MTP	0.96±0.22	PP2A-1	1.66±0.40
6	UBQ	1.40±0.28	UBQ	1.68±0.31	MTP	1.49±0.36	SKIP16	0.82±0.20	EF1α	1.00±0.26	PGK	1.84±0.47
7	SAND	1.16±0.30	SKIP16	1.37±0.33	EF1α	1.47±0.39	PP2A-2	0.84±0.21	PGK	1.11±0.29	UBQ	2.52±0.48
8	PP2A-2	1.19±0.30	ACT	1.21±0.33	PP2A-2	1.56±0.39	EF1α	0.91±0.24	SKIP16	1.16±0.30	EF1α	2.23±0.59
9	GAPDH	1.54±0.32	GAPDH	1.91±0.34	SKIP16	1.58±0.39	UBQ	1.70±0.33	PP2A-2	1.33±0.33	SKIP16	3.05±0.76
10	EF1α	1.28±0.34	MTP	1.52±0.35	PGK	1.59±0.42	ACT	1.92±0.51	GAPDH	2.66±0.62	ACT	3.48±0.94
11	TIP41	1.56±0.38	EF1α	2.56±0.66	ACT	1.66±0.45	GAPDH	2.46±0.53	ACT	3.27±0.90	GAPDH	8.43±1.76

RefFinder analysis

Finally, RefFinder was used to generate a comprehensive evaluation of candidate reference genes by integrating three distinct algorithms (geNorm, NormFinder and BestKeeper) [40]. The ranking orders determined by RefFinder are depicted in Table 5. Though some differences in ranking orders were found among four programs, the most stable genes were roughly identical. SAND was found to be the most stably expressed gene in set 1. PGK ranked the first in set 2, and PP2A-1 ranked the first in set 4. TIP41 was the most stable gene in set 3, set 5 and total samples. The least stable reference genes in set 1 and set 2 were UBQ and EF1α, respectively. ACT was the least stably expressed gene in set 3 and set 5, as well as GAPDH in set 4 and total samples.

Table 5

Stability ranking of eleven candidate reference genes based on RefFinder analysis.

Ranking	1	2	3	4	5	6	7	8	9	10	11
Set 1
GeNorm	GAPDH | SAND		EF1α	PP2A-2	TIP41	PGK	MTP	PP2A-1	UBQ	ACT	SKIP16
NormFinder	SAND	GAPDH	PP2A-2	PGK	PP2A-1	EF1α	TIP41	ACT	MTP	UBQ	SKIP16
BestKeeper	SKIP16	ACT	PP2A-1	PGK	UBQ	MTP	SAND	PP2A-2	GAPDH	EF1α	TIP41
Recommended comprehensive ranking	SAND	GAPDH	PP2A-2	PGK	PP2A-1	EF1α	ACT	SKIP16	TIP41	MTP	UBQ
Set 2
GeNorm	PGK | SAND		PP2A-1	TIP41	PP2A-2	GAPDH	MTP	UBQ	SKIP16	ACT	EF1α
NormFinder	PGK	SAND	PP2A-2	PP2A-1	TIP41	GAPDH	MTP	SKIP16	UBQ	ACT	EF1α
BestKeeper	SAND	PP2A-1	PGK	TIP41	PP2A-2	UBQ	SKIP16	ACT	GAPDH	MTP	EF1α
Recommended comprehensive ranking	PGK	SAND	PP2A-1	PP2A-2	TIP41	GAPDH	MTP	UBQ	SKIP16	ACT	EF1α
Set 3
GeNorm	PGK | SAND		TIP41	SKIP16	PP2A-2	EF1α	GAPDH	ACT	PP2A-1	MTP	UBQ
NormFinder	TIP41	SAND	PGK	SKIP16	EF1α	GAPDH	PP2A-2	PP2A-1	ACT	MTP	UBQ
BestKeeper	TIP41	UBQ	GAPDH	PP2A-1	SAND	MTP	EF1α	PP2A-2	SKIP16	PGK	ACT
Recommended comprehensive ranking	TIP41	SAND	PGK	SKIP16	GAPDH	EF1α	PP2A-2	PP2A-1	UBQ	MTP	ACT
Set 4
GeNorm	PP2A-1 | TIP41		PGK	PP2A-2	SAND	EF1α	MTP	SKIP16	UBQ	ACT	GAPDH
NormFinder	PGK	PP2A-1	TIP41	EF1α	SAND	MTP	PP2A-2	SKIP16	UBQ	ACT	GAPDH
BestKeeper	PP2A-1	SAND	TIP41	PGK	MTP	SKIP16	PP2A-2	EF1α	UBQ	ACT	GAPDH
Recommended comprehensive ranking	PP2A-1	PGK	TIP41	SAND	EF1α	MTP	PP2A-2	SKIP16	UBQ	ACT	GAPDH
Set 5
GeNorm	PP2A-1 | TIP41		UBQ	MTP	EF1α	PGK	SAND	PP2A-2	SKIP16	GAPDH	ACT
NormFinder	PGK	EF1α	TIP41	PP2A-1	MTP	PP2A-2	UBQ	SAND	SKIP16	GAPDH	ACT
BestKeeper	PP2A-1	TIP41	UBQ	SAND	MTP	EF1α	PGK	SKIP16	PP2A-2	GAPDH	ACT
Recommended comprehensive ranking	TIP41	PP2A-1	EF1α	PGK	UBQ	MTP	SAND	PP2A-2	SKIP16	GAPDH	ACT
Total
GeNorm	TIP41 | SAND		PP2A-2	MTP	PP2A-1	UBQ	PGK	EF1α	SKIP16	ACT	GAPDH
NormFinder	TIP41	PP2A-2	PGK	MTP	SAND	EF1α	PP2A-1	SKIP16	UBQ	ACT	GAPDH
BestKeeper	TIP41	SAND	MTP	PP2A-2	PP2A-1	PGK	UBQ	EF1α	SKIP16	ACT	GAPDH
Recommended comprehensive ranking	TIP41	PP2A-2	SAND	MTP	PGK	PP2A-1	EF1α	UBQ	SKIP16	ACT	GAPDH

Validation of selected reference genes by expression analysis of target genes

Previous studies have shown that LAZY1 (LA1) is the principal determinant of branch angle that mediate plant architecture [18,29,30,31-34]. PIN1, one of the PIN-FORMED (PIN) family members, encodes an auxin efflux carrier involved in the auxin redistribution in gravitropic response [35,36]. Melting curve and amplification efficiencies of these target genes are shown in S3 Fig. To verify the stability of the screened reference genes, the expression patterns of these target genes were analyzed using two most stable and one least stable reference genes based on the results of four programs. In set 1, when normalized using the most stable reference genes SAND, the relative expression level of CmLA1 increased in the erect bulk during the process of shoot development, while it basically remained unchanged in prostrate bulk among three developmental stages. The expression pattern of CmLA1 normalized with the second most stable reference gene GAPDH displayed similar trends with that normalized with SAND. However, when normalized by the least stable gene UBQ, the expression profile showed significant differences (Fig 5A). In set 2 and set 3, when normalized with the two most stable reference genes, the expression patters of CmLA1 were almost identical, with the transcript level three times higher in vertical strains and cultivars as compared to prostrate plants. Whereas, the expression pattern exhibited obvious discrepancies when normalized by the least stable gene (Fig 5B and 5C). During the process of gravitropic response in YT, when normalized with two most stable reference genes PP2A-1 and PGK, similar expression patterns were generated that relative expression levels of PIN1 in the upper and lower half of stem were almost the same. While, the transcription profile of PIN1 was contradictory to the above patterns when normalized using GAPDH, the least stable gene (Fig 5D). In FH (a vertical cultivar of ground-cover chrysanthemum ‘Fanhuasijin’, CNA20090874), the expression pattern of PIN1 normalized with TIP41, when subjected to gravistimulation, was almost the same as that normalized by PP2A-1. The transcript abundance of PIN1 in the lower half of stem was significantly higher than that in the upper half at 3, 9, 12, and 24 h after horizontal treatment. Nevertheless, no obvious differences were observed between the upper and lower half of stem when normalized by the least stable gene ACT (Fig 5E). Taken together, the significant discrepancies in expression level of PIN1 might be responsible for the differences in gravitropic response between YT and FH. These results verified the reliability of screened reference genes and revealed that the accuracy of qRT-PCR analysis could be altered when using the least stable reference gene.

Fig 5

The expression profiles of two target genes normalized by stable and unstable reference genes.

(A) The expression pattern of LA1 of prostrate and erect bulk at three different shoot development stages; (B) Expression levels of LA1 in prostrate and erect strains; (C) Expression levels of LA1 in new cultivars of ground-cover chrysanthemum; (D) The expression pattern of PIN1 during the process of gravitropic response in YT; (E) The expression pattern of PIN1 during the process of gravitropic response in FH. I, Stage I; II, Stage II; III, Stage III; P, prostrate; E, erect; HX, ‘Fukanhongxiu’, CNA20170986; BW, ‘Fukanbowu’, CNA20170990; NY, ‘Fukannongyun’, CNA20170985; CY, ‘Fukanchiyan’, CNA20170989; QY, ‘Beilinqiuyun’, CNA20170987; CL, ‘Fukanchenlu’, CNA20170988; L, the lower half of stem; U, the upper half of stem. Error bars represent SD, n = 3. *, P < 0.05, **, P < 0.01, Student t test.

Discussion

The significance of systematic evaluation of reference genes in shoot development and gravitropic response in chrysanthemum

The qRT-PCR has become a forceful method to analyze gene expression, however its accuracy mainly relys on suitable reference genes [1-3]. Previous studies had demonstrated that no universal reference genes could express steadily in various organisms and circumstances and that contradictory results could be generated while using unsuitable reference genes [4,5,7]. Several genes had been evaluated as suitable reference genes for chrysanthemum in previous studies focused on biotic and abiotic stress [8], photoperiodic treatments [38], cross-ploidy level comparisons [6], flower color [39], and floral development [9]. However, no detailed verification has been conducted on whether these reference genes were suitable for data normalizing during shoot development and gravitropic response. The traditional housekeeping gene ACT, which has been considered as a reliable reference gene for many years, was used for data normalization in real-time PCR during shoot graviresponse in chrysanthemum cultivar ‘Yuhuajinhua’ [23]. Whereas, in this research, ACT had poor performance in most of the sample sets, which was consistent with studies in chrysanthemum cultivar ‘Zhongshanzigui’ [8], Chrysanthemum lavandulifolium [38], and Lagerstroemia speciose [41]. Therefore, a systematic evaluation of reference genes is extremely important before analyzing the expression patterns of key genes involved in shoot development and gravitropic response in prostrate and erect chrysanthemum.

The high efficiency of candidate reference gene selection based on transcriptome data

The commonly-used approaches for selecting candidate reference genes are literature reviewing and public database searching for housekeeping genes related to certain biological processes such as glycolysis, cellular metabolism, protein synthesis and degradation [42]. Several studies have shown that the expression of these reference genes might vary to a great extent under different experimental conditions [43]. Thus, unreliable results may be generated when blindly selecting these reference genes as candidates without any other expression data [7,24]. Expression data of thousands of genes can be obtained from transcriptome sequencing, which is an important method for gene expression studies because of its high throughput, accuracy and efficiency [44,45]. Therefore, through analyzing the transcriptome data, the reference genes stably expressed in different cultivars, tissues, and under various experimental conditions can be screened out. This method has been successfully applied to several studies in plant species such as Prunus mume [2], seashore paspalum [40], and chrysanthemum [6,38,39]. In the present study, based on transcriptome data, eleven candidate reference genes with relatively stable expression (|log2Ratio| < 0.3) were screened out. Another four reference genes, including PSAA, F-box, TUA (α-tublin), and TUB (β-tublin), were eliminated because of their wide variation in expression. This effective approach to identify candidate reference genes based on transcriptome libraries provided a solid foundation for further evaluating the expression stability of these candidates.

The superiority of a comprehensive evaluation of stability of reference genes using four programs

Previous studies have revealed that there is no consensus on which kind of statistical program should be used to analyze the stability of reference genes because different programs based on distinct algorithms generate potentially conflicting results [9]. In the research of tobacco, two programs (geNorm and NormFinder) were used to evaluate the stability of reference genes, however, the ranking orders were inconsistent in stress-treated sample set [46]. This suggests that it is necessary to use at least three different algorithms to obtain reliable results. Multi-algorithm analysis had been performed to select suitable reference genes under different experimental situations in plants such as Chrysanthemum [9], Rhododendron molle [47], tree peony [48], seashore paspalum [40], and Lagerstroemia [41]. Therefore, in the present study, four computational programs (geNorm, NormFinder, BestKeeper, and RefFinder) were applied to evaluating stability of candidate reference genes. The results of geNorm and NormFinder were similar, but they exhibited quite differences from those of BestKeeper. In set 1, geNorm and NormFinder determined SAND and GAPDH as most stable genes and SKIP16 as the least stable gene, while SKIP16 was the most stable gene in BestKeeper analysis. But when we used RefFinder to integrate the ranking orders generated by these three algorithms, SAND and GAPDH were still the most stable reference genes. These results coincided with studies on flower development in Chrysanthemum [9] and Lagerstroemia [41]. Hence, a comprehensive analysis using four programs can generate more reliable reference genes.

The stability of candidate reference genes used in shoot development and gravitropic response in chrysanthemum

By interpreting results from four frequently-used methods (geNorm, NormFinder, BestKeeper, and RefFinder), several stable reference genes were identified in this study, including SAND, PP2A-1, and TIP41. The SAND family proteins play major roles in the downstream regulatory pathway of endocytic transport [49]. The current study determined SAND as the best candidate for normalization in different chrysanthemum genotypes during different shoot developmental stages. SAND have also been reported as an optimal reference gene in previous researches, such as in Petunia hybrida during leaf and flower development [50], in different citrus organs and following different biotic stresses [51], in vegetative tissues and organs during berry development [52], in Chrysanthemum during flower development [9], and in salt-treated seashore paspalum [40]. PP2A-1, which encodes the serine/threonine phosphatase, plays a prominent role in metabolism, DNA replication, transcription and translation, cell cycle, and signal transduction [53]. It was used as a superior reference gene in different cotton plant organs [54], different color lines during flower developmental stages of cineraria [55], and cold-treated seashore paspalum [40]. TIP41 (TAP42 INTERACTING PROTEIN OF 41 kDa), which modifies cell growth in response to nutrient status and environmental conditions [56], has been revealed as an appropriate reference gene in whole tomato developmental process [57]. In our study, PP2A-1 and TIP41 were stably expressed in set 4 and set 5, thus determining them as superior reference genes in gravitropic response of YT and FH. The traditional reference gene ACT was used for data normalization in qPCR analysis during graviresponse in Arabidopsis seedlings [58], peanut gynophores [20] and chrysanthemum shoot [23]. However, in the present study, the expression of ACT displayed unacceptable variation during shoot graviresponse. The variable expression of ACT may closely related to the involvement of cytoskeleton in the gravitropic bending growth [23]. ACT may regulate the asymmetric elongation of the upper and lower half of the bending part of the stem, so it inappropriate to use ACT as reference gene during shoot gravity response.

The reliability of comprehensive evaluation on reference genes

Expression patterns of target genes were found to vary significantly when normalized by stable and unstable reference genes, which led to conflicting results [9]. In this study, to further verify the reliability of the selected reference genes, the expression levels of target genes (LA1 and PIN1) were analyzed using two most stable and one least stable reference genes based on integrating results of four programs. In set 1, set 2, and set 3, when normalizing with the two most stable genes, the expression patterns of LA1 were similar, which was highly consistent with the results that the relative expression level of ZmLA1 in la1-ref mutant plants was obviously lower than that in wild maize [31]. However, the results were quite different when the least stable gene was used for normalization. During the process of stem graviresponse in YT or FH, transcription levels of PIN1 were almost the same when normalized with two most stable reference genes, while the expression pattern was contradictory to the above outcome when normalized using the least stable gene. The significant discrepancies in expression level of PIN1 might be responsible for the differences in gravitropic response between YT and FH. The weaker gravitropic response in YT might be due to the reduced asymmetric distribution of IAA which is caused by the absence of significant discrepancies in expression level of PIN1 between the upper and lower half of stem. These results validated the reliability of screened reference genes and also indicated the importance of selecting suitable internal control genes for qRT-PCR analysis.

Conclusions

This research provides the first systematic evaluation of reference genes for qRT-PCR analysis in shoot development and gravitropic response of prostrate and erect chrysanthemum. In general, TIP41 was the most stable gene for all the samples. SAND could be applied as a superior reference gene in different genotypes and during the process of stem development. The suitable reference genes for gravitropic response could be PP2A-1. The expression patterns of LA1 and PIN1 further verified the importance of selection of suitable reference genes for qRT-PCR analysis. Specific conclusions drawn from this study could provide more accurate and reliable qRT-PCR normalization for future studies on the expression patterns of genes regulating plant architecture of chrysanthemums.

Melting curve and melting peak of candidate reference genes.

(TIF)

Click here for additional data file.

Standard curve of candidate reference genes in the five-fold serial dilution.

(TIF)

Click here for additional data file.

Melting curve and standard curve of target genes.

(TIF)

Click here for additional data file.

FPKM and log2Ratio (III_P vs III_E) of candidate reference genes.

(DOCX)

Click here for additional data file.

Primer sequences for gene cloning of reference genes.

(DOCX)

Click here for additional data file.