Evaluation of suitable reference genes for gene expression analysis in the northern root-knot nematode, Meloidogyne hapla

The northern root-knot nematode (Meloidogyne hapla) is a critical pathogen with a wide host range. Quantitative real-time polymerase chain reaction (qRT-PCR) has been used to elucidate gene expression and function of M. hapla. Suitable reference genes are required to ensure accurate results of qRT-PCR for normalising gene expression. Eleven candidate reference genes of M. hapla were selected to evaluate gene expression stability under different conditions. The stability of candidate reference genes was ranked using RefFinder analysis, and the optimal number of reference genes was recommended with geNorm. Notably, the most stable reference genes were SDHA, Mdh, and RpS6 for all samples; SDHA and RpS6 were particularly stable during development stage treatments, whereas Mdh and RpS6 were appropriate for temperature and inorganic compound treatments. In contrast, the least stable reference genes were Actin1 during development stages and all other treatments, GAPDH for temperature treatments, and α-Tub for inorganic compound treatments. One target gene, Mh-Hsp90, was used to verify the selection of reference genes, results showed Mdh and RpS6 could be used as suitable reference genes for M. hapla, and Mdh plus RpS6 were better. Our finding contributes to further work on gene transcription analysis in M. hapla.


Introduction
Quantitative real-time polymerase chain reaction (qRT-PCR) is an important conventional method for measuring gene expression in molecular biology applications and has several advantages, including high sensitivity, wide dynamic range, and low cost [1][2][3][4]. However, experimental error can be caused by poor quality and low concentrations of RNA and cDNA [5][6][7][8][9][10]. In order to reduce error and achieve reliable results, reference genes, called housekeeping genes, are essential for normalising gene expression [11].

Nematode culture and treatments
The northern root-knot nematode M. hapla was maintained on susceptible tomatoes (L-402) in a greenhouse as described by Forge and MacGuidwin [34]. The eggs were extracted from tomato roots [26]. Second stage juveniles (J2) were collected 48 hr after egg hatching. Females were picked from diseased roots for the experiment. Development stage treatments. M. hapla eggs, J2 and females were transferred to 1.5-mL Eppendorf tubes and centrifuged. Pellets weighing approximately 20 mg were collected, immediately frozen in liquid nitrogen and stored at −80˚C for analysis.
Temperature treatments. Approximately 20 mg of centrifuged J2 was collected as described above, transferred to a 30 mm diameter Petri dish containing 4 mL sterile water. The samples were respectively exposed to a low temperature (4˚C) for 12 hr, preferred temperature (25˚C) for 12 hr, and high temperatures (38˚C and 40˚C) using a programmable cooling device (TEMI990, Shanghai, China) in a temperature-control chamber. The temperature was initially set at 34˚C and then increased in increments of 0.5˚C/min to 38˚C and 40˚C; samples were held for 30 min at the high temperatures. The samples were then cooled to 34˚C by decreasing the temperature in increments of 0.5˚C/min. After removal of the samples from the temperature chamber, the liquid supernatant was discarded, and the pellets were immediately frozen in liquid nitrogen and stored at −80˚C.
Inorganic compound treatments. Approximately 20 mg J2 was collected in a 1.5 ml tube and 500 uL inorganic compound (6 mM NH 4 HCO 3 , 0.77 mM FeCl 3 �6H 2 O, 0.16 mM CuCl 2 �2H 2 O, and 0.16 mM CuSO 4 �5H 2 O) was added to each tube [33]. The samples were incubated at 25˚C under dark conditions for 24 hr, rinsed for five times with RNase-free water, frozen, and stored at −80˚C. Three independent biological replicates were evaluated for each treatment.

Total RNA isolation and cDNA synthesis
Total RNA was extracted using a MiniBEST Universal RNA Extraction Kit (TaKaRa, Dalian, China) according to the manufacturer's instructions. The concentration and purity were determined twice for each RNA sample by NanoVue, and samples with an A 260 /A 280 ratio between 1.9 and 2.2 were used for cDNA synthesis. Five hundred nanogram of RNA was reverse-transcribed into cDNA in a final volume of 10 μL using a PrimeScript RT Master Mix (TaKaRa). The cDNA was serially diluted 10-fold (10×, 10 2 ×, 10 3 ×, 10 4 ×, and 10 5 × dilutions) to assess the amplification efficiency (E%) of primers and correlation coefficients (R 2 ) or 5-fold for qPCR analysis.

