Selection and validation of reliable reference genes for gene expression studies from Monilinia vaccinii-corymbosi infected wild blueberry phenotypes

Monilinia blight disease caused by Monilinia vaccinii-corymbosi (Reade) Honey (M.vc) causes severe damage and economic losses in wild blueberry growing regions. Molecular mechanisms regulating defence responses of wild blueberry phenotypes towards this causal fungus are not yet fully known. A reliable quantification of gene expression using quantitative real time PCR (qPCR) is fundamental for measuring changes in target gene expression. A crucial aspect of accurate normalisation is the choice of appropriate reference genes. This study evaluated the expression stability of seven candidate reference genes (GAPDH, UBC9, UBC28, TIP41, CaCSa, PPR and RH8) in floral tissues of diploid and tetraploid wild blueberry phenotypes challenged with M.vc. The expression stability was calculated using five algorithms: geNorm, NormFinder, BestKeeper, deltaCt and RefFinder. The results indicated that UBC9 and GAPDH were the most stable reference genes, while RH8 and PPR were the least stable ones. To further validate the suitability of the analyzed reference genes, the expression level of a pathogenesis related protein gene (i.e., PR3) was analysed for both phenotypes at four time points of infection. Our results may be beneficial for future studies involving the quantification of relative gene expression levels in wild blueberry species.

Scientific RepoRtS | (2020) 10:11688 | https://doi.org/10.1038/s41598-020-68597-9 www.nature.com/scientificreports/ has become an effective approach to examine and validate the changes in gene expression owing to its accuracy, specificity and sensitivity [8][9][10][11] . However, the reliability of the technique is largely influenced on the normalization strategy such as the usage of one or more stable reference genes. Ideally, the reference genes used for normalisation should have a uniform expression regardless of the experimental conditions 11,12 . Hence, the selection and validation of stable reference genes for each experimental condition is a prerequisite for performing qRT-PCR studies [13][14][15] .
In the present study, seven reference genes (GAPDH, CaCSa, TIP41, UBC28, PPR, UBC9 and RH8) were selected as candidate reference genes based on previous reports in Vaccinium spp 16,17 . Their expression stabilities in floral tissues of diploid and tetraploid wild blueberry phenotypes challenged with Monilinia vaccinii-corymbosi was evaluated at four time-points after infection. Five different statistical software programs (geNorm 8 , NormFinder 18 , BestKeeper 19 , delta Ct 20 and RefFinder 21 ) were used to analyse the stability of the candidate reference genes and to select the most appropriate ones. This study will lay a foundation for future gene expression research in wild blueberry.

