Selection and validation of reference genes for quantitative expression analysis of miRNAs and mRNAs in Poplar

Background Quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) is a rapid and sensitive approach to identify miRNA and protein-coding gene expression in plants. However, because of the specially designated reverse transcription and shorter PCR products, very few reference genes have been identified for the quantitative analysis of miRNA expression in plants, and different internal reference genes are needed to normalize the expression of miRNAs and mRNA genes respectively. Therefore, it is particularly important to select the suitable common reference genes for normalization of quantitative PCR of miRNA and mRNA. Results In this study, a modified reverse transcription PCR protocol was adopted for selecting and validating universal internal reference genes of mRNAs and miRNAs. Eight commonly used reference genes, four stably expressed novel genes in Populus tremula, three small noncoding RNAs and three conserved miRNAs were selected as candidate genes, and the stability of their expression was examined across a set of 38 tissue samples from four developmental stages of poplar clone 84K (Populus alba × Populus glandulosa). The expression stability of these candidate genes was evaluated systematically by four algorithms: geNorm, NormFinder, Bestkeeper and DeltaCt. The results showed that Eukaryotic initiation factor 4A III (EIF4A) and U6-2 were suitable for samples of the callus stage; U6-1 and U6-2 were best for the seedling stage; Protein phosphatase 2A-2 (PP2A-2) and U6-1 were best for the plant stage; and Protein phosphatase 2A-2 (PP2A-2) and Oligouridylate binding protein 1B (UBP) were the best reference genes in the adventitious root (AR) regeneration stage. Conclusions The purpose of this study was to identify the most appropriate reference genes for qRT-PCR of miRNAs and mRNAs in different tissues at several developmental stages in poplar. U6-1, EIF4A and PP2A-2 were the three most appropriate reference genes for qRT-PCR normalization of miRNAs and mRNAs during the plant regeneration process, and PP2A-2 and UBP represent the best reference genes in the AR regeneration stage of poplar. This work will benefit future studies of expression and function analysis of miRNAs and their target genes in poplar. Electronic supplementary material The online version of this article (10.1186/s13007-019-0420-1) contains supplementary material, which is available to authorized users.


Introduction
Given its high sensitivity, quantitative accuracy, low cost and specificity, quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) has become the most common and widely used technique for quantifying miRNA expression and mRNA transcript levels among different tissues and experimental conditions in plants [1,2]. However, the accuracy of qRT-PCR is easily affected by several factors, including the quality of RNA samples, reverse transcription efficiency, cDNA quality and amount, and differences in extraneous tissue and cell activities [2][3][4]. To avoid bias in qRT-PCR analysis, validation of suitable reference genes for data normalization is an elementary prerequisite for each

