Evaluation of potential reference genes for real-time qPCR analysis in a biparental beetle, Lethrus apterus (Coleoptera: Geotrupidae)

Sample collection

Samples were collected near Dorogháza, northern Hungary (47°59′29″N, 19°53′36″E) on 16th April, 4th May and 28th May in 2015, which dates corresponded to the beginning, middle and end of the breeding season of Lethrus apterus, respectively. Sample collection was approved by the Nothern Hungarian Inspectorate for Environment Protection and Nature Conservation (No. 9007-8/2014). The first sampling date represents the period of mate choice, while the second and third samplings were done during the period when parents were collecting leaves for the offspring. On each sampling dates, head and thorax samples were collected from eight males and eight females. All tissues were removed from the head capsule and muscle samples were taken from the thorax. Samples were collected in the field in less than five minutes after euthanizing the individuals. Each head and thorax sample was put immediately into separate eppendorf tubes which already contained 600 µl RNAlater^® Stabilization Solution (Thermo Fisher Scientific, Waltham, MA, USA), then stored at −20 °C in the laboratory in order to inhibit RNase enzyme activity until RNA extraction.

RNA extraction and cDNA synthesis

Total RNA was isolated from each samples using TRIzol^® Reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. The extracted RNA was eluted in 15–30 µl RNase-free water, depending on the pellet size. Yield of RNA was quantified by NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). To eliminate genomic DNA, samples were treated with RQ1 RNase-Free DNase (Promega, Madison, WI, USA) just before the reverse transcription. First strand cDNA was synthesized from 1 µg DNA-free RNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA).

Reference gene selection and primer design

Using a draft genome of Lethrus apterus (unpublished data) eleven reference genes, which were already described as stable reference genes in other arthropods, were selected (Table 1). We manually designed primers (Table 2) using the web-based Sequence Manipulation Suite (Stothard, 2000) and Multiple Primer Analyzer (Thermo Fisher Scientific, Waltham, MA, USA) in order to avoid the forming of possible secondary structures of the primers. To check the specificity of primer pairs and to determine optimal annealing temperature, PCR reactions were performed in 10 µl volumes containing the following components: 10x buffer, 2 mM MgCl₂, 0.2 mM dNTP, 0.02 U/µL Taq DNA polymerase enzymes (DreamTaq Green, Thermo Fisher Scientific, Waltham, MA, USA), 0.2 µM forward and 0.2 µM reverse primer and 0.1 µg cDNA. PCR conditions were optimized by determining the optimal annealing temperature using temperature gradient ranging from 54 °C to 62 °C for primer binding. In this study, we used ABI Veriti^® 96-Well Thermal Cycler (Applied Biosystems, Foster City, CA, USA). Cycling conditions consisted of a denaturing step at 95 °C for 2 min followed by 40 cycles at 95 °C for 30 sec, at a temperature gradient (54 °C, 56 °C, 58 °C, 60 °C or 62 °C) for 30 s and at 72 °C for 90 s, and finally at 72 °C for 10 min. PCR amplicons were run on 1% agarose gel stained with GelRed™ (Biotium, Fremont, CA, USA).

Real-time quantitative PCR

RT-qPCR was performed on a QuantStudio 12K Flex Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using SYBR^® Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) and ROX Passive Reference Dye (Affymetrix, Santa Clara, CA, USA). Amplifications were carried out under the following conditions: initial denaturation at 95 °C for 10 min followed by 40 cycles of 10 sec at 95 °C and for 1 min at the optimal annealing temperature. This was followed by a melting curve analysis in which the temperature raised from 65 °C to 95 °C in sequential steps of 0.05 °C for 1 s. Three technical replicates were performed for each biological sample, and the average cycle threshold (Ct) values of triplicates were calculated. Furthermore, no-template control was done in order to check whether primer-dimers or contamination with amplified PCR product were detectable. Five 5-fold serial dilution was made from cDNA samples to create a standard curve, and the amplification efficiency was determined for each candidate gene. The efficiency (E) values were calculated according to the equation: E = (10^{(−1∕slope)} − 1) × 100, where slope is the slope of the standard curve (Radonić et al., 2004).

Statistical analysis of raw Ct values

In order to examine the differences between sample groups, random intercept mixed-effects models were used with sample id as a random factor for each gene. Significance of fixed terms was investigated by likelihood ratio tests. For likelihood ratio tests models were fitted using Maximum Likelihood estimation. The analyses were carried out using “lme4” (Bates et al., 2014) and “car” (Fox & Weisberg, 2011) packages in the R statistical environment version 3.3.2 (R Core Team, 2016).

