Identification of suitable reference genes for expression studies in pomegranate under different biotic and abiotic stress conditions

Pomegranate (Punica granatum L.) is an important economic fruit crop, facing many biotic and abiotic challenges during cultivation. Several research programs are in progress to understand both biotic and abiotic stress factors and mitigate these challenges using gene expression studies based on the qPCR approach. However, research publications are not available yet to select the standard reference gene for normalizing target gene expression values in pomegranate. The most suitable candidate reference gene is required to ensure precise and reliable results for qPCR analysis. Eight candidate reference genes' stability was evaluated under different stress conditions using different algorithms such as ∆Ct, geNorm, BestKeeper, NormFinder, and RefFinder. The various algorithms revealed that EFA1 and 18S rRNA were common and most stable reference genes (RGs) under abiotic and wilt stress. Whereas comprehensive ranking by RefFinder showed GAPDH and CYPF were the most stable RGs under combined biotic (pooled samples of all biotic stress) and bacterial blight samples. For normalizing target gene expression under wilt, nematode, bacterial blight, and abiotic stress conditions both GAPDH and CYPFreference genes are adequate for qPCR. The above data provide comprehensive details for the selection of a candidate reference gene in various stresses in pomegranate


Introduction
Pomegranate (Punica granatum L.) is an important fruit crop of the world's subtropical and tropical regions. It has been used in traditional medicine to cure many diseases across different cultures and civilizations [1]. At the global level, India is the leading producer of pomegranate, with a yearly production of 2,442 thousand tonnes grown in 209 thousand hectares. However, pomegranate production has been hampered by various biotic and abiotic factors leading to severe yield loss in quality and quantity. Major biotic and abiotic factors include bacterial blight caused by Xanthomonas citri pv. punicae [2], wilt caused by Ceratocystis fimbriata [3], drought [4], and cold [5,6], respectively. Understanding the plant's underlying mechanisms to these stresses with detailed molecular studies through gene expression analyses is essential.
The gene expression analysis requires the most sensitive, accurate, and reproducible measurements for specific mRNA sequences. Real-time qPCR is emerged as the most sensitive yet powerful technique to study mRNA abundances and quantifying gene expression [7,8]. The gene expression studies required good reference genes for consistent and reproducible results. These genes are essential for maintaining the cell's basic functioning and are stably expressed regardless of developmental stage or environmental conditions. Presently, the normalization of quantitative gene expression is carried out by several traditional reference genes, including 18S [9], Actin [10], GAPDH [11], and CYP [11]. However, in pomegranate reference genes such as 18S rRNA, Ribosomal protein S, named PgRPSΙΙ, and GAPDH has been used for gene normalization in various gene expression studies [12][13][14]. However, this may misinterpret the gene expression because all reference genes may not constantly express in different tissues under diverse environmental conditions. Both abiotic and biotic stresses are the major factors affecting gene expression, and wide variation exists in the expression of RGs under stress conditions [15]. Therefore, it is essential to identify the most stable reference genes to measure differentially expressed target genes.
Recently Pomegranate chloroplast and whole-genome have been sequenced, although a significant portion of the genome has been annotated [16,17], most of the genes remain functionally uncharacterized. In spite of being widely cultivated in India and other parts of the world, investigation on the molecular basis of biotic and abiotic stress is limited. It is a prerequisite to functionally characterize stress-associated genes by identifying the stable reference genes for different stress. Thus, the present study identifies the candidate reference gene for biotic stress such as bacterial blight, wilt, nematodes (Fig. 1), and abiotic stress such as drought and cold stress in pomegranate.

Plant material
Three varieties of pomegranate viz., Bhagwa, Nana, and Daru were grown using stem cuttings approximately 12 cm long in plastic pots filled with vermicompost, soil, and sand (1:1:3) and maintained at 29 ± 2 °C at a relative humidity of 60-80% in the greenhouse. One-year-old plants were used to carry out experiments for biotic and abiotic stress. Bhagwa, Nana, and Dharu were used for the bacterial blight experiment and only Bhagwa was used in other experiments. Experimental units were arranged in a completely randomized design with five replications.

Bacterial blight infection
Pure and a single colony of XAP (Accession no. KX702398.1) culture was grown on NGA medium (0.5% proteus peptone, 0.3% beef extract, and 0.25% sucrose) for 48 h at 2 ± 0.5 °C. The bacterial concentration was set for 0.3 OD at A 600 nm using Bio-Spectrometer (Eppendorf AG, Hamburg, Germany). The XAP was spray inoculated on pomegranate plants Bhagwa, Nana, and Dharurespectively, as per described protocol [18]. Leaves samples were collected 10 days post pathogen inoculation, for reference gene identification.