Open Access
Plant Methods *Correspondence: lumz@caf.ac.cn † Fang Tang and Liwei Chu contributed equally to this work 1 State Key Laboratory of Tree Genetics and Breeding, Key Laboratory of Tree Breeding and Cultivation of the National Forestry and Grassland Administration, Research Institute of Forestry, Chinese Academy of Forestry, Beijing 100091, China Full list of author information is available at the end of the article experimental condition in different tissues or species [5]. However, no single reference gene can be universal under all experimental situations, even including the most stable reference gene(s) reported [6,7]. Therefore, optimal reference genes should be validated for different species, tissues or specific treatments.
MicroRNAs (miRNAs) are endogenous ~ 22 nt small noncoding RNAs that guide the cleavage or repress the translation of their target mRNAs by approximate basepairing rules [8,9] or mediate mRNA decay by directing rapid deadenylation of mRNAs [10]. In plants, miR-NAs are master regulators in controlling developmental processes and in response to biotic and abiotic stress responses [11][12][13][14][15]. Due to its short sequence (only ~ 22 nt in length), the quantification of miRNAs by qRT-PCR requires extending the length of mature miRNAs using stem-loop primers [7,16] or adding poly(A)-tails [17][18][19][20]. This extension requirement causes different internal reference genes to be commonly used for normalization in qRT-PCR of miRNAs and mRNAs. Some housekeeping genes, such as actin 7 (ACT 7), eukaryotic initiation factor 4A III (EIF4A), polyubiquitin (UBQ), glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and protein phosphatase 2A-2 (PP2A-2), were widely adopted for gene expression analysis in diverse plants as reference genes [21][22][23]. Several noncoding RNA and small RNA, such as 5.8S ribosomal RNA (5.8S rRNA) and U6 small nuclear RNA (U6 snRNA), are commonly used as reference genes for miRNA quantity [24,25]. Hurteau developed a modified universal reverse transcription PCR protocol, in which mature miRNAs could be polyadenylated by poly (A) polymerase and reverse transcribed into cDNA using oligo-dT primers [17]. Then, mRNAs and miRNAs could be specifically amplified and quantified at same transcriptional level, and the relative quantification of a miRNA and its predicted mRNA target can be both assessed precisely [17]. In this case, it is particularly important to select a suitable reference gene for normalization in quantitative PCR of miRNA and mRNA. As a typical model woody plant, Populus has many advantages in basic research, such as rapid and perennial growth, moderate genome size, biomass-related traits and facile transformation [26]. Completion of the genomic sequence for Populus trichocarpa (black cottonwood) [27] has led to the development of genomic and molecular resources, and the ideal genetic transformation system provides a powerful genetic analysis tool for dissecting adaptive traits in poplar [28,29]. Poplar clone 84K (Populus alba × Populus glandulosa) is now commonly used for gene functional studies because it is easier to obtain transgenic plants through Agrobacterium tumefaciens-mediated leaf discs [30][31][32]. The regeneration of transgenic plants involves callus induction, shoot differentiation, seedling culture and plant growth. This process is time consuming, requiring 2-3 months for tissue culture seedling and approximately 3-4 months for plant growth in soil. Therefore, we perform transgenic identification and gene expression analysis of early regenerated shoots and/or roots even small seedlings to reduce the time for identifying gene function in transgenic plants. In addition, the expression levels of miRNAs or target genes in different tissues are also required for analysis between transgenic and normal plants. However, most miRNA and mRNA expression levels vary greatly in different developmental stages and tissues, so it is necessary to identify a more stably expressed miRNA or gene as the internal reference to normalize expression using qRT-PCR.
In this study, we have tested 18 genes and noncoding RNAs for candidate reference genes. The expression stability of these genes was validated across a set of 38 tissue samples from four developmental biological processes of 84K poplar using a modified universal reverse transcription PCR protocol [17]. The cycle threshold (Ct) values of candidates were used to evaluate the expression stability using four algorithms: geNorm, NormFinder, Bestkeeper and DeltaCt. U6-1, EIF4A and PP2A-2 were the top three most appropriate reference genes for qRT-PCR of miR-NAs and mRNAs (including miRNA target genes) during the plant regeneration process, and PP2A-2 and UBP were the best combination as reference genes in the AR regeneration stage of poplar.

Verification of amplification and efficiency of the primers
A total of 12 protein-coding genes and 6 small noncoding RNAs were used as candidate reference genes for quantitative detection of miRNAs and mRNAs. The qRT-PCR primer sequences and amplicon characteristics of these candidate genes in 84K poplar are presented in Table 1. The PCR amplification specificities were confirmed by melting curves (Additional file 1: Fig. S1), agarose gel electrophoresis (Additional file 1: Fig. S2) and sequencing (Additional file 2: Fig. S3), which demonstrated the specific product of expected size and sequence. The qRT-PCR products ranged from 49 to 147 bp, and the sequence similarity between 84K poplar and Populus trichocarpa ranged from 97 to 100% despite belonging to different species. Therefore, the primers of these reference genes could also be used in other poplar species. To evaluate the amplification efficiency of pair-primers, the standard curves were obtained using a set of 10-fold diluted cDNA templates. The amplification efficiency (E) of the 18 candidate reference genes ranged from 96.16 to 116.69% and the regression coefficient (R 2 ) varied between 0.941 and 1.000 (Table 1). These results suggest that the primers of all the candidate reference genes exhibit high amplification efficiency and specificity in the qRT-PCR system.

Expression profile of candidate reference genes during plant regeneration
To evaluate the stability of the reference genes across all experimental samples, the transcript abundance of the 18 candidate reference genes was assessed based on mean Ct values. The average Ct values for the 18 candidate reference genes ranged from 21 to 33, with most values between 26 and 27 across all samples. miR482 had the highest expression level with a Ct value of 20.65 cycle, whereas UBP was the lowest abundantly expressed gene with Ct values up to 31.47 (Fig. 1). The Ct values of EIF4A (26.50 ± 0.55) and U6-2 (27.27 ± 0.59) with minimum SD data indicate that these genes are the most stable genes in all the samples. The next most stable genes include  Table S1).