Determination of reference gene expression stability

In order to determine the expression stability of the selected reference genes, we used the following methods: geNorm (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Ørntoft, 2004), BestKeeper (Pfaffl et al., 2004), delta Ct method (Silver et al., 2006) and RefFinder (Xie et al., 2012). For the analyses with the geNorm and NormFinder procedures, the average Ct values were transformed to relative quantities by dividing sample values by the lowest average Ct value. For calculations by BestKeeper, delta Ct method and RefFinder, the untransformed average Ct values were used. All calculations, except the ones done by the web-based RefFinder, were carried out in R with “NormqPCR” package (Perkins et al., 2012).

geNorm calculates the expression stability value M by assessing the mean pairwise expression ratio for each candidate gene against all the other candidates (Vandesompele et al., 2002). The basic assumption of this method is that the expression ratio between two reference genes is identical across the samples. The lower the M value the more stable the expression of the candidate reference gene. Stepwise exclusion of the genes with the highest M value results in the selection of the two most stably expressed reference genes in the tested samples both sharing the same M value. Vandesompele et al. (2002) also suggest not to accept candidate genes as stably expressed reference genes with M value higher than 1.5. Moreover, the procedure determines the normalisation factor by taking the geometric mean of the expression levels from the most stable genes and then additively recalculating with each of the next most stable gene. The pairwise variation, V_n∕n+1 between two sequential normalisation factors is then calculated in order to determine the effect of each newly added gene to the normalisation factor. The optimum number of genes is the lowest number of genes with V_n∕n+1 less than 0.15 (Vandesompele et al., 2002).

NormFinder determines the stability of the candidate reference genes by measuring the intra- and intergroup variation between user specified groups (e.g., male and female groups or treated and control groups) first. Stability values for each candidate gene are then calculated by adding the two sources of variation. The lowest stability value means the most stable expression (Andersen, Jensen & Ørntoft, 2004).

BestKeeper calculates, for each candidate reference gene across the samples, the geometric mean, the arithmetic mean, the minimal and the maximal Ct values, in addition to the average absolute deviation from the arithmetic mean. Genes with the lowest average absolute deviation can be considered as stably expressed reference genes. BestKeeper Index is calculated as the geometric mean of the Ct values of the candidate reference genes. Inter-gene relations are estimated by performing pairwise correlation analyses of all possible reference gene pairs. Furthermore, correlation between the expression level of each candidate gene and the BestKeeper Index is calculated, describing the relation between the index and the contributing genes by the Pearson correlation coefficient, coefficient of determination and the corresponding p-value (Pfaffl et al., 2004).

The delta Ct method compares relative expression of pairs of candidate genes within each sample in order to identify the stably expressed housekeeping genes. If the ΔCt value of the two genes fluctuates when analysed in different samples, it means that one or both genes are variably expressed. If the ΔCt value remains constant, both genes are stably expressed among the samples (Silver et al., 2006).

Each procedure mentioned above uses different algorithms to calculate an expression stability value which represents the suitability of the candidate genes as reference genes, therefore the ranking of the examined genes according to the methods may vary. The web-based tool RefFinder (Xie, 2012) was used in order to combine our results and rank the candidate genes. This user-friendly program integrates the four methods mentioned above. Using the ranking from each program, it assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall ranking. The lowest rank indicates the most stably expressed gene (Xie et al., 2012).

For each analysis, except for NormFinder, seven sample groups were used: all samples irrespective of body part or sex; head samples irrespective of sex; male head samples; female head samples; thorax samples irrespective of sex; male thorax samples; female thorax samples. The calculation by NormFinder requires subgroup specification, therefore, body regions were set as subgroups for the analysis of all samples. In order to investigate the effect of sexes, male and female subgroups were specified for the analysis of head and thorax sample groups separately. In this way, three sample groups, each divided into two subgroups, were analysed by NormFinder.

Transcriptional profiling of candidate reference genes

Before the evaluation of expression stability of the eleven candidate genes, specificity of each primer pair was checked on 1% agarose gel which showed single products with the expected sizes. Moreover, gene-specific amplification was confirmed by single melting curve peaks. These results indicate that no primer-dimers or nonspecific amplification products were formed. Additionally, no fluorescent signals were detected in the negative control during the RT-qPCR. Each amplicons were sequenced and annotated to the sequences from which the primer design was based in order to check that the correct genes were amplified. The sequences are available in File S1. The efficiency of the eleven candidate genes ranged from 99.33 to 112.70%. The efficiency values and other basic information of the RT-qPCR required based on the guideline of Bustin et al. (2009) are included in Table 2.