Nematode infection
The root-knot nematode M. incognita, was used for the experiment, according to [19]. Nematode population was maintained on the tomato plants (Lycopersicon esculentum) and Non -sterile egg masses were collected in a 2 ml tube containing 500 ml of sterile distilled water. The suspension was centrifuged at 1000×g for 5 min. The supernatant collected was re-suspended in 5 ml sterile distilled water and ready to use. Pomegranate plant (cv. Bhagwa) was root inoculated with1000 J2 stage M. incognita and control plants were inoculated with water and then maintained for approximately 25 to 30 days under greenhouse conditions. For gene expression studies, 30 days old root and leaf samples were collected.

Wilt infection
Ceratocystis fimbriata (UHS-CF5, Acc. No KU877189) grown on Potato dextrose broth at 25 ± 0.5 °C, with a photoperiod of 12 h for 7 days. Before inoculation of the pathogen, plants were wounded on the xylem using cork bore to facilitate the pathogen entry. The pathogen was then artificially inoculated by soil drenching using a 30 ml inoculum with spore concentrations of 4 × 10 8 spores/ml. Plants treated with water were maintained as control. For gene expression analysis, 45 days old root and leaf samples were collected.

Cold stress
One-year-old Bhagwa plants were subjected to cold stress. Cold treatment was performed by keeping plants at 4 °C for 48 h in the cold growth chamber. Leaves were sampled at 42 h post-stress imposition of the stress and immediately flash feezed in liquid nitrogen and stored at -80 ℃ until use.

Water stress
In the greenhouse, six-month-old pomegranate plants of Bhagwavariety were subjected to water-limited conditions. All plants were maintained at field capacity before the imposition of water stress. Control plants were maintained at field capacity throughout the experiment according to the protocol described by [20], whereas water-stressed plants were subjected to a water-limited condition for 12 days by withholding until the plants showed a wilting symptom. Leaves were sampled at 12 days post-stress imposition and immediately frozen in liquid nitrogen and stored at − 80 ℃ until use.

RNA isolation and first-strand cDNA synthesis
Total RNA was isolated using RNeasy plant mini kit (Sigma, USA) according to pavan et al. [21], for root sample incubation period was extended for an additional 15 min. Further, post RNA isolation samples were treated with RNase-free DNase I (Thermo Fisher, MA, USA) to remove DNA contamination. The purity of RNA was measured by taking the observance at 260/280 nm using NanoDrop Spectrophotometer (ND-1000 Thermo Fisher, MA, USA). cDNA synthesis was performed taking 1 μg of total RNA in a 20 μl reaction using maxima first strand cDNA Synthesis Kit (ThermoFisher, K1672, USA) following manual instruction. Constructed cDNA sample was diluted 10 times with nuclease-free water and stored at -20 °C for qPCR analysis.

RT-qPCR
All the quantitative reverse transcription-polymerase chain reaction (RT-qPCR) reactions were performed using SYBR Green I technology, on StepOne plus a real-time PCR machine (Applied Biosystems, USA). A cocktail mixture for each PCR reaction was constituted of 1X Power UP SYBR Green master mix (Thermo Fisher Scientific), 0.5 µM primers (forward and reverse primers), and 2 µL of tenfold diluted cDNA to the reaction volume of 10 µL. The following PCR condition was used: initial denaturation for 30 s at 95 °C, followed by 40 cycles of 95 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s. All the qPCR analyses were carried out in three biological and two technical replicates and the total mean was considered for analysis. To determine gene-specific PCR efficiency (E), tenfold serial dilutions (1-1000) of bulked cDNA were used for generating a standard curve for each gene and each stress condition. Amplification efficiency (E) of each primer was determined by linear regression model according to the equation E (%) = (10 − 1/ slope − 1) × 100% [22]. The efficiency of the PCR products was calculated by measuring the CT to a specific threshold of each gene [23]. The correlation coefficients (R 2 ) values were acquired from the standard curve.