Expression stability of candidate reference genes during plant regeneration
To detect the expression stability of 18 candidate reference genes more accurately, four software programs, including geNorm, NormFinder, BestKeeper and Delta-Ct method, were used for statistical analysis. These candidate reference genes were evaluated by each program and ranked from the most to the least stably expressed genes. Then, the geometric mean of each gene was calculated and reordered using RankAggreg software. Data of each reference gene from samples of different developmental stages were analyzed separately and then integrated together.
The M-values of 18 candidate reference genes calculated by geNorm were all less than 1.5 in different tissues at plants developmental stages. U6-1 and U6-2 exhibited the highest stable expression with least M-values of 0.064 (Fig. 2a) and 0.082 (Fig. 2b) at the callus and seedling stage, respectively. PP2A-2 and U6-1 both exhibited the highest stable expression with the lowest M-value of 0.1456 at the plant stage (Fig. 2c), and EIF4A and U6-1 exhibited the lowest M-value at 0.236 in all samples of these three stages (Fig. 2d). Overall, U6-1 and EIF4A could be chosen as the most sable reference genes because of their lower M-value in various tissues at different developmental stages, whereas miR403 and miR171 were the least stable genes with increased M-values (from 0.577 to 0.955) in all samples. In the subsets of different developmental stages, all the pairwise variation values of V2/3 were less than 0.15, which suggests that the combined use of the two most stable reference genes would be most effective for normalizing gene expression analysis (Fig. 3).
The stability values of the candidate reference genes constructed by NormFinder are presented in Table 2. At the callus stage, Histone was the most stable gene with the lowest stability value followed by EIF4A and U6-2. U6-1, U6-2 and EIF4A exhibited more stable expression at the seedling stage, whereas UBP, bHLH and DNAJ exhibited increased stability in various tissues at the plant stage. In all samples, bHLH, PP2A-2 and U6-1 were the top three stable genes. miR403, miRN482 and miR171 were the least stable genes during the plant regeneration process.
The rankings using the Delta Ct method are presented in Table 4. EIF4A, U6-2 and U6-1 were the top three stable genes at the callus and seedling stages. PP2A-2 was the most stable gene at the plant stage and in all samples during these three developmental stages, followed by U6-1. This finding indicates that EIF4A, PP2A-2 and U6-1 might be the most stably expressed genes as determined using the Delta Ct method.
To obtain a consensus regarding the most stable reference genes as recommended by the four methods, the geometric mean of four algorithms corresponding rankings for each candidate gene were calculated using the RankAggreg software. EIF4A, U6-2 and U6-1 were ranked as the top three stable reference genes in samples from the callus developmental stage. U6-1 was the most stable gene at the seedling and plant stages, and U6-2 and PP2A-2 ranked second at the seedling stage and plant stage respectively. The expression values of miR403, miR171 and miR482 were extremely variable in all tissues at different stages. Based on these results, U6-1, EIF4A and PP2A-2 are the best combination of reference genes in all samples of different developmental stages. EIF4A and U6-2 are the most stable reference genes for the samples from callus stage. U6-1 and U6-2 are the best reference genes for the seedling stage, and U6-1 and PP2A-2 suitable for the plant stage.   (Table 5).