The primer design and qRT-PCR method
The eleven candidate reference genes were AK, Actin1, EF1-α, GAPDH, Mdh, Pur, RpS6, TAF, α-Tub, Ubp, and SDHA. The EST sequences of candidate genes and the mRNA sequences of the target gene (Mh-Hsp90) were obtained from the GenBank database. Twelve pairs of specific primers were designed by Primer Premier 5 according to the design parameters of qPCR primers with 55-65˚C melting temperature, 18-23 bp primer length, 30-55% GC content and 90-260 bp product length.
qPCR was performed using SYBR Premix Ex Taq II (TaKaRa) following the manufacturer's protocol on a Bio-Rad CFX-96 real-time PCR system (Bio-Rad, Hercules, CA, USA). Each 10-μL qPCR mixture included 5 μL SYBR Premix Ex Taq II, 1 μL diluted cDNA, 0.4 μL of each primer (10 μM), and 3.2 μL ddH 2 O. The reaction conditions were as follows: 95˚C for 30 s, followed by 40 cycles of 95˚C for 5 s and 60˚C for 30 s, and a final melt curve from 65˚C to 95˚C with a 0.5˚C increment. Each treatment included three technical and biological replicates.
Four data analysis methods, the ΔCq method, geNorm [21], Normfinder [22], and Bestkeeper [23], and the website tool RefFinder were used to evaluate the best reference genes. The data were directly analysed with BestKeeper, but converted into relative quantities for geNorm and Normfinder via the formula 2 -ΔCq , ΔCq = the corresponding Cq value − minimum Cq [36]. The geNorm program calculated the expression stability value (M) for each gene, and the pairwise variation (V n / V n + 1 ) to find the optimal number of reference genes. Using stability values, each candidate gene was estimated with NormFinder software. The lowest M-value indicated the highest stability. The gene with the highest stability was ranked with BestKeeper based on the standard deviation value and coefficient of variation value. Finally, a web-based analysis tool RefFinder was used to comprehensively estimate the best candidate genes from each program [37][38][39].

Analysis of reference gene validation
The expression level of the heat-shock protein 90 gene (Hsp90) involved in the regulation of environmental stress, is increased in nematodes under abiotic stresses, such as heat shock stress and inorganic compound stress [40,41]. Therefore, the expression levels of Hsp90 in M. hapla (Mh-Hsp90) were evaluated to validate the identified reference genes. The qRT-PCR primer pairs for Mh-Hsp90 were 5'-TTGCTAAATCTGGCACGAAGG -3' (forward) and 5'-ATGAAGGAACCACCAGCAGA -3' (reverse).

Primer amplification efficiency and specificity of candidate reference genes
The eleven candidate reference genes with accession numbers, primer pair sequences, annealing temperatures, PCR product lengths, the Tm of products, amplification efficiencies (E%), and correlation coefficients (R 2 ) are listed in Table 1. The E% values for the eleven candidate reference genes ranged from 93.6% to 106.7%, and the R 2 values reached 0.99 (Table 1). A single peak in the melting curve showed specific amplification of all primers (Fig 1).

Cq value analysis of candidate reference genes
Cq values indicate the expression levels of reference genes. The distributions of all Cq values for all samples are shown in Fig 2. The Cq values ranged from 14.91 to 29.34 for the eleven candidate reference genes, and the mean values for Actin1 and TAF were 18.59 and 25.7, respectively. Low Cq values indicate high expression levels. Among the eleven reference genes, Actin1 showed high expression, whereas TAF showed low expression.