Data mining and statistical analysis to determine the expression stability
Standard curves were generated from a fivefold series dilution sample from pooled cDNA as a template for each primer, in three replicates. The melting curves of each primer pairs are shown in Fig. 1. The stabilities of selected reference genes were evaluated using Microsoft Excel-based software tools, NormFinder [24], BestKeeper [25], geNorm [26]. BestKeeper program calculates the stability of the expressed gene using raw data i.e., ct values, the geNorm software program analysis based on the stability of the expressed values (M), and ranks the gene order. Lower the M value (stability value), higher is the stability of any gene expression, M value less than 1.5 considered as ideal reference genes in all stress conditions. NormFinder calculates the stability of expressed genes based on intra-and intergroup variations and gives the gene rank order by combining these values, lower the rank higher is the stability of any expressed gene. RefFinder performed comprehensive analysis using the data from geNorm (M values), NormFinder (Stability values), BestKeeper (CV and SD), and ΔCt values.

Primer specificity and PCR amplification efficiency for each candidate reference gene
To identify appropriate reference genes for pomegranate under biotic and abiotic stress, eight candidate reference genes such as EFA1, CYPF, UBQ, βTUB, αTUB, GAPDH, 18S, and PPCT were selected based on the previous studies on various plant species. The selected reference gene sequences were obtained from the National Centre for Biotechnology Information (NCBI, USA). The products of these eight reference genes were associated with a wide variety of biological functions of plants. To check the primer specificity for these candidate reference genes, we performed melting curve analysis using qPCR and agarose gel electrophoresis. The PCR product of each primer pair exhibited a single peak in the melting curve and a single band with the expected size after agarose gel electrophoresis (Figs. 2 and S1). The amplification efficiency for the eight candidate reference genes ranged from 90.65% (GAPDH) to 103.04% (PPCT), and correlation coefficient (R2) values ranged from0.981 (GAPDH) to 0.999 (EF1-α and 18S rRNA) (Table S1).

Delta-Ct (ΔCt)
The stability of reference gene using delta-Ct analysis, which was represented by the standard deviation (SD). The ranking of reference genes based on the ΔCt method is presented in Table 1. According to the ΔCt analysis, EFA1 was the most stable reference genes with low SD under all the stress conditions, except under wilt infection. 18S rRNA showed stable expression under wilt infection with a standard deviation of 1.27.

Best keeper
Data on reference gene expression was reanalyzed using BestKeeper to identify stable reference genes [25]. The analysis revealed that GAPDH was the most stable gene in both biotic stress and bacterial blight infection with a standard deviation (S.D) of 0.22 and 0.73 respectively ( Table 2). Pairwise correlation between the reference gene was performed by the best keeper algorithm to calculate the stability of individual genes. Alpha tubulin showed the most stable expression under abiotic stress with the least S.D of 0.05.

NormFinder
NormFinder is excel based application, which calculates expression stability for individual reference genes based on variance analysis. A higher stability value of the reference gene indicates that they are unstable, whereas a lower stable value indicates greater gene stability. According to Norm Finder the stability analysis (Table 3), EFA1 occupied the top position with a stability value of 0.07 and 0.12 in abiotic and bacterial blight. CYPF, 18S, and βTUB were the two most stable genes under biotic, wilt, and nematode respectively. geNorm geNorm algorithm calculated the average gene stability (M value) of eight candidate reference genes. Similar to Nor-mFinder, the lower the M value, the higher the gene stability (Fig. 2). GAPDH and EFA1 were most stable under biotic stress, 18S and CYPF stable under bacterial blight, BTUB and CYPF stable under wilt infection, EFA1 and ATUB stable under nematode infection, BTUB and EFA1stable under abiotic stress (Table 4). To determine the optimal number of reference genes for qPCR normalization for different stress in pomegranate, the average pairwise variation was calculated between two consecutive normalization factors (Vn/Vn + 1, where n represents the reference gene number) across individual stress. According to Vandesompele and colleagues, recommend Vn/Vn + 1 cut-off should be less than 0.15, although this threshold should not be regarded as too stringent. The results exhibited that the Vn/Vn + 1 values for different stress were below the cut-off value of  (Fig. 3), except for biotic stress, which indicated that two reference genes would be sufficient for normalization of qPCR analysis in pomegranate.

Comprehensive ranking by RefFinder
RefFinder is a web-based interface, which calculates the geometric mean for each reference gene based on ranks calculated by other approaches. According to RefFinder  . 3 The optimal number of reference genes for accurate normalization was calculated by geNorm during different biotic and abiotic stress. Pairwise variation (Vn/Vn + 1) analysis of eight candidate reference genes analyzed in five different conditions expression stability of all reference genes were sorted according to the geometric mean (Table 5). Reference genes with the lowest geometric mean are considered the most stable. GAPDH ranked the most stable gene for the biotic stress sample group with a stability value of 1.86. 18S rRNA for wilt infected samples with a stability value of 1.46. Under bacterial blight infection, CYPF showed the least stability value (1.86). For abiotic stress, EFA1 was found to be highly stable with a stability of 1.18.