Reference genes validation
To validate the stability of the above reference genes, the expression profiles of miR166a, which is highly expressed in leaves and xylem [33], and PtHB4, which expressed specifically in shoot apex, cambium and xylem [34,35], were measured and normalized with the most stable and the least stable reference genes. The three top ranking reference genes (U6-1, EIF4A, Fig. 3 The optimal number of reference genes for accurate normalization calculated by geNorm during plant regeneration. Pairwise variation (V n / V n+1 ) analysis of 18 candidate reference genes analyzed in four sample subsets. callus, callus and shoots induced by 84K leaves; seedling, various tissues from 84K seedling stage; plant, various tissues from 84K plant stage; ALL, all samples from callus, seedling and plant developmental stages PP2A-2) either alone or combination and two unstable reference genes (5.8 s and miR403), were used as reference genes for RT-qPCR analysis at callus, seedling and plant stages. As shown in Fig. 4, the relative expression profiles of miR166a and PtHB4 normalized with U6-1, EIF4A, PP2A-2, U6-1 + EIF4A and U6-1 + PP2A-2    exhibited perfect consistency in seventeen tissues. Compared with tender leaves from 0.5-month-old seedlings, miR166a was highly expressed in stems, mature leaves and all samples at the callus stage (Fig. 4), whereas PtHB4 was specifically expressed in stems and samples at the callus stage (Fig. 4). However, the expression values of miR166a and PtHB4 normalized with 5.8s were increased, and the values normalized to miR403 were reduced compared with than other reference genes (Fig. 4). For example, the relative expression value of miR166 in mature leaves at 3 months was approximately 1500 when normalized with U6-1. However, the expression value was 4500 when normalized with 5.8 s and 18 when normalized with miR403. Because the combined use of two most stable reference genes would be suitable for normalizing gene expression analysis at the AR developmental stage, PP2A-2, UBP, bHLH, PP2A-2 + UBP and PP2A-2 + bHLH were selected as stable reference genes, and 5.8s and miR403 were selected as the two least stable reference genes for RT-qPCR analysis of the AR developmental stage. As shown in Fig. 5, the relative expression values of miR166 and PtHB4 normalized with PP2A-2 were higher than values normalized with bHLH, but similar expression patterns were noted. Therefore, the combination of PP2A-2 and bHLH could neutralize the expression values normalized with PP2A-2 and bHLH. miR166 was highly expressed in AR-60H, in which the AR callus regenerated and expanded. PtHB4 exhibited a higher expression value in AR-18H when AR induction had begun. The relative expression profiles of miR166 and PtHB4 normalized with 5.8s and miR403 were reduced compared with values normalized with other reference genes, and different expression patterns were noted (Fig. 5). Overall, the combination of, PP2A-2 and UBP or PP2A-2 and bHLH should be the best reference gene set for normalization of qRT-PCR at the AR developmental stage.

