Identification of Appropriate Reference Genes for Normalizing miRNA Expression in Citrus Infected by Xanthomonas citri subsp. citri

MicroRNAs (miRNAs) are short noncoding RNA molecules that regulate gene expression at the posttranscriptional level. Reverse transcription-quantitative PCR (RT-qPCR) is one of the most common methods used for quantification of miRNA expression, and the levels of expression are normalized by comparing with reference genes. Thus, the selection of reference genes is critically important for accurate quantification. The present study was intended to identify appropriate miRNA reference genes for normalizing the level of miRNA expression in Citrus sinensis L. Osbeck and Citrus reticulata Blanco infected by Xanthomonas citri subsp. citri, which caused citrus canker disease. Five algorithms (Delta Ct, geNorm, NormFinder, BestKeeper and RefFinder) were used for screening reference genes, and two quantification approaches, poly(A) extension RT-qPCR and stem-loop RT-qPCR, were used to determine the most appropriate method for detecting expression patterns of miRNA. An overall comprehensive ranking output derived from the multi-algorithms showed that poly(A)-tailed miR162-3p/miR472 were the best reference gene combination for miRNA RT-qPCR normalization in citrus canker research. Candidate reference gene expression profiles determined by poly(A) RT-qPCR were more consistent in the two citrus species. To the best of our knowledge, this is the first systematic comparison of two miRNA quantification methods for evaluating reference genes. These results highlight the importance of rigorously assessing candidate reference genes and clarify some contradictory results in miRNA research on citrus.


Introduction
Citrus is one of the most important fruit crops due to its high nutritional and economic value. Citrus crops are cultivated in a wide range of regions with a total fruit production over 130 million tons in 2015, ranking first in quantity in world fruit production (FAO, 2017). Citrus production, however, has encountered several serious diseases, including citrus canker and Huanglongbing [1]. Citrus canker is caused by X. citri subsp. citri (Xcc) with symptoms of water-soaked eruptions, circular lesions, and pustule-like lesions on all plant tissues, resulting in heavy economic losses [2,3].
For the use of the poly-A tail method, each polyadenylation and reverse transcription reaction was performed with 1 µg of total RNA in a final volume of 10 µL using a Mir-X™ miRNA First Strand Synthesis kit (Takara/clontech, Mountain View, CA, USA) where there were 5 µL mRQ Buffer (2×) and 1.25 µL mRQ Enzyme. The reaction tubes were incubated in a thermocycler for 1 h at 37 • C and terminated at 85 • C for 5 min to inactivate the enzymes.
For the stem-loop reverse transcription method, the experimental RNA (up to 1 µg), 1 µL (1 µM) stem-loop specific primer (Table S1), and 1 µL (50 µM) Oligo(dT) 15 primer (Promega, Madison, WI, USA) were mixed in nuclease-free water for a final volume of 5 µL per RT reaction and incubated at 70 • C for 5 min followed by ice-cooling. cDNA synthesis was carried out using Improm-II reverse transcription system according to the manufacturer's instruction (Promega). Briefly, the reaction thermal profile was set with an initial temperature of 16 • C for 30 min followed by 60 cycles of 30 • C for 30 s, 42 • C for 30 s and 50 • C for 1 s. After 60 cycles, the reaction was terminated by incubating the samples at 85 • C for 5 min [28].