Raw Ct values ranged from 11.66 (TUB1a) to 30.12 (ARF4) (Fig. 1). The mean and standard deviation (SD) of the Ct values across all samples were calculated for each gene (Table 3). Since the mean Ct values ranged between 15 and 30 for all the candidate reference genes, all of them were analysed further (Kozera & Rapacz, 2013). ARF1 had the least variable expression level with the lowest SD value (SD = 1.85), while ARF4 had the most variable expression level (SD = 3.04). Low average Ct values indicate high expression level in TUB1a and EF2 (Ct_mean = 15.11), on the other hand, high Ct values of ARF1 (Ct_mean = 20.77) indicated low expression.

Based on the likelihood ratio tests, sex had no significant effect on the expression level of the candidate genes. However, significant effect of body region was found in case of six genes: GAPDH, EF1a, ARF1, ARF4, RP4 and L10. The interaction of sex and bodypart had no significant effect on the expression level of the candidates, except for RPS8 (Table 4).

Expression stability of candidate reference genes

Based on geNorm analysis for all samples, eight candidate genes had an M value below the threshold of 1.5 (Table S1). The results show that the lowest M value was 0.390 for RPS8 and L7A. Among the head samples irrespective of sex, all of the tested genes except L10 had an M value below 1.5, and RPS8 and RP18 were co-ranked as the most stable genes from the candidates (M = 0.304). Furthermore, the same two genes had the lowest M value considering male and female head samples separately (M = 0.264 for females and M = 0.346 for males). In case of the thorax samples irrespective of sex, eight genes had an M value below the threshold. RPS8 and L7A were the most stable candidate gene pair with an M value of 0.358. In thorax samples collected from females, RPS8 and RP18 were the most stable genes as well with an M value of 0.222. However, in thorax samples of males, RPS8 and L7A were ranked as the best reference gene pair (M = 0.288).

According to NormFinder, L7A was the most stable gene when calculating with all samples divided into groups of head and thorax samples (Table S2). The second and third genes were RP4 and RPS8, indicating that these are also worth considering as reference genes. In the case of specifying males and females as subgroups within head and thorax samples, L7A was found again to be the most stably expressed gene among the candidate ones. In both head and thorax samples, L7A was followed by similar ranking order: EF2, RP4, RP18 and RPS8 as second, third, fourth and fifth genes, respectively.

Based on BestKeeper, across all samples ARF1 had the lowest mean absolute deviation (MAD) value; however, L7A had the highest correlation r value (Table S3). In the group of head samples irrespective of sex, ARF1 had the lowest MAD value, while among the thorax samples irrespective of sex, L10 was the most stable according to the MAD value. This was surprising as the other programs ranked this gene consistently as one of the least stable genes. On the other hand, in both head and thorax samples, L7A had the highest r value. In head samples of males RPS8 (MAD = 0.827), and of females ARF1 (MAD = 1.122) were the most stable candidate genes. In both male and female head samples, L7A had the highest correlation r value. Considering male thorax samples L10 had the lowest MAD value (MAD = 1.576), while in female thorax samples EF2 was ranked as the most stable with MAD value 1.236. L7A had the highest r value in female thorax samples, however, in male thorax samples EF2 had the highest correlation r value.

According to the delta Ct method, L7A was the most stable gene among the candidates overall with the stability value always below 1.0 (Table S4).

Finally, the candidate genes were evaluated by RefFinder to combine the results of individual methods (Table 5). Using all samples irrespective of body region and sex, and separately the head and thorax samples irrespective of sex, L7A was ranked first, as the most stably expressed gene among the candidate reference genes. In head samples, RP18 was co-ranked with L7A as the most stable reference genes. In thorax samples, RP4 was ranked on the second place. In female head and thorax samples, RP18 was the most stably expressed gene of the candidates. Considering head samples of males, RPS8 was ranked on the first place, while in thorax samples of males, L7A was ranked as the best reference gene. L7A was ranked on the second place in all subgroups, with the exception of male thorax samples, where EF2 was the second best reference gene according to RefFinder.

Optimal number of reference genes

To determine the minimal number of genes necessary for normalisation, the V-value was computed by geNorm. The results demonstrated that across all samples V_2∕3 was the first V-value lower than the cut-off value of 0.15 (Fig. 2). Considering separately the head and thorax samples, V_2∕3 was again lower than 0.15. Separate analyses of female and male samples within head and thorax groups showed that V_2∕3 was also the first value below the threshold in all cases (results not shown). Therefore, two stably expressed reference genes are sufficient for normalisation in any case of sample classification.