Discussion
Numerous studies have performed reference gene validation experiments for mRNA and miRNA qRT-PCR [36,37]. For example, some house-keeping genes, such as ACT 7, UBQ, GAPDH and TUB, were widely used for gene expression analysis of mRNAs in diverse plants, and several noncoding RNAs, such as U6 snRNA and 5.8S rRNA, were typically chosen for normalizing miRNA quantification data. As miRNA research continues to expand, the potential use of miRNAs as reference genes has attracted increasing attention, and some small RNAs were obtained from plant species, such as grapevine [7], wheat [38], peach [39], soybean [40], and tea [41] demonstrating that miRNAs are more stable than the currently used reference genes under specific conditions. However, there are no reports on universal reference genes for the quantitative expression of both mRNA and miRNA, which may be necessary to calibrate both miRNA and its target gene(s) in a given sample. In general, the reverse transcription of miRNA is different from other types of RNA because of the shorter length of miRNA. Therefore, quantitative expression analysis of miRNAs and mRNAs requires different internal reference genes for normalization of their respective transcripts [15,42]. If miRNA and mRNA expression data are normalized with universal internal reference genes, quantitative PCR must be performed at the same transcriptional level. Fortunately, Hurteau developed a modified universal reverse transcription PCR protocol that is designed to specifically amplify and quantify mRNAs and miRNAs from the same sample [17]. The modified technique involves the enzymatic addition of a poly A tail to non-poly(A) RNAs followed by reverse transcription using a universal RTprimer. Then, the transcript-specific forward primers can be used to amplify noncoding RNA (including miRNA) and mRNA from the same sample [17]. To obtain suitable universal reference genes for normalization of miRNA and mRNA expression, the expression stability of the selected candidate reference genes in this study was validated across 38 tissue samples from four developmental stages in 84K poplar using a modified universal reverse transcription PCR protocol. Four algorithms (geNorm, NormFinder, DeltaCt, Bestkeeper) were used to minimize the bias for the evaluation of the 18 candidate reference genes. Discrepancy was observed in gene stability ranking and validation generated by the four different algorithms above. For example, at the plant stage, PP2A-2 was ranked first by geNorm and DeltaCt, whereas it was ranked fifth by Nor-mFinder and Bestkeeper. Histone was ranked among the top four stable genes by BestKeeper in all samples at the seedling and plant stages but was ranked in the middle or bottom position by geNorm, NormFinder and Del-taCt. This apparent divergence is probably due to the statistical algorithms used to calculate stability. Genorm and Normfinder software have similar algorithms, which calculate the ΔCt values and the stability of each internal reference gene according to the minimum Ct value in all the samples [43,44]. Therefore, the alteration of the Fig. 4 Relative expression of miR166 and PtHB4 using validated reference genes for normalization under plant regeneration process. The validated reference gene(s) used as normalization factors were one (U6-1, EIF4A, PP2A-2 alone) or two (U6-1 + EIF4A and U6-1 + PP2A-2) of three best stable reference genes and two most unstable reference genes (5.8s and miR403) Fig. 5 Relative expression of miR166 and PtHB4 using validated reference genes for normalization under AR regeneration. PP2A-2, UBP, bHLH, PP2A-2 + UBP and PP2A-2 + bHLH were selected for the stable reference genes, and 5.8 s and miR403 were the two least stable genes for normalization RT-qPCR analysis minimum Ct value will have a great influence on the stability of expression of candidate reference genes. Besides that, BestKeeper program can calculate the standard deviation (SD) and variation coefficient (CV) values between the Ct values of each internal reference gene and the average of all Ct values, so the stability of candidate reference gene depending on the dispersion degree of Ct values [45]. Thus, it has been recommended that more than two algorithms should be used for reference gene stability evaluation [43]. In this study, RankAggreg software was used to generate the final overall ranking of the tested reference genes based on the geometric mean of the weights of every gene calculated by each program. Through this comprehensive ranking analysis, U6-1, EIF4A and PP2A-2 ranked as top three for all samples at the callus, seedling and plant stages. Thus, these genes would be the most suitable internal reference genes for the quantification of mRNAs and noncoding RNAs during the plant regeneration process of poplar.
However, the candidate reference genes selected for samples during the AR developmental stage were almost completely different. UBP ranked among the top three most stable genes in AR developmental stage but near the bottom at the callus stage using the four algorithms. In contrast, miR482 was the most stable gene in the AR developmental stage identified by Bestkeeper, but it was a highly unstable gene in other developmental stages. This is mainly because the development of adventitious roots is quite different from other biological processes, which include root primordium initiation, callus differentiation, adventitious root emergence and elongation process [46,47]. During this process, tissue morphology and structure have undergone dramatic changes, and gene expression patterns varied tremendously [48][49][50]. Therefore, the genes expressed stably in other tissues or developmental stages may not have stable expression values during AR developmental stages. Comparatively, PP2A-2, UBP and bHLH exhibited more stable expression in the AR developmental stage; thus, these genes could be suitable for normalization. These results indicate that different sets of internal reference genes may be assigned for different tissues and developmental stages even in the same species. Furthermore, 5.8S, which is a commonly used reference gene for miRNA, and commonly used miR-NAs (miR171, miR482 and miR403) exhibited extremely unstable expression and were ranked among the four least stable genes during all developmental stages in poplar, indicating these miRNAs were more unstable than protein-coding genes and exhibited obvious tissue-specific expression.
Traditionally, reference genes are typically cellular maintenance genes that play housekeeping roles in basic cellular components and functions [23]. In this study, with the exception of two house-keeping genes (EIF4A and PP2A-2), two noncoding RNA (U6-1 and U6-2) and a nontraditional reference gene (UBP) were also confirmed as the most suitable internal reference genes at the seedling and AR developmental stages. EIF4A and PP2A are typically used as reference genes for quantitative PCR and exhibit stable expression in different tissues, developmental stages or biotic and abiotic stress conditions in a number of species. For example, PP2A and EIF4A were reported as the best reference genes for all samples of various tissues and abiotic stress conditions in Sorghum [51], and PP2A was suitable for Switchgrass [52]. EIF4α was also ranked as a stably expressed gene under most of the experimental conditions tested in Carica papaya [53], different tissue/organs and fruit developmental stages in Litsea cubeba [54], and different tissues under abiotic stresses in Pennisetum glaucum [55]. PP2A housekeeping genes were superior references for normalization of gene expression data in different cotton plant organs [56], different color lines of cineraria during flower developmental stages [57] and diurnal and developmental time-course in lettuce [58]. In addition, U6 is also one of the most commonly used reference genes in miRNA qRT-PCR and had the most stable expression in reference gene selection studies [7,25,41]. In our study, U6-1 and U6-2 exhibited the most stable expression values at the callus and seedling stages, and U6-1 ranked among the top two at the plant stage. Therefore, as a noncoding RNA, U6 can be used not only as a good reference gene for RT-qPCR of miRNA alone but also a universal internal reference gene for mRNA quantification. The UBP gene encodes a heterogeneous nuclear RNA binding protein (hnRNP) that is involved in the regulation of pre-mRNA maturation at different levels and pre-mRNA splicing [59]. UBP is not a traditional internal reference gene, but its expression value did not change in various tissues and developmental stages in poplar [60]. UBP was the most stably expressed gene in the AR regeneration stage, and the stability of UBP was ranked at the top and middle positions at the plant and seedling stages. This finding indicates that UBP might be a reliable reference gene used in the AR and plant developmental stages in poplar.