Selection of Candidate Reference Genes and Primer Design
Our previous study identified the transcription levels of conserved miRNAs from sweet orange and Ponkan mandarin infected with Xcc using the high-throughput sequencing. NEB Next Ultra small RNA Sample Library Prep Kit for Illumina (New England, NEB, USA) was employed for generating all 12 sequencing libraries. The Illumina HiSeq2500 platform was used to sequence the library, thereby generating a total of 260,117,346 raw reads. The transcription levels of 20 conserved miRNAs (including published reliable references miR171 [32], miR166b [31], miR159a [49], miR396a [50] and miR166b [51] were obtained and analyzed. Eight putative reference RNA targets (miR472, miR428b, miR396a, miR166b, miR3954, miR160, miR162-3p and miR403) with a median expression and the smallest standard deviations were selected from high-throughput sequencing data for identifying the most suitable normalizer for miRNA qPCR in citrus infected by Xcc.
For the stem-loop method, the primers (Table S1) for miRNAs cDNA synthesis and forward miRNA primers for qPCR were designed according to the report of Chen et al. (2005) [28] or by using the online tool (http://150.216.56.64/miRNAprimerDesigner.php). Universal reverse primer (Table 1) was also based on the report of Chen et al. (2005) [28]. For Poly-A tail primer design, miRNA specific 5' primer was used as the entire sequence of mature miRNA (21-23 nt). The 3' primer for qPCR was the mRQ 3' Primer supplied with the Mir-X™ miRNA First-Strand Synthesis kit.
Non-coding RNA (U4, U5, U6 and snoR14) primers (Table 1) were obtained from previous research [52]. Eight protein-coding traditional reference genes (ACTIN1, ACTIN2, UBC28, TUA5, EF1α, TUB4, GAPDH and PP2A) reported to be good potential candidates in previously published studies [53][54][55][56] also included in this study. The sequences of citrus protein-coding traditional candidate reference genes and genes coding for the CsLOB1 transcription factor were extracted from citrus genome data (http://citrus.hzau.edu.cn/orange/index.php). Primer design was made according to a standard set of design criteria (e.g., primer Tm, length and GC content), which generated a unique, short PCR product of the expected length and sequence [57].

Reverse Transcription-Quantitative PCR Analysis
The reaction mixture with a total volume of 20 µL contained 2 µL diluted cDNA, 0.4 µL of each 10 µM primer, 7.2 µL ddH 2 O and 10 µL SYBR Advantage Premix (Takara). Four biological replicates for each sample together with two technical replicates for each well were performed. The RT-qPCR conditions were as follows: denaturation at 95 • C for 10 s, followed by 40 cycles of denaturation at 95 • C for 5 s, and annealing at 60 • C for 30 s in CFX96 Touch™ Real-Time PCR detection system (Bio-Rad, Hercules, CA, USA). After amplification, a melting curve analysis (60-95 • C with at increments of 0.5 • C) was generated for each reaction to verify the specific amplification.
The four common statistical approaches (Delta-Ct method, geNorm, NormFinder and BestKeeper) might yield contradictory results from the same data. However, RefFinder could make comprehensive stability ranking of candidate genes according to the results from these softwares. Eight miRNA genes were analyzed in a comparative way to determine the most reliable normalizer, data from stem-loop RT and poly(A) extension were imported to RefFinder individually. To evaluate the potentiality of these candidate genes as miRNA reference gene, we also analyzed their expressional stability by combining data from protein coding RNA.

Reference Gene Validation
A citrus canker-induced gene Lateral Organ Boundaries (CsLOB1) was analyzed in C. sinensis by RT-qPCR [63], CsLOB1 primers was designed according to the previously mentioned methods, and the same conditions and criteria as for the reference genes were used in CsLOB1. To calculate the relative fold changes in the gene expression of sweet orange leaves after Xcc inoculation, the relative expression of the target was calculated using the comparative 2 −∆∆ Ct method and normalized to the selected reference genes.

Verification of Primer Efficiency for the Candidate Reference Genes
The correlation coefficient (r) represents whether the plotted data points fit the regression line generated by the standard curve. As a good linear performance, r should be more than 0.99. In this study, r of the standard curve ranged from 0.9730 to 0.9998 (Table 1). The amplification efficiency (E) refers to amplification rate of a template during the reaction. The theoretically ideal E value is 100% indicating that the target cDNA template duplicates exponentially in every PCR cycle [64]. The amplification efficiency of the stem-loop primers varied from 86.8% to 104.7%, and the poly(A)-based reaction E value ranged from 91.5% to 103.6% (Table 1). In practical terms, primer efficiency between 90% and 110% was thought to be acceptable for qPCR experiments [65]; thus, all the primer pairs were deemed sufficiently designed. All of the non-coding RNAs showed low amplification efficiency, which was not fit for qPCR experiment.

Ct Value Distribution of Candidate Reference Genes
The cycle threshold (Ct) values of most genes ranged from 15 to 30, the non-coding RNA (U4, U5 and U6) showed relatively high expression levels in both sweet orange and Ponkan, with the Ct values lower than 15 ( Figure 1). St-miR160 had a highest Ct value of 36.87, which was not recommended (15 ≤ Ct ≤ 30) for a RT-qPCR [66]. St-miR472 showed the narrowest Ct ranges (21.61 to 23.75) in sweet orange ( Figure 1A), whereas the non-coding RNA gene U5 was the most variable one. Moreover, the Ct value of the stem-loop primer was more variable than poly(A) extension microRNA in all the tested samples. The variability represents a difference between max Ct and min Ct. St-miR162-3p and St-miR403 showed high variability in sweet orange, with their Ct values ranging from 25.83 to 31.31 and 21.61 to 26.97, respectively ( Figure 1A). Similarly, St-miR160 showed the highest Ct value in Ponkan ( Figure 1B), while Ct variance of miR162-3p and miR472 was minimal. In Ponkan plants where the Ct value ranged from 15 to 20, the expression was relatively stable.

Reference Gene Stability Analyzed by geNorm, NormFinder, BestKeeper, Delta-Ct Method, and Refinder
To analyze the expression stability of candidate reference genes during different infection stages of citrus leaves, four statistical algorithms (Delta-Ct method, BestKeeper, NormFinder, and geNorm) and a web-based tool (RefFinder) were used to rank the stability of 28 reference genes including 8 miRNA genes reverse transcribed by stem-loop method ( Table 2).
All the Ct raw data were directly input to the program BestKeeper. The Ct value of genes with standard deviation (SD) < 1 was considered to be stable enough or SD (± x-fold) < 2 was also a statistical parameter to determine the reference genes. A matrix of pairwise comparisons and r calculation was performed using BestKeeper. Genes with higher expression stability should have a r closer to 1. Based on the BestKeeper analysis, U5, TUA5, GAPDH and st-miR162-3p could not pass the filter with SD values of 1.44, 1.17, 1.74 and 1.12, respectively, which were much higher than 1 in sweet orange (Table 3). On the other hand, st-miR472 and miR160 (SD: 0.42 and 0.48) were confirmed to be the most stably expressed. Meanwhile, the BestKeeper calculation also provided two best reference genes miR162-3p and miR472 with low SD values (0.25 and 0.26). GAPDH was the least stably-expressed gene for the two citrus species with r values over 25 (Tables 3 and 4).
The geNorm is a popular algorithm to determine the most stable reference genes from a set of tested candidate reference genes in a given sample panel. We identified miR166b, UBC28, miR472 and miR162-3p as the most stable genes in sweet orange according to the average expression stability value M using geNorm software (Figure 2A). In geNorm, the lower M value represents the higher stability. It is noteworthy to note that all the reference genes in sweet orange were deemed unstable under citrus canker infection according to the commonly proposed cutoff value of M ≤ 0.5. Using a M threshold value of ≤ 0.5, miR472, miR162-3p, miR427b, miR403 and miR160 had the lower M value (0.2790, 0.2790, 0.3930, 0.4450 and 0.4720, respectively) indicating these genes were the most stably expressed genes in Ponkan under canker infection ( Figure 2B). The protein coding gene GAPDH was the least stably expressed gene in two citrus species.

Reference Gene Stability Analyzed by geNorm, NormFinder, BestKeeper, Delta-Ct Method, and Refinder
To analyze the expression stability of candidate reference genes during different infection stages of citrus leaves, four statistical algorithms (Delta-Ct method, BestKeeper, NormFinder, and geNorm) and a web-based tool (RefFinder) were used to rank the stability of 28 reference genes including 8 miRNA genes reverse transcribed by stem-loop method ( Table 2).
All the Ct raw data were directly input to the program BestKeeper. The Ct value of genes with standard deviation (SD) < 1 was considered to be stable enough or SD (± x-fold) < 2 was also a statistical parameter to determine the reference genes. A matrix of pairwise comparisons and r calculation was performed using BestKeeper. Genes with higher expression stability should have a r closer to 1. Based on the BestKeeper analysis, U5, TUA5, GAPDH and st-miR162-3p could not pass the filter with SD values of 1.44, 1.17, 1.74 and 1.12, respectively, which were much higher than 1 in sweet orange (Table 3). On the other hand, st-miR472 and miR160 (SD: 0.42 and 0.48) were confirmed to be the most stably expressed. Meanwhile, the BestKeeper calculation also provided two best reference genes miR162-3p and miR472 with low SD values (0.25 and 0.26). GAPDH was the least stably-expressed gene for the two citrus species with r values over 25 (Tables 3 and 4).
The geNorm is a popular algorithm to determine the most stable reference genes from a set of tested candidate reference genes in a given sample panel. We identified miR166b, UBC28, miR472 and miR162-3p as the most stable genes in sweet orange according to the average expression stability value M using geNorm software ( Figure 2A). In geNorm, the lower M value represents the higher stability. It is noteworthy to note that all the reference genes in sweet orange were deemed unstable under citrus canker infection according to the commonly proposed cutoff value of M ≤ 0.5. Using a M threshold value of ≤ 0.5, miR472, miR162-3p, miR427b, miR403 and miR160 had the lower M value (0.2790, 0.2790, 0.3930, 0.4450 and 0.4720, respectively) indicating these genes were the most stably expressed genes in Ponkan under canker infection ( Figure 2B). The protein coding gene GAPDH was the least stably expressed gene in two citrus species.
Additionally, geNorm is also based on a pairwise comparison model to determine the optimal number of reference genes for accurate data normalization. Average pairwise variation Vn/n + 1 with a cut-off score of 0.15 determined the minimum number of genes necessary for data normalization, below which the inclusion of an additional reference gene had no significant contribution to the normalization. For reference genes expression in sweet orange, the analysis revealed that the V3/4 value was 0.168, the normalization factor should preferably contain at least the four best reference genes. While the V4/5 value was 0.126 ( Figure 3A), the fifth additional gene was not necessarily needed to normalize the expression, which confirmed the geNorm results. In all of the Ponkan sample sets, the V2/3 value was considerably smaller than the default threshold value of 0.15, suggesting that the two most stable reference genes were miR162-3p and miR472, which should be appropriate for normalization ( Figure 3B).
We further used NormFinder for identifying appropriate reference genes. NormFinder is a Visual Basic Application for Excel. Like geNorm, the software calculates the expression stability value (M) for all candidate genes. Using an analysis of variance (ANOVA)-based comparison model, NormFinder can estimate the variation of the intra-and inter-group. Based on these two statistical methods, the lowest M always suggests the most reliable pair of reference genes. In the sweet orange sample set, the NormFinder analysis revealed that the most stably expressed genes were miR162-3p and miR472, followed by miR160 ( Figure 4A). Furthermore, the analysis indicated GAPDH as the least stably expressed gene. The algorithm also selected an optimal pair of reference genes in Ponkan, and the most stable ones were miR403 and miR428b with an M value of 0.28 and 0.45 ( Figure 4B).  TUB4  U6  TUB4  TUB4  26  TUA5  TUA5  TUA5  TUA5  TUA5  ACTIN1  U5  ACTIN1  ACTIN1  ACTIN1  27 U5    Additionally, geNorm is also based on a pairwise comparison model to determine the optimal number of reference genes for accurate data normalization. Average pairwise variation Vn/n + 1 with a cut-off score of 0.15 determined the minimum number of genes necessary for data normalization, below which the inclusion of an additional reference gene had no significant contribution to the normalization. For reference genes expression in sweet orange, the analysis revealed that the V3/4 value was 0.168, the normalization factor should preferably contain at least the four best reference genes. While the V4/5 value was 0.126 ( Figure 3A), the fifth additional gene was not necessarily needed to normalize the expression, which confirmed the geNorm results. In all of the Ponkan sample sets, the V2/3 value was considerably smaller than the default threshold value of 0.15, suggesting that the two most stable reference genes were miR162-3p and miR472, which should be appropriate for normalization ( Figure 3B). Ponkan. The pairwise variation (Vn/Vn+1) was analyzed based on geNorm algorithm to determine the optimal number of reference genes for accurate normalization. We proposed 0.15 as a threshold value, which suggested that adding one more gene into the combination of reference genes is not required.
We further used NormFinder for identifying appropriate reference genes. NormFinder is a Visual Basic Application for Excel. Like geNorm, the software calculates the expression stability value (M) for all candidate genes. Using an analysis of variance (ANOVA)-based comparison model, NormFinder can estimate the variation of the intra-and inter-group. Based on these two statistical methods, the lowest M always suggests the most reliable pair of reference genes. In the sweet orange  Additionally, geNorm is also based on a pairwise comparison model to determine the optimal number of reference genes for accurate data normalization. Average pairwise variation Vn/n + 1 with a cut-off score of 0.15 determined the minimum number of genes necessary for data normalization, below which the inclusion of an additional reference gene had no significant contribution to the normalization. For reference genes expression in sweet orange, the analysis revealed that the V3/4 value was 0.168, the normalization factor should preferably contain at least the four best reference genes. While the V4/5 value was 0.126 ( Figure 3A), the fifth additional gene was not necessarily needed to normalize the expression, which confirmed the geNorm results. In all of the Ponkan sample sets, the V2/3 value was considerably smaller than the default threshold value of 0.15, suggesting that the two most stable reference genes were miR162-3p and miR472, which should be appropriate for normalization ( Figure 3B). Ponkan. The pairwise variation (Vn/Vn+1) was analyzed based on geNorm algorithm to determine the optimal number of reference genes for accurate normalization. We proposed 0.15 as a threshold value, which suggested that adding one more gene into the combination of reference genes is not required.
We further used NormFinder for identifying appropriate reference genes. NormFinder is a Visual Basic Application for Excel. Like geNorm, the software calculates the expression stability value (M) for all candidate genes. Using an analysis of variance (ANOVA)-based comparison model, NormFinder can estimate the variation of the intra-and inter-group. Based on these two statistical methods, the lowest M always suggests the most reliable pair of reference genes. In the sweet orange The pairwise variation (Vn/Vn+1) was analyzed based on geNorm algorithm to determine the optimal number of reference genes for accurate normalization. We proposed 0.15 as a threshold value, which suggested that adding one more gene into the combination of reference genes is not required. sample set, the NormFinder analysis revealed that the most stably expressed genes were miR162-3p and miR472, followed by miR160 ( Figure 4A). Furthermore, the analysis indicated GAPDH as the least stably expressed gene. The algorithm also selected an optimal pair of reference genes in Ponkan, and the most stable ones were miR403 and miR428b with an M value of 0.28 and 0.45 ( Figure 4B). To obtain the most suitable reference genes, we evaluated the candidates by the Delta Ct method. This approach was employed by comparing the mean difference of 'gene pairs' and the variation in Ct values within each sample. The ranking of the reference genes from this method was similar to that obtained by the NormFinder algorithm. Delta Ct method analysis found that miR162-3p and miR472 had lowest StdDev, followed by miR160 which suggested that miR162-3p and miR472 were two top-ranked, stable reference genes among all of the studied genes in sweet orange ( Figure 5A), which also confirmed the results from NormFinder analysis. For Ponkan (Figure 5B), miR403 To obtain the most suitable reference genes, we evaluated the candidates by the Delta Ct method. This approach was employed by comparing the mean difference of 'gene pairs' and the variation in Ct values within each sample. The ranking of the reference genes from this method was similar to that obtained by the NormFinder algorithm. Delta Ct method analysis found that miR162-3p and miR472 had lowest StdDev, followed by miR160 which suggested that miR162-3p and miR472 were two top-ranked, stable reference genes among all of the studied genes in sweet orange ( Figure 5A), which also confirmed the results from NormFinder analysis. For Ponkan (Figure 5B), miR403 followed by miR482b and miR162-3p, were found to be the most stably expressed genes with low StdDev (0.09, 1.04 and 1.04 respectively). To obtain the most suitable reference genes, we evaluated the candidates by the Delta Ct method. This approach was employed by comparing the mean difference of 'gene pairs' and the variation in Ct values within each sample. The ranking of the reference genes from this method was similar to that obtained by the NormFinder algorithm. Delta Ct method analysis found that miR162-3p and miR472 had lowest StdDev, followed by miR160 which suggested that miR162-3p and miR472 were two top-ranked, stable reference genes among all of the studied genes in sweet orange ( Figure 5A), which also confirmed the results from NormFinder analysis. For Ponkan (Figure 5B), miR403 followed by miR482b and miR162-3p, were found to be the most stably expressed genes with low StdDev (0.09, 1.04 and 1.04 respectively). RefFinder is a web-based tool integrating the above four computational statistical approaches to produce a comprehensive evaluation based on the geometric mean for each candidate gene. In the end, all untransformed raw data values were imported to RefFinder, then the software output ranks RefFinder is a web-based tool integrating the above four computational statistical approaches to produce a comprehensive evaluation based on the geometric mean for each candidate gene. In the end, all untransformed raw data values were imported to RefFinder, then the software output ranks reference genes according to the expression of reliability. When the stabilities from all the samples were combined, the output results agreed with the results of NormFinder software and Delta Ct methods. An overall comprehensive ranking output revealed that miR162-3p ( Figure 6A), followed by miR472, miR160 and miR166b (Geomean: 2.3, 2.78, 3.08, and 3.74, respectively) were the most stably expressed genes in canker-infected sweet orange, while U5 and GAPDH were the least stably expressed. The RefFinder algorithm showed that miR162-3p and miR403 constitute the best combination of two genes with comprehensive ranking values of 1.86 and 2.21 in all of the Ponkan samples ( Figure 6B). Again, U5 and GAPDH were ranked as among the least stably expressed genes, which had the highest Geomean. reference genes according to the expression of reliability. When the stabilities from all the samples were combined, the output results agreed with the results of NormFinder software and Delta Ct methods. An overall comprehensive ranking output revealed that miR162-3p ( Figure 6A), followed by miR472, miR160 and miR166b (Geomean: 2.3, 2.78, 3.08, and 3.74, respectively) were the most stably expressed genes in canker-infected sweet orange, while U5 and GAPDH were the least stably expressed. The RefFinder algorithm showed that miR162-3p and miR403 constitute the best combination of two genes with comprehensive ranking values of 1.86 and 2.21 in all of the Ponkan samples ( Figure 6B). Again, U5 and GAPDH were ranked as among the least stably expressed genes, which had the highest Geomean. The top three most stable genes were miR162-3p, miR160 and miR472 for the poly(A) extension method in sweet orange samples (Table 5). For the stem-loop method, st-miR396a, st-miR162-3p and st-miR166b were ranked as among the most stable genes, while st-miR3954 was the least stable genes under stress conditions. According to RefFinder, the comprehensive ranking poly(A) tail miRNA from the most stable to the least stable gene for Ponkan samples was as follows: miR160 > miR162-3p > miR472 > miR403 > miR428b > miR166b > miR396a > miR3954. While the comprehensive ranking of stem-loop miRNA for Ponkan samples was as follows: st-miR403 > st-miR162-3p > st-miR166b > The top three most stable genes were miR162-3p, miR160 and miR472 for the poly(A) extension method in sweet orange samples (Table 5). For the stem-loop method, st-miR396a, st-miR162-3p and st-miR166b were ranked as among the most stable genes, while st-miR3954 was the least stable genes under stress conditions. According to RefFinder, the comprehensive ranking poly(A) tail miRNA from the most stable to the least stable gene for Ponkan samples was as follows: miR160 > miR162-3p > miR472 > miR403 > miR428b > miR166b > miR396a > miR3954. While the comprehensive ranking of stem-loop miRNA for Ponkan samples was as follows: Taken together, UBC28, PP2A and EF1a were identified to be the most stable protein coding reference genes in all sweet orange samples. The protein coding reference genes for accurate transcript normalization in Ponkan were similar to those of sweet orange.

Validation of Candidate Reference Genes
To confirm the stability of the best ranked candidate reference genes, the expression pattern of CsLOB1 was examined in response to Xcc. The results clearly indicated that CsLOB1 was significantly upregulated after Xcc infection using the combination of miR162-3p/miR472 for normalization ( Figure 7A). Correspondingly, there was little difference in mock at each sampling day. However, when the poor reference gene GAPDH was used as an internal control, the expression of CsLOB1 showed a slight fluctuation in mock, but the expression differences between treatment and mock were significantly reduced ( Figure 7B).

Validation of Candidate Reference Genes
To confirm the stability of the best ranked candidate reference genes, the expression pattern of CsLOB1 was examined in response to Xcc. The results clearly indicated that CsLOB1 was significantly upregulated after Xcc infection using the combination of miR162-3p/miR472 for normalization ( Figure 7A). Correspondingly, there was little difference in mock at each sampling day. However, when the poor reference gene GAPDH was used as an internal control, the expression of CsLOB1 showed a slight fluctuation in mock, but the expression differences between treatment and mock were significantly reduced ( Figure 7B).

Discussion
A large body of evidence has shown that miRNA-mediated posttranscriptional regulation plays a fundamental role in plant responses to stressful factors, including pathogen infection, salt, drought, heat, cold, nutrient deficiency and metal toxicity [81][82][83][84][85][86]. A great deal of effort has been focused on the identification and functional analysis of plant miRNAs in response to the infection of pathogens, Figure 7. Relative expression of CsLOB1 in C. sinensis after Xcc infection or mock using the selected candidate reference genes: miR162-3p and miR160 (A) and GAPDH (B) for normalization. Means ± standard deviation (SE) were calculated from four discs per treatment (n = 4) per given day.
Northern blot, microarray and high-throughput sequencing have been used to analyze miRNA expression. However, RT-qPCR is still the most commonly used method. In RT-qPCR analysis, reference genes play a vital role in quantifying the expression level of target genes. Even a slight change in miRNA expression may affect targeting mRNAs for cleavage or repressing translation [96]. Improper reference genes may lead to inaccuracy or false results and incorrect conclusions [97,98]. Recently, there is an increased reporting on the evaluation of plant reference genes under different conditions [30,50,99,100]. However, limited information is available on miRNA reference genes in phytopathogenic research. Several studies showed miRNA in citrus plants [44,45,85], but there is no report on appropriate reference genes for normalizing miRNA expression in citrus canker-infected plants.
The present study identified appropriate reference genes for quantifying miRNA expression in two citrus species infected by Xcc. To achieve reliable results without bias, we selected 28 likely reference genes, including protein-coding RNA and non-coding RNA, and evaluated their expression levels using two different methods in the two citrus species, differing in resistance to Xcc. We then analyzed the performance of the reference genes by five methods and identified appropriate reference genes for each species. miR162-3p, miR160 and miR472 performed best as reference genes for sweet orange plants infected by Xcc, and miR160, miR162-3p and miR472 were the most stable genes for Ponkan attached by the same pathogen. To further confirm the reference genes, we used LOB as the target gene from sweet orange to verify their accuracy. Among the three reference genes, miR160 is known to mediate the interaction between auxin and cytokinin by suppressing the levels of ARF transcription factors, which plays a critical role in various aspects of plant development [101][102][103]. DICER-LIKE 1 (DCL1) is the main enzyme processing miRNA precursors into mature miRNAs and incorporated into the AGO1-RISC, which is considered indispensable for plant development [104]. Meanwhile, miR162 maintains a proper level of DCL1 transcripts by targeting DCL1 mRNA for a negative feedback regulation [104][105][106].
A single reference gene has been commonly used to normalize RT-qPCR results [107]. However, accumulating evidence has suggested that a single reference gene tends to show higher expression variability due to different biological factors [108]. Conversely, the use of multiple reference genes improved RT-qPCR statistical veracity [60,108]. In this study, geNorm software suggested that three reference genes are needed for a more accurate normalization for sweet orange infected by Xcc. As for Ponkan, two reference genes are found to be sufficient for accurately quantifying target gene expression levels ( Figure 3). Ponkan is a canker-resistant genotype and showed minor physiological changes and metabolic disorders than the susceptible sweet orange plant used in this study. All Ponkan samples remained relatively stable during Xcc infection processes at different infection stages. Consequently, Ponkan required fewer reference genes for normalizing qPCR results. Nevertheless, our results documented that the use of miR162-3p, miR160 and miR472 should accurately normalize the expression levels of the target gene in both sweet orange and Ponkan during Xcc infection.
GAPDH was considered to be a stable reference gene for qPCR normalization in different plants. GAPDH encodes key enzymes catalyzing the conversion of glyceraldehyde-3-phosphate to 1,3-biphosphoglycerate in the presence of NAD+ and inorganic phosphate in the glycolytic pathway [109][110][111]. Wu et al. (2014) suggested that GAPDH was the most reliable reference gene in navel orange (C. sinensis L. Osbeck) fruit during five different ripening stages [25]. In contrast, our results indicated that GAPDH was unstable, and was not an ideal reference gene for gene expression analyses, which was consistent with the reports of the others [112][113][114][115][116], that GAPDH was an unstably expressed gene in different varieties of oranges under different conditions. In another study, a total of 22 GAPDH genes were identified in wheat. Members of this family were found to notably respond to abiotic stresses [117]. In other words, some of these GAPDHs may participate in multiple functions except for metabolism roles in plants, such as involving abiotic stress resistance. Thus, caution is needed to use the GAPDH gene as a reference gene in RT-qPCR analysis.
CsLOB1 is transcription factor in the family of the Lateral Organ Boundaries Domain gene. All strains of Xcc can encode transcription activator-like (TAL) effectors which bind to effector binding elements in the promoter of CsLOB1 and induce the expression of disease susceptibility genes [63]. In this study, the identified reference genes (miR162-3p and miR472) showed highly reliable for normalizing CsLOB1. The expression profile of CsLOB1 was consistent with two suitable reference gene normalizations, which was also in agreement with the study of Hu [63]. When GAPDH was used as an internal control, CsLOB1 had no significant change in 24h, and mock showed fluctuating patterns of expression (Figure 7). These inaccurate results indicated that GAPDH would be unreliable for RT-qPCR analysis in the citrus bacterial canker response.
Analyzing miRNA expression profiles has inherent difficulties including (1) short nucleotide (18-24 nt); (2) heterogeneous GC content; (3) no sequence feature [e.g., poly(A)]; (4) target sequence interference; and (5) highly homologous within family [18]. Thus, poly(A) extension RT-qPCR and stem-loop RT-qPCR have been used for analyzing miRNA expression. In the stem-loop RT-qPCR, the reverse transcription of individual mature miRNAs uses a specific stem-loop primer. Due to the use of miRNA-specific primers, this approach is highly specific, thus effectively reducing background noise. This method, however, is time consuming and labor intensive. In the poly(A) extension RT-qPCR, miRNAs are first tailed by E. coli Poly(A) Polymerase and then reverse transcribed by oligo(dT) with a universal primer-binding sequence. This method is a universal reverse transcription, the product of the reaction contains cDNA from non-coding RNA (miRNA) and mRNA with a universal anchor primer. Compared to stem-loop RT, the Poly(A) extension method can analyze the transcript level of miRNA, miRNA target genes and other related genes with the same background. In the present study, we evaluated the two methods by comparing reference genes' stability in different citrus samples through two miRNA quantification methods. Our results showed that the average Ct value of poly-A tail qPCR was lower than the stem-loop RT method.
Moreover, the efficiency of PCR amplification for the stem-loop primer was lower than the poly-A tail, and some of them were even below 90%. Accordingly, the poly-A tail qPCR has higher amplification efficiency and better sensitivity than the stem-loop RT-qPCR. In forensic casework, poly-A tail extension also exhibited apparently more amplification than stem-loop RT for a body fluid identification of miR-451 and miR-205 [118]. In addition, different statistical algorithms revealed that the rank of the stability of the reference genes was in the following order: Poly(A) extension microRNA > stem-loop microRNA > protein coding RNA > non-coding RNA. In the present study, the poly(A) extension method achieved consistent results in two citrus species. Conversely, the top three stable stem-loop microRNAs were different in sweet orange and Ponkan, except st-miR162-3p. A study of Triticum dicoccoides showed that the results of stem-loop RT-qPCR were only partially consistent with data obtained from sequencing or microarray [119]. Similar results have been reported in Triticum aestivum L., where five miRNA genes in different wheat tissues were analyzed by two miRNA qPCR approaches, as well as deep sequencing. The deep sequencing results had higher correlation coefficients with results obtained from poly(A) RT-qPCR, but not with the results derived from stem-loop RT-qPCR [50].

Conclusions
The selection of reference genes is critically important for the accurate quantification of target gene expression in RT-qPCR analysis. To identify appropriate miRNA reference genes for normalizing the level of miRNA expression in C. sinensis and C. reticulata infected with a citrus canker pathogen, five algorithms: Delta Ct, geNorm, NormFinder, BestKeeper and RefFinder were used for screening reference genes. Poly(A) extension RT-qPCR and stem-loop RT-qPCR were performed to determine the most appropriate method for detecting expression patterns of miRNA. Our results showed that poly(A) RT-qPCR is a more reliable method than stem-loop RT-qPCR, and miR162-3p and miR472 genes in poly(A) RT-qPCR are the most stable internal control genes for miRNA profiling analysis of citrus infected by Xcc.
Author Contributions: Conceptualization, D.P. and J.C.; methodology, D.P. and J.C.; formal analysis, W.C.; investigation, S.L. and Y.Y.; resources, S.X. and G.C.; writing-original draft preparation, S.L. and J.C.; writing-review and editing, S.L. and J.C.; project administration, W.S. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors would like to thank Agricultural "Five New" Project, demonstration and promotion of key technologies for promoting ecological quality and high efficiency production (Fujian Development and Reform Commission, 2018) for supporting this study.