Discussion
In the post-genomic era, various omics approaches have empowered us to understand the molecular mechanisms underlying the defence response of plants subjected to various stress. qPCR is the most widely used technique to study the molecular mechanisms of plant stress responses by analysing gene expression. However, gene expression studies are affected by several factors such as variation in RNA integrity from sample to sample, the variation in reverse transcriptase efficiency, and the quantity of cDNA template in each PCR reaction. Hence, normalization of target gene expression level against a stably expressed reference gene can eliminate all kinds of variation across samples [25]. The accuracy of the qPCR results strongly depends on selecting one or more reference genes that are stably expressed across different tissue and experimental conditions. Thus, in the present study, eight commonly used reference genes (EFA1, CYPF, UBQ, βTUB, αTUB, GAPDH, 18S, and PPCT) were used for evaluating their expression stability under both biotic (Bacterial blight, wilt, and nematode) and abiotic (cold and moisture stress) stress conditions in pomegranate. The selection of reference genes based on the literature that these genes are primarily involved in the growth and development of the plant and they are typically expressed at constant levels under different environmental conditions in different tissues.
The primer pairs specificity of the RG's was confirmed by melt curve by qPCR and agarose gel electrophoresis analysis (Figs. 1 and S1). The amplification efficiency of all RGs exhibited between 92.1 to 103 and r2 ≥ 0.99, considered acceptable [27,28]. Expression stability of all RGs was estimated using different statistical algorithms such as ∆ Ct, BestKeeper, NormFinder, and geNorm. These methods use distinct statistical algorithm procedures and have their advantages and drawbacks [29]. However, geNorm has been considered the best algorithm for identifying stable RG and it is highly sensitive for co-regulation compared to norm finder [30,31]. According to geNorm, EFA1 was found to be a stable gene under biotic stress and abiotic stress. Similarly, in ∆Ct method, EFA1 exhibited high stability across different stress conditions than other RGs in pomegranate. EFA1 was considered stable RG in Eucalyptus globulus during cold acclimation, and it has a relatively close taxonomic relationship with pomegranate based on comparative genomic analysis [32]. EFA1 also reported being a stable RG in different crops such as Cucumber [33], Tobacco under abiotic stress [34], Maize for abiotic and hormonal treatments [35], Perlmillet under abiotic stress [36]. CYP2 showed stability of 0.03 under H 2 O 2 stress using normfinder [37]. UBQ showed high stability with M value 0.16 and 0.06 under drought and cold stress respectively in white clove [38]. Ref-Finder provides consensus expression stability values for RGs by combining other statistical algorithm methods such as ∆Ct, BestKeeper, NormFinder, and geNorm. Comprehensive ranking by RefFinder indicated that EF-1α was the most suitable reference gene for nematode infection and abiotic stress, CYPF for bacterial blight, 18S for wilt infection, and GAPDH for biotic stress. In virus infection of tomato, GAPDH (stability 22.4) and ubiquitin (stability 20) were the highly stably expressed RG, when analysed using bestkeeper [39]. Interestingly, 18S was the most stable RG ranked by all other methods under wilt infection, except geNorm. 18S rRNA was also reported to be stable under different stress conditions in crops such as eggplant [40] in tea plants (Camellia sinensis) exposed to Zn stress [41]. However, few studies have reported that due to a higher abundance of 18S rRNA transcripts than the target gene in the tissue, it is difficult to subtract the baseline value in qRT-PCR accurately [42][43][44]. Though the reference gene's stable expression is important in qPCR normalization, 18S rRNA cannot be considered a bad gene for normalization [26]. It can be used as the reference gene for those target genes with high expression under certain environmental conditions [42]. Using a single RG for normalization will affect the accuracy of calculating the target gene's expression [45]. Hence, using two or more RGs or using a single, highly stable RG for normalization would increase the accuracy of knowing the target gene's exact expression in a given sample. The geNorm algorithm determined the optimal number of reference genes required for normalization in pomegranate under different stress conditions. With a threshold of 0.15, two reference genes were sufficient for normalization under wilt, bacterial blight, nematode, and abiotic stress. More than two RGs are required to normalize under biotic stress (pooled samples of all biotic stress). In biotic stress, the sample's complexity may result in higher variability in the expression of RGs. However, using more RGs would help to reach optimum results.