Conclusion
The purpose of this study was to identify the most appropriate reference genes for qRT-PCR of miRNAs and mRNAs during poplar regeneration and development. The expression stability of 18 candidate genes was validated and evaluated across 38 tissue samples from four developmental stages of 84K poplar using four algorithms. The results demonstrated that EIF4A and U6-2 were suitable for samples of callus stage, U6-1 and U6-2 were best for seedling stage samples. PP2A-2 and U6-1 were best for the plant stage, and U6-1, EIF4A and PP2A-2 were the top three reference genes during the plant regeneration process in poplar. In addition, PP2A-2 and UBP or PP2A-2 and bHLH were the best combination as reference genes in the AR regeneration stage. This work will benefit future studies of expression and function analysis of miRNAs and their target genes in poplar.

Plant materials and tissue harvesting
The seedlings and plants of 84K poplar (Populus alba × Populus glandulosa) were grown in a tissue culture room under long-day conditions (16 h light/8 h dark) at 25/22 °C (day/night). Based on the Agrobacterium-mediated leaf disk transformation method [31,32], the regeneration processes of poplar plants can be divided into three growth stages, including callus induction and shoots differentiation, seedlings on culture medium and plants in soil. Callus and shoot samples at different developmental stages were collected during shoot regeneration processes of 84K poplar [61], including callus induction stage at 4 days (Ca-4D), callus proliferation stage at 6 days (Ca-6D), callus expansion stage at 10 days (Ca-10D), callus transition stage at 12 days (Ca-12D), shoot emergence stage at 16 days (Shoot-16D) and the shoot elongation stage at 20 days (Shoot-20D). The samples from seedlings on the culture medium included the leaves (L-0.5M), stems (S-0.5M) and roots (R-0.5M) of 0.5-month-old seedlings as well as the leaves (L-1M), stems (S-1M) and roots (R-1M) of 1-month-old seedlings. The tissues from 3-month-old plants grown in soil included the shoot apical meristem (SAM-3M), unexpanded leaves (YL-3M), the first and second expanded leaves (ML-3M), the first to thirteenth internodes (N1-3M ~ N13-3M), the stem base (Base-3M), the roots (R-3M), and the root tips (RT-3M). The samples from the adventitious root (AR) developmental process included the AR induction stage at 18 h (AR-18H), the AR callus regeneration stage at 42 h and 60 h (AR-42H and AR-60H), the AR emergence stage at 4 and 5 days (AR-4D and AR-5D) and the AR elongation stage at 7 and 11 days (AR-7D and AR-11D) [61,62]. Samples with three replicates were collected, immediately frozen and stored in liquid nitrogen.