Expression stability of candidate reference genes under different treatments
Five methods (ΔCq method, geNorm, NormFinder, BestKeeper, and RefFinder) were used to evaluated the stability of eleven candidate reference genes. Each reference gene was subjected to ten treatments, and the stability of them analysed individually ( Table 2). In addition, these ten treatments were divide into four groups for more comprehensive analysis: "Development stage (egg, J2 and female), Temperature treatments (4, 25, 38 and 40˚C), Inorganic compound treatments (CuSO 4 �5H 2 O, FeCl 3 �6H 2 O, CuCl 2 �2H 2 O, and NH 4 HCO 3 ), and All treatments (composed of all the treatments sets). The ranks of the eleven genes for groups were calculated and shown in Table 2.
Development stage treatments. Mdh was the stable gene used by ΔCq method, geNorm and RefFingder in egg treatments; The stable gene was Actin 1 in female treatments and α-Tub in J2 treatments through ΔCq method, NormFinder and RefFingder. For the groups, SDHA, RpS6 and Mdh was identified as the stable gene by the ΔCq method and Normfinder; Mdh and SDHA were identified as the most stable genes by geNorm; and Ubp, EF1-α, and α-Tub were identified as the most stable genes by BestKeeper. Therefore, combining all four rankings by RefFinder, SDHA, RpS6, and Mdh were considered the most stable genes, and α-Tub, AK, and Actin1 were considered the least stable genes ( Table 2). The optimal number of reference genes defined by geNorm shown in Fig 3. The V 2/3 values was less than 0.15 among egg, female, J2and development stage treatments. Therefore, the best reference gene combination was Mdh and TAF, GAPDH and SDHA, GAPDH and TAF, Mdh and SDHA, respectively.
Temperature treatments. Pur and Mdh were the most stable gene in 4˚C treatment by ΔCq method, NormFinder, BestKeeper and RefFingder. Mdh was the most stable gene in 38 and 40˚C treatments by ΔCq method, NormFinder, and RefFingder; For the temperature treatments, EF1-α and α-Tub were identified as the least stable genes by the other three analysis programs except BestKeeper. Mdh and RpS6 were the most stable genes identified by the ΔCq method and geNorm analysis, and Pur, AK, and Mdh showed stable expression in Normfinder analysis. In contrast, α-Tub, Ubp, and EF1-α were identified as stable genes using Best-Keeper analysis. Taken together, these results suggested that Mdh and RpS6, were the most stable genes, whereas EF1-α, and GAPDH were the least stable genes by RefFinder (  (Fig 3).       Table 2). V 2/3 values of less than 0.15 indicated that two reference genes, SDHA and Mdh, were sufficient for normalisation in all treatments (Fig 3).

Validation of stable reference genes
The stable genes Mdh and RpS6, the unstable gene GAPDH and Ubp, and the combined group of Mdh + RpS6 were selected to normalise the expression of Mh-Hsp90 under temperature (40˚C and 4˚C) and inorganic compound (CuSO 4 �5H 2 O) treatments. Generally, the use of multiple reference genes presents more accurate normalization of the gene expression. Results showed in Fig 4, the relative expression levels of Mh-Hsp90 were similar when normalised using Mdh, RpS6, and Mdh + RpS6, but different when normalised using GADPH and Ubp.

Discussion
qRT-PCR is a powerful technique with high sensitivity and specificity and enables gene expression analysis within a large dynamic range [1,42,43]. The amplification efficiency (E%) of reference genes should be similar to that of the target gene, which are essential for improving the accuracy of gene expression. Additionally, an optimal reference gene should show moderate and stable expression levels in all test samples. In our study, the amplification efficiencies (E%) of stable candidate reference genes were similar to that of Mh-Hsp90 on the same treatments. Moreover, except for Actin1, α-Tub, and TAF, all other genes showed moderate Cq values, and all candidate reference gene were specifically amplified. The values of candidate reference genes were slight difference between programs (geNorm, NormFinder, Bestkeeper) under different treatments, that caused by the difference of algorithms employed [22]. Therefore, it is better to evaluate reference gene using multiple methods, and then determined the suitable  Evaluation of suitable reference genes for gene expression analysis in Meloidogyne hapla reference genes with the geometric mean of comprehensive ranking for all programs which generated by RefFinder. The stable reference genes have been identified in Caenorhabditis elegans that were tba-1, Y45F10D.4 and pmp-3 for studying nanoparticle-induced genetic response [44], cdc-42, pmp-3 and Y45F10D.4 for normalizing 5 sod expression levels [45]. In this paper, the Mdh and RpS6 were reliable gene in M. hapla, consistent with the RPS15 in Helicoverpa armigera [37], RPS4 in turbot [46], RpS6 in Macrobrachium olfersii [47] for different developmental stages, and similar with the RPS20 for Sesamia inferens [7], RPS15 and RPS27 for H. armigera [37] under different temperature treatments. EF1-α, and GAPDH were not good reference genes in M. hapla for temperature stress, while EF1-α was useful reference gene for parsley [39]. According to these findings, there were no absolute reference genes for different species and treatments. While, the qRT-PCR relied on accurate normalization of stable reference genes. Therefore, the stability of reference genes should be validated for different experimental condition before use.
Inorganic compounds affect the survival of J2 and the hatch rate of egg-masses for M. hapla [33]. The responses of M. hapla to inorganic compound stress are still unknown. Analysis of gene expression will enable researchers to study the effects of inorganic compound stresses on M. hapla. In this study, we identified suitable reference genes (RpS6 and Mdh) for normalising gene expression. These findings are expected to facilitate further analyses of the mechanisms of inorganic compound stress in M. hapla.
In conclusions, this work validated that RpS6, Mdh, SDHA and Pur could be used as suitable reference genes for normalising qRT-PCR data in M. hapla under different treatments, and the combination of RpS6 + Mdh were better. This study provides a basis for future studies of gene function in M. hapla.