Results
Selection of reference genes, amplification efficiency and specificity. A total of seven candidate reference genes were selected to identify suitable RGs for gene expression studies using qPCR in wild blueberry. Additionally, PR3 gene was used to validate the accuracy of identified RGs in wild blueberry phenotypes infected with M. vaccinii-corymbosi. The specificity of the analysed primer pairs was confirmed via detecting single fragment of the expected size on 2% agarose gel electrophoresis and a single peak with no signals on the negative controls in the melt curve analysis ( Table 1, Supplementary Fig S1 & S2). All the tested RGs had efficiency (E %) values ranged from 95 to 105%, with regression coefficient (R 2 ) varying from 0.994 to 0.999 (Table 1, Supplementary Fig S3). The results showed that all the primer pairs were suitable for RT-qPCR analysis.
Expression profiling of candidate reference genes. The raw quantification cycle (C q ) values were used to quantify the expression levels of candidate reference genes where lower Cq values mean higher expression levels. Cq values for each of the seven candidate reference genes in V. myrtilloides and V. a. f. nigrum are listed in Supplementary Table S1 (there were no Cq values in the negative controls), and a box and whiskers plot were used to describe the raw Cq value distribution (Fig. 1) Stability ranking of the candidate reference genes. Expression stabilities of the seven candidate reference genes were determined using geNorm, NormFinder, ΔCq, and BestKeeper and their overall stabilities were ranked by RefFinder across all the timepoints and phenotypes. geNorm analysis. The expression stability rankings based on the M-values for the seven candidate reference genes on wild blueberry phenotypes is displayed ( Table 2). The gene with the lowest M-value (cut-off 1.5) was the most stable reference gene in terms of gene expression and vice versa. The M-values for the tested genes in all samples and groups were lower than the default limit of 1.5. Genorm analysis for V. myrtilloides illustrated that CaCSa (0.385) and GAPDH (0.386) were the more stable genes and PPR (0.888) exhibited the least stability. On contrary, UBC9 (0.239) and GAPDH (0.239) exhibited high expression stability in V. a. f. nigrum and RH8 (0.432) were the least stable one (Table 2). Among the total samples, GAPDH and UBC9 were the most stable genes with M values of 0.487 and 0.512, whereas RH8 exhibited least stability (0.701). Finally, the pairwise variation (V n /V n+1 ) for both the phenotypes and entire group resulted in V 2/3 < 0.15 ( Fig. 2) which indicated that two reference genes were sufficient for accurate normalisation of RT-qPCR data.
NormFinder analysis. The expression stability of the seven candidate genes was analyzed using NormFinder, which is an excel based mathematical tool that measures gene expression stability by comparing the variation within and between user-defined sample groups. NormFinder ranks the control genes on the basis of their stability value (SV), lower value indicates higher gene expression stability and vice versa. The NormFinder algorithm results agreed with GeNorm analysis ( Table 2). NormFinder selected CaCSa (SV = 0.043) and UBC9 (SV = 0.047) as the most stable genes for V. myrtilloides, whereas, GAPDH (0.014) and PPR (0.047) for V. a. f. nigrum. For overall analysis, UBC9 (0.080) and CaCSa (0.104) exhibited the most stable genes and RH8 the least stable one (0.195).
BestKeeper analysis. BestKeeper ranks the stabilities of the candidate reference genes based on their standard deviation (SD) and the coefficient of variation (CV). Genes with SD > 1 were considered unacceptable reference genes. Based on the results from the BestKeeper analysis for V. myrtilloides, all genes except RH8 were calculated to have an SD value lower than 1 ( Table 2). The rankings by BestKeeper analysis for V. myrtilloides showed that the most stable reference genes were CaCSa (0.62) followed by UBC28 (0.75) and UBC9 (0.75). For V. a. f. nigrum UBC28 (0.42) and GAPDH (0.43) were observed as the most stable ones ( Table 2). In the total dataset, UBC28 and CaCSa were the most stable genes, with SD values of 0.59 and 0.62 respectively (p < 0.001). For all the three groups RH8 exhibited the least stability. RefFinder analysis. A comprehensive ranking was performed to confirm the stability ranking of the seven candidate reference genes. In the comparative analysis of reference genes from V. myrtilloides and V. a. f. nigrum, UBC9 and UBC28 were ranked as the most stable RGs (Fig. 3). However, due to the differences in ploidy as well as defense response level, a comprehensive ranking order of reference genes was generated for each phenotype. The expression of CaCSa, UBC9 and UBC28 were found to be the most stable for gene expression normalization in V. myrtilloides. While, GAPDH, UBC9 and PPR ranked as the best for V. a. f. nigrum. By contrast, RH8 and TIP41 ranked as the least stable genes for both the phenotypes.

Validation of the best and least ranked reference genes.
To validate the effectiveness of the selected reference genes, the expression pattern of PR3 (MK292725) was analyzed in wild blueberry-Monilinia pathosystem. The relative transcript abundances of PR3 gene was normalized to two most stable genes (UBC9 and www.nature.com/scientificreports/ GAPDH) for V. a. f. nigrum and (UBC9 and UBC28) for V. myrtilloides as resulted from the comprehensive analysis.
In V. myrtilloides, the tolerant phenotype, data normalizations using the two most stable reference genes (UBC9 and UBC28) resulted in consistent PR3 expression pattern with gradual increase in expression over time points (Fig. 4A). However, in V. a. f. nigrum, the susceptible phenotype, the expression was high at day 0, then reduced the expression and not significantly elicited compared to control.
The expression data was also normalized using the least stable gene (RH8) reported for the two phenotypes (Fig. 4B). For both V. a. f. nigrum and V. myrtilloides no remarkable expression observed for all the analysed timepoints. These results indicated that the least stable gene, RH8 failed to standardize the expression data effectively. Our results confirm that using different reference genes for normalisation causes great differences among the expression patterns.