Total RNA extraction and cDNA synthesis
Total RNA was extracted using the LC sciences Total RNA Purification Kit (#TRK-1001, LC sciences, USA), which purifies all sizes of total RNA, including mRNA, ribosomal RNA, miRNA and other small RNA (20-200 nt), according to the previous methods with some modification [63]. The powder ground from 50 mg sample in liquid nitrogen was immediately transferred into a 2.0 ml RNase-free tube and added 600 μl extraction buffer with 6% Plant RNA Isolation Aid (Ambion, #Am9690). After shaking vigorously, the mixture was incubated in the ice for 15 min and then centrifuged at 12,000 rpm for 10 min at room temperate, by which the yield of total RNA could be improved. The remaining processes followed the manufacturer's instructions. The integrity of total RNA was further assessed by 1.5% agarose gel electrophoresis, and the RNA concentration and purity were determined by NanoDrop ™ 8000 Spectrophotometer (Thermo Fisher Scientific, USA). Only RNA samples with an A260/ A280 ratio between 1.9 and 2.1 and A260/A230 greater than 1.80 were used for cDNA synthesis. Then, 1.5 μg of total RNA was polyadenylated with ATP by poly (A) polymerase (PAP) at 37 °C for 1 h in a 20-μl reaction mixture using the Poly(A) Tailing Kit (#AM1350, Invitrogen, USA). Then, 10 μl (750 ng) of the E-PAP-treated total RNA was reverse transcribed with a poly(T) adapter universal reverse transcription (RT)-primer (5′-AAC GAG ACG ACG ACA GAC TTT TTT TTT TTT TTTV-3′) using SuperScript III reverse transcriptase Kit (#18080-051, Invitrogen, USA) following the manufacturer's instruction. The cDNA was diluted 20-fold with nuclease-free water for qRT-PCR.

Quantitative real-time PCR analysis
Quantitative real-time PCR was conducted using LightCycler ® 96 Plates and performed on the LightCycler ® 480 System (Roche Molecular Systems, Germany). The reaction mixture contained 10 μl KAPA SYBR FAST qPCR Master Mix (# K4601, KAPA Biosystems, USA), 2 μl 20-fold diluted cDNA, 0.4 µM of each forward and reverse primer ( Table 1) and ddH 2 O in a final volume of 20 μl. Amplifications were performed with the following program: 95 °C for 3 s; 40 cycles of 95 °C for 10 s, 60 °C for 30 s, 72 °C for 3 s; and melting curve analysis conditions (95 °C for 5 s, 65 °C increased to 95 °C with temperature increment of 0.11 °C every 1 s). No-template reactions were used as negative controls, and each sample was assessed in four technical replicates. Using a series of 10-fold diluted cDNA as templates, the standard curves were generated for each candidate reference gene. The correlation coefficient (R 2 ) and slope were obtained from the linear regression model created by the LightCycle 480 system, and the PCR amplification efficiency (E) was calculated using the formulas E = 10 −1/slope − 1.

Stability analysis for the candidate reference genes
To visualize the expression stability of the 18 candidate reference genes, the raw cycle threshold (Ct) values from different tissues and developmental stages were produced and calculated statistically by box plots and five different programs and algorithms, including geNorm [43], Nor-mFinder [44], BestKeeper [45], the Delta CT method [66,67] and the RankAggreg software [68]. The geNorm algorithm can calculate the average expression stability value (M), which is defined as the average pairwise variation in a particular gene with all other potential reference genes. The threshold of M value was set as 1.5, and genes with lowest M values exhibit the most stable expression. Additionally, geNorm also calculates the pairwise variation (V n + V n+1 ) value that determines the optimal number of reference genes for accurate normalization with a cut-off value of V n+1 < 0.15 [43]. The NormFinder program uses an ANOVA-based model to consider intra-and intergroup variation in expression levels to calculate a stability value (SV) for expression, and a lower SV indicates increased stability [44]. Bestkeeper is an excel-based tool that determines the stability ranking of reference genes based on the coefficient of variance (CV) and the standard deviation (SD) of the average Ct values. The most stable gene exhibits the lowest CV ± SD value, and genes with SD higher than 1 were considered unacceptable and were excluded [45]. The Delta Ct method compares the relative expression of 'pairs of genes' within each sample. The stability of the reference gene is ranked according to a 'process of elimination' technique, by which genes can be compared against one another and either selected or eliminated on the basis of ΔCt among samples [66]. Finally, the raw Ct values of each gene were used to calculate the comprehensive ranking of reference genes using the RankAggreg software [68], which was based on the ranking of candidate references obtained from the four programs mentioned above. The program assigns an appropriate weight to an individual gene and calculates the geometric mean of the weight, providing an overall comprehensive ranking.

Validation of identified reference genes
To examine the expression stability of potential reference genes, the relative expression levels of miR166 and its target gene PtHB4 were analyzed in various tissues and AR developmental stages in poplar. The relative expression levels were normalized separately to the most stable and least stable reference genes analyzed by the four algorithms. The qRT-PCR amplification conditions of miRNA and genes were the same as described above. The relative expression levels of these genes were calculated according to the 2 −∆∆Ct method and presented as foldchange [69].
of Horticulture and Forestry Sciences, Huazhong Agricultural University, Wuhan 430070, China.