1. Verification of primer specificity and PCR efficiency
In the present study, 11 traditional housekeeping genes were chosen as candidate reference genes used in qRT-PCR studies namely, 18S ribosome RNA(18S rRNA), Actin(ACT), Cyclophilin(CYP), Elongation factor 1-α(EF1a), Elongation factor 1-β(EF1b), Glyceraldehyde-3-phosphate dehydrogenase(GAPDH), Malate dehydrogenease (MDH), Ubiquitin-conjugating enzyme(UBC), Tubulin-α(TUBa), Tubulin -β(TUBb), and Ubiquitin-A(UBQ) in various conditions(Supplemental Table 1). The primers sequences of 11 candidate reference were describe in Table 1. The melting curve analysis following qRT-PCR presented a single peak, indicated that the specificity of primers (Supplementary Fig. 1). The amplification efficiency of primers were assessed by qRT-PCR after 45 cycles from all sample under different gradient concentrations (Supplementary Fig. 2), ranged from 91.9-102.5%, and the correlation coefficient (R2) near 0.99 (Table 1).
Table 1
Primers sequences of the candidate reference genes
Gene name
|
Genes symbol
|
Gene ID
|
Primers sequences(5′-3′)
|
Primer Efficiency
|
R2
|
18S ribosome RNA
|
18S
|
MN650750
|
F: CTGAGAAACGGCTACCACATC
|
102.51
|
0.963
|
R: CGAAGAGCCCGGTATTGTTATT
|
ACT-β
|
ACT
|
MN627205
|
F: GAGGAGAGGATCTGTCGTAAA
|
91.96
|
0.998
|
R: GATAACAAGGAGAGGCCAAAG
|
Cyclophilin
|
CYP
|
MN627206
|
F: AACAAGTTCGCCGATGAG
|
95.05
|
0.999
|
R: GTCTTCGCAGTGCAAATAAAG
|
Elongation factor 1-α
|
EF1a
|
MN627207
|
F: CCCACTTCAGGGTGTTTAC
|
97.99
|
0.994
|
R: CGAAGGTGACAACCATACC
|
Elongation factor 1-β
|
EF1b
|
MN627208
|
F: GGCTGCTAAAGCATCTACAA
|
94.41
|
0.998
|
R: CGAACAGCTTCCTCTAGTTTC
|
Glyceraldehyde 3-phosphate dehydrogenase
|
GADPH
|
MN627209
|
F: GTTGGTGACTGTAGGTCAAG
|
93.73
|
0.999
|
R: AGGTCCAACACTCGGTTA
|
Malate dehydrogenease
|
MADH
|
MN627210
|
F: GCTGGTCTCATCTATTCTTTCC
|
98.61
|
0.999
|
R: CGTCCAACTTCTTCCTTGAG
|
Tubulin-α
|
TUBa
|
MN627211
|
F: GGCTTGTGTCTCAGGTTATT
|
94.9
|
0.999
|
R: GTGGATATGGGACCAAGTTAG
|
Tubulin-β
|
TUBb
|
MN627212
|
F: ATATGAGGATGAGGAGGAAGG
|
98.35
|
0.999
|
R: CCCATAATCAGCCACTGTAAA
|
Ubiquitin
|
UBQ
|
MN627213
|
F: GTGGATGTTGATGGATGAAAC
|
98.68
|
0.999
|
R: GTACTTACAGAGCGTCCTTAC
|
Ubiquitin-conjugating enzyme
|
UBC
|
MN627214
|
F: CCACCCAAGGTAGCATTTAG
|
95.95
|
0.982
|
R: CTGGGCTCCATTGTTCTTTA
|
2. Expression profile of the candidate reference genes
Cycle threshold (Ct) value of qRT-PCR is the number of cycles when the fluorescence signal intensity reach a set threshold level of detection, and it represents the transcript levels of genes in test samples. The transcript level of 11 reference genes exhibited a relatively high variation from all the experimental set as show in Supplemental Table 2. The median values of Ct ranged from 16.73 for 18S rRNA to 26.41 for GAPDH, most values lying between 20 and 23 across all samples. 18S rRNA showed the most high expression level with the lowest mean Ct value (17.66±2.62),followed by CYP(21.05±2.39), MDH(21.97±2.26), ACT(22.82±2.20), UBQ(23.05±3.08), UBC(23.25±2.39), EF1a(23.79±2.98), EF1b(23.80±2.41), TUBb(25.40±2.21), GAPDH(27.69±2.86). The expression of TUBa(27.94±3.22) candidate gene was found to be the lowest level from all the test sample. The genes with higher SD of Ct values indicated more variable expression compared to these with lower SD. ACT showed the smallest variation in gene expression with lowest SD(22.82±2.20), while TUBa(27.94±3.22) appearance the most variable levels of expression(Figure 1).
3. Expression stability analysis of candidate reference genes
In order to select the best reference genes, RefFinder online tool was used to rank the expression stability values. This method include comparative ΔCt, geNorm, NormFinder and BestKeeper, four popular algorithms, determinates a comprehensive expression stability ranking combined with above four data generated from different algorithms(Fig. 2; Table 2).
Table 2
|Expression stability ranking of the 11 candidate reference genes.
Method/Rank
|
1
|
2
|
3
|
4
|
5
|
6
|
7
|
8
|
9
|
10
|
11
|
(A)Ranking order under cold stress(better-good-average)
|
Delta CT
|
ACT
|
GAPDH
|
CYP
|
EF1b
|
MDH
|
TUBb
|
UBC
|
EF1a
|
18S
|
UBQ
|
TUBa
|
BestKeeper
|
UBC
|
TUBb
|
EF1b
|
ACT
|
CYP
|
MDH
|
GAPDH
|
TUBa
|
EF1a
|
UBQ
|
18S
|
NormFinder
|
ACT
|
GAPDH
|
CYP
|
EF1b
|
TUBb
|
MDH
|
UBC
|
EF1a
|
18S
|
TUBa
|
UBQ
|
geNorm
|
ACT / EF1b
|
GAPDH
|
CYP
|
MDH
|
TUBb
|
UBC
|
18S
|
EF1a
|
UBQ
|
TUBa
|
Comprehensive Ranking
|
ACT
|
EF1b
|
GAPDH
|
CYP
|
UBC
|
TUBb
|
MDH
|
EF1a
|
18S
|
TUBa
|
UBQ
|
(B)Ranking order under heat stress(better-good-average)
|
Delta CT
|
EF1b
|
ACT
|
18S
|
MDH
|
EF1a
|
CYP
|
TUBb
|
GAPDH
|
UBQ
|
TUBa
|
UBC
|
BestKeeper
|
TUBb
|
TUBa
|
MDH
|
ACT
|
18S
|
EF1b
|
GAPDH
|
CYP
|
EF1a
|
UBC
|
UBQ
|
NormFinder
|
EF1b
|
ACT
|
EF1a
|
18S
|
MDH
|
TUBb
|
GAPDH
|
CYP
|
UBQ
|
TUBa
|
UBC
|
geNorm
|
ACT / EF1b
|
MDH
|
18S
|
TUBb
|
CYP
|
EF1a
|
GAPDH
|
UBQ
|
TUBa
|
UBC
|
Comprehensive Ranking
|
EF1b
|
ACT
|
MDH
|
TUBb
|
18S
|
EF1a
|
TUBa
|
CYP
|
GAPDH
|
UBQ
|
UBC
|
(C)Ranking order under NaCl stress(better-good-average)
|
Delta CT
|
ACT
|
CYP
|
EF1b
|
GAPDH
|
UBC
|
MDH
|
EF1a
|
UBQ
|
18S
|
TUBb
|
TUBa
|
BestKeeper
|
18S
|
MDH
|
CYP
|
EF1b
|
UBC
|
ACT
|
GAPDH
|
EF1a
|
UBQ
|
TUBb
|
TUBa
|
NormFinder
|
CYP
|
ACT
|
EF1b
|
GAPDH
|
EF1a
|
UBC
|
UBQ
|
MDH
|
18S
|
TUBb
|
TUBa
|
geNorm
|
ACT / EF1b
|
UBC
|
MDH
|
CYP
|
GAPDH
|
UBQ
|
EF1a
|
18S
|
TUBb
|
TUBa
|
Comprehensive Ranking
|
ACT
|
CYP
|
EF1b
|
MDH
|
UBC
|
GAPDH
|
18S
|
EF1a
|
UBQ
|
TUBb
|
TUBa
|
(D)Ranking order under PEG stress (better-good-average)
|
Delta CT
|
ACT
|
UBC
|
EF1b
|
MDH
|
CYP
|
EF1a
|
UBQ
|
18S
|
GAPDH
|
TUBa
|
TUBb
|
BestKeeper
|
ACT
|
UBC
|
EF1b
|
CYP
|
MDH
|
TUBa
|
TUBb
|
18S
|
EF1a
|
GAPDH
|
UBQ
|
NormFinder
|
UBC
|
ACT
|
EF1b
|
CYP
|
MDH
|
EF1a
|
UBQ
|
18S
|
GAPDH
|
TUBa
|
TUBb
|
geNorm
|
ACT / EF1b
|
UBC
|
CYP
|
MDH
|
EF1a
|
UBQ
|
GAPDH
|
18S
|
TUBa
|
TUBb
|
Comprehensive Ranking
|
ACT
|
UBC
|
EF1b
|
CYP
|
MDH
|
EF1a
|
UBQ
|
18S
|
TUBa
|
GAPDH
|
TUBb
|
(E)Ranking order under ABA treatment (better-good-average)
|
Delta CT
|
ACT
|
CYP
|
TUBb
|
EF1b
|
UBC
|
18S
|
MDH
|
EF1a
|
UBQ
|
GAPDH
|
TUBa
|
BestKeeper
|
MDH
|
UBC
|
ACT
|
EF1b
|
TUBb
|
CYP
|
18S
|
EF1a
|
GAPDH
|
UBQ
|
TUBa
|
NormFinder
|
ACT
|
CYP
|
EF1b
|
TUBb
|
18S
|
UBC
|
EF1a
|
UBQ
|
MDH
|
GAPDH
|
TUBa
|
geNorm
|
MDH / UBC
|
ACT
|
TUBb
|
CYP
|
EF1b
|
18S
|
EF1a
|
UBQ
|
GAPDH
|
TUBa
|
Comprehensive Ranking
|
ACT
|
UBC
|
MDH
|
CYP
|
TUBb
|
EF1b
|
18S
|
EF1a
|
UBQ
|
GAPDH
|
TUBa
|
(F)Ranking order under GA treatment (better-good-average)
|
Delta CT
|
UBQ
|
EF1b
|
MDH
|
EF1a
|
UBC
|
ACT
|
GAPDH
|
TUBb
|
CYP
|
TUBa
|
18S
|
BestKeeper
|
EF1b
|
MDH
|
UBC
|
ACT
|
UBQ
|
EF1a
|
GAPDH
|
CYP
|
TUBb
|
TUBa
|
18S
|
NormFinder
|
UBQ
|
EF1a
|
EF1b
|
GAPDH
|
MDH
|
ACT
|
UBC
|
TUBb
|
CYP
|
TUBa
|
18S
|
geNorm
|
MDH / UBC
|
ACT
|
EF1b
|
UBQ
|
EF1a
|
GAPDH
|
TUBb
|
CYP
|
TUBa
|
18S
|
Comprehensive Ranking
|
EF1b
|
UBQ
|
MDH
|
UBC
|
EF1a
|
ACT
|
GAPDH
|
TUBb
|
CYP
|
TUBa
|
18S
|
(G)Ranking order under MeJA treatment (better-good-average)
|
Delta CT
|
ACT
|
CYP
|
TUBb
|
UBQ
|
MDH
|
GAPDH
|
18S
|
UBC
|
EF1b
|
EF1a
|
TUBa
|
BestKeeper
|
UBC
|
MDH
|
ACT
|
TUBb
|
CYP
|
EF1b
|
18S
|
GAPDH
|
UBQ
|
TUBa
|
EF1a
|
NormFinder
|
ACT
|
CYP
|
TUBb
|
UBQ
|
GAPDH
|
MDH
|
18S
|
EF1b
|
UBC
|
EF1a
|
TUBa
|
geNorm
|
ACT / MDH
|
UBC
|
TUBb
|
CYP
|
EF1b
|
UBQ
|
18S
|
GAPDH
|
EF1a
|
TUBa
|
Comprehensive Ranking
|
ACT
|
MDH
|
CYP
|
TUBb
|
UBC
|
UBQ
|
GAPDH
|
EF1b
|
18S
|
EF1a
|
TUBa
|
(H)Ranking order under ETHYLENE treatment (better-good-average)
|
Delta CT
|
UBC
|
CYP
|
ACT
|
UBQ
|
MDH
|
TUBb
|
18S
|
EF1b
|
GAPDH
|
TUBa
|
EF1a
|
BestKeeper
|
MDH
|
UBC
|
ACT
|
EF1b
|
CYP
|
TUBb
|
GAPDH
|
UBQ
|
18S
|
TUBa
|
EF1a
|
NormFinder
|
CYP
|
UBC
|
UBQ
|
ACT
|
MDH
|
TUBb
|
18S
|
GAPDH
|
EF1b
|
TUBa
|
EF1a
|
geNorm
|
ACT / UBC
|
MDH
|
CYP
|
TUBb
|
EF1b
|
UBQ
|
GAPDH
|
18S
|
TUBa
|
EF1a
|
Comprehensive Ranking
|
UBC
|
ACT
|
CYP
|
MDH
|
UBQ
|
TUBb
|
EF1b
|
18S
|
GAPDH
|
TUBa
|
EF1a
|
(I)Ranking Order of root (Better–Good–Average)
|
Delta CT
|
CYP
|
EF1a
|
GAPDH
|
TUBa
|
EF1b
|
ACT
|
UBC
|
UBQ
|
MDH
|
TUBb
|
18S
|
BestKeeper
|
18S
|
TUBa
|
CYP
|
EF1a
|
ACT
|
UBQ
|
GAPDH
|
EF1b
|
UBC
|
TUBb
|
MDH
|
NormFinder
|
EF1a
|
CYP
|
GAPDH
|
TUBa
|
EF1b
|
ACT
|
UBC
|
UBQ
|
MDH
|
TUBb
|
18S
|
Genorm
|
CYP / TUBa
|
EF1a
|
GAPDH
|
ACT
|
EF1b
|
UBC
|
UBQ
|
MDH
|
TUBb
|
18S
|
Comprehensive ranking
|
CYP
|
EF1a
|
TUBa
|
GAPDH
|
ACT
|
EF1b
|
18S
|
UBQ
|
UBC
|
MDH
|
TUBb
|
(J)Ranking Order of stem (Better–Good–Average)
|
Delta CT
|
ACT
|
MDH
|
UBC
|
TUBb
|
EF1b
|
18S
|
GAPDH
|
CYP
|
UBQ
|
EF1a
|
TUBa
|
BestKeeper
|
UBC
|
MDH
|
TUBb
|
ACT
|
EF1b
|
CYP
|
18S
|
GAPDH
|
UBQ
|
EF1a
|
TUBa
|
NormFinder
|
ACT
|
TUBb
|
18S
|
MDH
|
EF1b
|
UBC
|
GAPDH
|
UBQ
|
CYP
|
EF1a
|
TUBa
|
Genorm
|
ACT / UBC
|
MDH
|
TUBb
|
EF1b
|
CYP
|
18S
|
GAPDH
|
UBQ
|
EF1a
|
TUBa
|
Comprehensive ranking
|
ACT
|
UBC
|
MDH
|
TUBb
|
EF1b
|
18S
|
CYP
|
GAPDH
|
UBQ
|
EF1a
|
TUBa
|
(K)Ranking Order of leaf (Better–Good–Average)
|
Delta CT
|
CYP
|
EF1b
|
ACT
|
UBQ
|
EF1a
|
MDH
|
GAPDH
|
TUBb
|
18S
|
TUBa
|
UBC
|
BestKeeper
|
ACT
|
MDH
|
CYP
|
TUBb
|
EF1b
|
UBC
|
18S
|
EF1a
|
UBQ
|
GAPDH
|
TUBa
|
NormFinder
|
CYP
|
EF1b
|
ACT
|
EF1a
|
UBQ
|
GAPDH
|
MDH
|
TUBb
|
18S
|
TUBa
|
UBC
|
Genorm
|
ACT / MDH
|
EF1b
|
CYP
|
TUBb
|
UBQ
|
EF1a
|
GAPDH
|
18S
|
TUBa
|
UBC
|
Comprehensive ranking
|
ACT
|
CYP
|
EF1b
|
MDH
|
UBQ
|
EF1a
|
TUBb
|
GAPDH
|
18S
|
UBC
|
TUBa
|
(L)Ranking Order of all sample (Better–Good–Average)
|
Delta CT
|
ACT
|
EF1b
|
CYP
|
MDH
|
EF1a
|
GAPDH
|
UBQ
|
TUBb
|
18S
|
UBC
|
TUBa
|
BestKeeper
|
ACT
|
MDH
|
UBC
|
CYP
|
EF1b
|
TUBb
|
18S
|
UBQ
|
EF1a
|
GAPDH
|
TUBa
|
NormFinder
|
ACT
|
CYP
|
EF1b
|
MDH
|
EF1a
|
GAPDH
|
UBQ
|
TUBb
|
18S
|
UBC
|
TUBa
|
Genorm
|
ACT / MDH
|
EF1b
|
CYP
|
EF1a
|
UBQ
|
GAPDH
|
18S
|
TUBb
|
UBC
|
TUBa
|
Comprehensive ranking
|
ACT
|
MDH
|
EF1b
|
CYP
|
EF1a
|
UBQ
|
GAPDH
|
UBC
|
TUBb
|
18S
|
TUBa
|
ΔCt algorithm analysis
The comparative ΔCt algorithm ranks the potential candidate reference genes expression stabilities based on the standard deviation values(Supplemental Table 3) [17]. As determined by ΔCt, ACT was consistently identified as the most stable gene through the total samples(Supplemental Fig. 3). CYP was the most stable reference gene in macadamia root and leaf, while ACT was the most stable reference gene in the stem in the different organs dataset from macadamia. In the subset of cold stress, NaCl stress, PEG stress, ABA treatment and MeJA treatment, ACT was still considered as the best stable gene in the 11 candidate reference genes, and TUBa was the least stable gene under the cold stress, NaCl stress, PEG stress, ABA treatment and MeJA treatment. Under the heat stress EF1b ranked top than ACT gene, and was determined the best stable reference gene. UBQ was ranked as the stable reference gene, while 18S was the least stable reference gene under the GA treatment. In the ethylene treatment, UBC was the most highly ranked, and EF1a performanced the least stable gene(Fig. 2; Table 2, Supplemental Table 3).
geNorm Analysis
The geNorm program used the M-value to assess stability of the reference gene expression. The M-value was calculated at each step during stepwise exclusion of the least stable reference gene until two best gens were obtained[18]. The top two stable candidate reference gene were ACT / MDH for all samples(Supplemental Fig. 4), leaf and MeJA treatment, ACT and EF1b under cold stress, heat stress, NaCl stress and PEG stress, MDH / UBC in the GA treatment and ABA treatment, ACT and UBC for stem and ethylene treatment show the most stability. Additionally, as show in Fig. 2, Table 2 and Supplemental Table 4, the least stably expressed was TUBa under cold stress, NaCl stress, ABA treatment, MeJA treatment, stem and all samples, UBC under heat stress and leaf, TUBb under PEG stress, 18S under the GA treatment and root, EF1a under the ethylene treatment.
NormFinder Analysis
The NormFinder algorithm ranking candidate reference gene stability based on a variance estimation approach[19]. The NormFinder ranking orders are similar to the geNorm analysis. According this analysis(Figure S4; Table 2; Supplemental Table 5), the most stable gene were ACT for cold stress, ABA treatment, MeJA treatment, stem and all samples, EF1b for heat stress, CYP for leaf, ethylene treatment and NaCl stress, UBC for PEG stress, UBQ for GA treatment, EF1a for root. (Supplemental Fig. 5)
BestKeeper Analysis
The BestKeeper algorithm determined the stabilities of candidate reference genes based on the CV ± SD values[20]. The stability rankings by BestKeeper was quite different from the ranking order determined by ΔCt, geNorm and NormFinder algorithms(Supplemental Fig. 6). ACT was considered as the best stable gene under all samples, PEG stress, and leaf. TUBb in heat stress, UBC for cold stress, MeJA treatment and stem subset, 18S for NaCl stress and root, MDH for ABA treatment and ethylene treatment, EF1b for GA treatment were all ranked in top stable genes in different situation. Additionally, 18S and TUBa were determined as the least stable reference gene in the most samples by BestKeeper algorithm(Fig. 2; Table 2; Supplemental Table 6).
Comprehensive Ranking
In order to obtain a consensus result of the best stable reference genes as recommended by the four above approached according to the RefFinder[21], the geometric mean of four algorithms corresponding rankings for each candidate gene were calculated(Table 2; Fig. 2; Supplemental Tables 7 and 8). ACT was ranked the first stable gene in all samples, cold stress, NaCl stress, PEG stress, ABA treatment, MeJA treatment, stem tissue and leaf, EF1b for heat stress and GA treatment, CYP for root, UBC for ethylene treatment (Fig. 2; Table 2).
4. Validation of the reference genes with the SAD gene
To further validate the selected reference genes, the relative expression level of a target gene SAD in macadamia under different experimental conditions were evaluated using qRT-PCR. It was normalized using the most stable reference genes and two least stable reference genes as the internal control both singly and combination under different treatment subsets (Fig. 3). The relative expression abundance of SAD showed substantial divergence when normalized to different kinds of reference genes. In this study, the water treatment was designed as a negative control, when the most stable genes ACT and CYP were employed respectively and together to normalized expression levels of SAD gene, the expression abundant was not significantly effect; when the least stable genes TUBa were used to normalized expression level of target gene, the expression abundant was up-regulation dramatically. When the EF1a was used as the reference, the expression pattern of SAD gene was similar with the most stable genes normalization results. The expression level of SAD gene gradually declined under the cold stress using ACT and EF1b as the internal control, while the expression level was increased dramatically at 6h and 24h when the least stable reference genes TUBa was used, and increased sharply at 24h for UBQ normalized as internal reference. Under heat stress, the expression pattern of SAD gene were similar when the most stable reference genes (ACT and EF1b) and the least stable reference genes( UBC and UBQ) were used as controls, but the highest expression level at 12h was 2.5-fold higher by using UBQ normalizations than using the most stable reference gene ACT normalization. In PEG stress, the expression level of SAD gene was up-regulated at 24h and then down-regulated at 48h when using ACT and UBC genes normalization, while the expression appeared overestimated when normalized using GAPDH and underestimated for TUBb. Similarly, in NaCl stress, the expression profiles of the SAD showed same trends when stable reference genes ACT and CYP were used as internal control, when the least stable reference gene TUBb and TUBa were used as reference control, the expression level of target gene was significantly overestimated and underestimated respectively. Under ET and ABA treatment, the relative expression trends of SAD gene showed a similarly trend, appeared slightly descent at 6h and then increased at 24h when the most stable genes were used in combination and independent for normalization, while appeared down-regulated at 12h and 24h when normalized by using the least stable reference genes. In MeJA stress, a similar expression pattern was observed when either ACT or MDH was used alone and combination for normalization. While the expression levels of SAD gene were overestimated when the least stable genes TUBa and EF1a were used for normalization. In GA stress, the relative expression of SAD gene increased at 6h when using EF1b and 18S as internal control, while declined slightly at 6h when normalized by using TUBa and UBQ as reference genes. Above all, the result showed that ACT and EF1b were more suitable for different treatment subset and TUBa failed standardized the expression data effectively.
In addition, to validate the selected reference genes, 4 genes involved fatty acid accumulation including FATA, FAD, Oleosin and SAD in qRT-PCR were compared with the expression pattern in RNA-seq results. In the qRT-PCR results, Actin was used as the normalizer to quantify the expression. As show in the Fig. 4, all the selected genes revealed similar expression file in the both methods. Therefore, our result provide the reliable reference gene in Macadamia.