Discussion
RT-qPCR has emerged as a powerful tool to study transcript abundance of a specific gene in distinct biological samples owing to its precision, accuracy and sensitivity 11 . However, accurate normalization of gene expression remains a major criterion during qPCR analysis, as various steps during qPCR analysis can introduce variations arising from RNA extractions, cDNA synthesis, PCR procedure and sample loading 13 . Using a stably expressed reference gene is a prerequisite to obtain accurate interpretation of the transcript abundance results. Ideally, a reference gene should have constant expression in samples irrespective of the experimental conditions, developmental stages or species 22,23 . Several studies have reported the observation of variable expression of traditional reference genes in different plant species, as they have differential expression under different experimental conditions [24][25][26] . For example, Czechowski et al. 27 reported that Arabidopsis have demonstrated variations in www.nature.com/scientificreports/ reference gene stability under different experimental conditions. Thus, the most appropriate reference genes should be properly evaluated and confirmed in all biological samples. Many reliable reference genes have been determined in Vaccinium spp 16,17,28 however, this study is the first to assess reference genes in wild blueberry-monilinia pathosystem. Also, the uniqueness of this study is that it was performed under field conditions. Field conditions differ from controlled environmental conditions, as plants are exposed to adverse environmental factors and multiple stresses. As stated by Samarth and Jameson 29 , the selection of appropriate reference genes from field conditions is more complex than controlled environmental studies. Tashiro et al. 30 also pointed out the necessity of reference gene validation in studies involving non-model plant species from heterologous plant population. In this study, we assessed seven RGs for their use as internal www.nature.com/scientificreports/ controls in gene expression studies of the wild blueberry phenotypes upon monilinia blight infection under field conditions. In the present study, four algorithms (geNorm 8 , NormFinder 18 , BestKeeper 19 and delta CT 20 ) was used to evaluate the best suited reference gene in wild blueberry phenotypes, V. myrtilloides and V. a. f. nigrum, for studying gene expression pattern during Monilinia blight infection. Several studies indicated the use of multiple statistical algorithms which not only minimize the errors associated with reference gene selection but ensures a more reliable evaluation [31][32][33][34] . This is probably due to the differences in algorithm programs exhibited by each method 23,33 . Moreover, variability in the stability ranking of reference genes based on algorithms has been reported in several studies 35,36 . Interestingly, the ranking order generated by using four programs roughly the same for both the phenotypes, with the first three most stable genes differing only on the ranking order. In addition, the pairwise variation determined using geNorm could be used as an indicator for estimating the optimal number of reference genes for normalization. In this study, the pairwise variation values for both V. myrtilloides and V. a. f. nigrum was V 2/3 < 0.15, indicates that two reference genes were sufficient for gene expression normalization.
Even though the ranking order of the candidate reference genes may have differed, all the statistical programs showed UBC9 as the most stable reference gene in the analysed phenotypes having different ploidy level. In this study, GAPDH and UBC9 were observed as most stable reference genes for V. angustifolium. f. nigrum and UBC9, www.nature.com/scientificreports/ UBC28 and CaCSa for V. myrtilloides. GAPDH was used for normalization in studies with V. myrtillus-Botrytis pathosystem 37 . Vashisth et al. 16 and Die and Rowland 17 ranked reference genes on vegetative and reproductive organs of rabbiteye and southern highbush blueberry. Vashisth et al. 16 reported CaCSa, RH8, and UBC28 as the most stably expressed gene in southern highbush blueberry across multiple organs analysed. Die and Rowland 17 also reported RH8, CaCSa, PPR, GAPDH and UBC9 (M = 0.483 and CV = 0.210) as the most stably expressed reference genes for floral bud tissues evaluated in highbush blueberry. However, in our study, RH8 and TIP41 were found to be the least stable reference gene for all the assessed phenotypes. Several studies have reported variable expression levels in reference genes among closely related species 38 . Moreover, several studies demonstrated variations in expression profiles of reference genes in different pathosystems 33,39 . PR3 was used as a target gene to validate the credibility of the selected reference genes. PR3 which include chitinases of classes Ia, Ib, II, IV, VI, and VII are important weaponry of plants against pathogens 40 . According to Thomma et al. 41 , PR3 which belong to the pathogenesis-related (PR) protein family play an important role in plant defense response to necrotrophic pathogens. Several studies reported the up-regulation of PR genes against many phytopathogenic fungi [42][43][44] . When comparing the most stable and least stable reference genes, the expression of target gene was consistent and upregulated in V. myrtilloides (UBC9 and UBC28) and V. a. f. nigrum (UBC9 and GAPDH) even though differences in expression observed between the phenotypes. However, no response observed in both phenotypes when analysed using the least stable gene (RH8). Our results agree with the findings of Cardot et al. 45 , where elicitation of chitinase genes observed in tolerant than in susceptible varieties.
In conclusion, this is the first study in which a set of candidate reference genes was analysed in terms of their expression stability in wild blueberry phenotypes infected with Monilinia vaccinii-corymbosi. Five different statistical algorithms showed slight differences in the final ranking of reference gene, however by combining and analysing the data together, we demonstrated that UBC9 is the most stably expressed transcript in wild blueberry phenotypes regardless of ploidy level. www.nature.com/scientificreports/ Methods plant material, M.vc inoculations and experimental design. The wild blueberry diploid (V. myrtilloides) and tetraploid (V. angustifolium f. nigrum) phenotypes were selected from a commercial wild blueberry field, NS, Canada. Three biological replicates were selected for each phenotype and each replicate was separated into two, 0.5 × 1 m sample areas. One day before inoculation, one sample area within each replicate was sprayed with the fungicide Proline (a.i. prothioconazole) at a rate of 315 ml product·ha -1 using a CO 2 powered, Bell spray Inc. hand-held research sprayer with 2 m boom with 4 Tee Jet Visiflow 8003VS nozzles at a pressure of 220 kPa to serve as control plots. Monilinia vaccinii-corymbosi inoculum was prepared from four-week old M.vc cultures isolated from monilinia blighted shoots and mummy berries. Floral buds at F3 stage (floral bud scale separation and appearance of new growth) was tagged for each phenotype and inoculum (2 × 10 5 ascospores·mL −1 ) was sprayed at all angles until runoff. Each phenotype sample area was immediately covered with 2 mm plastic film and row cover to provide incubating conditions (100% RH), which created conditions required for Monilinia infection 46 . After 72 h, the plastic film and row cover were removed and floral bud tissue from 15 random stems in each plot (control and inoculated) was harvested for RNA extraction and immediately flash frozen in liquid and stored at − 80 °C until further use. Floral bud tissues were collected as day 0 (before inoculation), 3, 6 and 10 days after inoculation.
total RnA isolation and cDnA synthesis. Total  Quantitative real-time pcR (qpcR). qPCR assay was performed using a CFX Connect Real-time Detection System (Biorad, US) to analyze the specific expression of reference/target gene. Each PCR reaction mixture (10 µl) contained 2 µl of diluted cDNA (20-fold dilution), 5 µl SsoAdvanced SYBR Green Supermix (Biorad), and 1 µl (10 nM) of each forward and reverse primer. The amplification conditions were as follows: an initial denaturation at 95 °C for 180 s, followed by 40 cycles at 95 °C for 10 s, 60 °C for 20 s. Each run was completed with a melting curve analysis (65-95 °C with at increments of 0.5 °C) to verify the specificity of the amplification. The cycling conditions were based on the method described by Petriccione et al. 33 . A no-template control (NTC) was included for each gene assay to confirm the absence of non-specific products 47 .
Determination of reference gene expression stability. The Cq value (quantification value) of each reference gene under the four different time points for both V. myrtilloides and V.a.f nigrum was recorded using the qPCR system. Four widely used software: geNorm 8 , NormFinder 18 , BestKeeper 19 , and ΔCq 20 method was used to rank the expression stability of the reference genes. Finally, we used RefFinder 21 , a web-based userfriendly comprehensive tool, which integrates all four algorithms providing an overall ranking of the used genes. The GeNorm algorithm, which is a module of qbase + software package (Biogazelle), was used to evaluate the candidate reference genes based on their expression stability values (M-values) and pairwise variations (Vn/ Vn + 1) 48 . The default set value was 1.5; gene with the lowest M-value was the most stably expressed one. The computed pairwise variation (Vn/n + 1), was used to determine the optimal number of reference genes required for normalisation of the data. A (Vn/n + 1) value < 0.15 indicated the appropriate number of reference genes required for analysis 8 .
For Normfinder, the raw Cq values were converted into relative quantities (RQ) using the formula RQ = 2 (Cq min − Cq sample) , where Cq min is the lowest Cq value across the sample pool. Normfinder evaluated the expression stability of candidate reference genes at inter-group and intra-group levels. Ideally, the two genes with the lowest stability values were the most appropriate genes to be used for normalisation 18 .
BestKeeper was performed using the original Microsoft Excel-based formulas 19 . It calculates the standard deviation of the Cq value between the whole data set, and the gene with the lowest standard deviation (SD) is proposed as most suitable. The comparative ΔCq method manually compares relative expression of pairs of genes within each sample.