There is no magic bullet: the importance of testing reference gene stability in RT-qPCR experiments across multiple closely related species

Study insects and rearing regime

We used four closely related species in the genus Schistocerca, maintained as laboratory colonies in the Department of Entomology at Texas A&M University. The four species were the Central American locust, S. piceifrons, and three sedentary species, S. americana (Drury), S. serialis cubense (Saussure), and S. nitens (Thunberg). For conciseness, we use the species epithet (piceifrons, americana, cubense, and nitens) throughout the paper. The piceifrons colony was established from an outbreak population in Yucatan, Mexico collected in October 2015, and imported under a USDA permit (USDA APHIS PPQ P526P-15-03851). The americana colony was established from a population in Brooksville, Florida, collected in September 2010. The cubense colony was established from a population in Islamorada in the Florida Keys collected in January 2011. Finally, the nitens colony was established from a population in Terlingua, Texas, collected in May, 2015. For this study, the colonies of these four species were reared under crowded and isolated conditions. For the isolated condition, nymphs were isolated as hatchlings and reared in separate plastic cages (10.16 × 10.16 × 25.4 cm) with one transparent side and connected to a filtered positive airflow. For the crowded condition, both cage size and the number of specimens depended on the species. For piceifrons, about 200 nymphs were kept in a large cage (40.64 × 34.29 × 52.07 cm) in a USDA approved quarantine facility. For americana and cubense, 150–200 nymphs were kept in larger cages (30.48 × 35.56 × 50.8 cm). For nitens, over 50 individuals were kept together in a small cage (30.48 × 35.56 × 50.8 cm) because we observed that the insects would die in the density used for other species. In both density conditions, the insects were reared at 12 h of light and 12 h of darkness at 30 °C, and were fed daily Romaine lettuce and wheat bran. We reared the insects until they molted to the last nymphal instars to conduct qPCR experiments.

RNA extraction and cDNA synthesis

Sample collection, RNA extraction and RNA quality assessment were performed as previously described (Wang et al., 2020). Briefly, crowded-reared and isolated-reared female nymphs were marked with a ceramic marker on the abdomen after molting to their last nymphal instar, and were dissected around 72 h later. Only specimens that molted before 10 AM were used, and all dissections occurred between 8 and 9 AM. After removing gut tissues, head and thorax were dissected using sterilized dissection scissors. Both tissues were preserved in RNAlater (ThermoFisher Scientific) at −20 °C, following the manufacturer’s guideline. A total of 10 nymphs/density/species were dissected. Half of these was used for RNA sequencing, and the other half was used for the qPCR experiment. Tissues were placed in MagNA Lyser Green Beads (Roche) sample tubes and were homogenized for 30 s in one mL of Trizol (ThermoFisher Scientific) using a MagNa Lyser instrument (Roche) at 6,500 rounds per minute. Whole-tissue RNA was extracted using a Trizol-chloroform extraction, followed by clean-up with a RNeasy mini kit (Qiagen) using an on-column DNAse treatment with a RNase-free DNAse set (Qiagen). RNA concentrations were measured with a Denovix DS-11 spectrophotometer; 260/280 and 260/230 values were above 2.0 for all samples. Additionally, microcapillary electrophoresis with a Fragment Analyzer (Agilent Technologies RNA) was used to analyze RNA integrity of samples destined for RNA sequencing. Only those samples with a RNA Quality Number (RQN) over 3.9 were used. Due to differences in 28S ribosomal RNA structure compared to other eukaryotic species, RQN values for insects are often lower than what is generally considered a valid threshold in mammalian samples (Escobar & Hunt, 2017; Macharia, Ombura & Aroko, 2015; Winnebeck, Millar & Warman, 2010). Samples used for qPCR were diluted to a concentration of 100 ng*ml⁻¹ and subsequently used to synthesize cDNA using the iScript cDNA synthesis kit (Bio-Rad) following the manufacturer’s guideline.

RNA sequencing and transcriptome assembly

Both RNA sequencing and transcriptome assembly used for the present study were previously described in detail in Wang et al. (2020). Briefly, we generated RNA-seq data by performing paired-end sequencing (150 bp) on 8.5 lanes on an Illumina HiSeq4000 (San Diego, CA). Library preparation, sequencing, and read formatting was performed at Texas A&M’s AgriLife Research Genomics and Bioinformatics Service. After initial processing of raw data, raw reads were imported into a personalized Galaxy environment (Afgan et al., 2018) on a supercomputing cluster of the High-Performance Research Computing group of Texas A&M University (Ada, https://hprc.tamu.edu) for further processing, trimming, and quality control. FastQ Screen (Wingett & Andrews, 2018) was used to filter out sequences from the following potential contamination sources: (UniVec core (June 6, 2015), PhiX (NC_001422.1), Illumina adapters, Gregarina niphandrodes genome (GNI3), Encephalitozoon romaleae genome (ASM28003v2), Escherichia coli genome (K12), Methylobacterium sp., Bosea sp., Bradyrhizobium sp., Klebsiella pneumoniae, Sphingomonas sp., Rhodopseudomonas sp. and Propionibacterium acnes). We used Trinity (Grabherr et al., 2011) for de novo assembly, which was performed separately for each species and tissue. We further filtered the resulting assemblies using CD-hit-EST (Fu et al., 2012; Li & Godzik, 2006) and Transrate (Smith-Unna et al., 2016). We assessed transcriptome quality using Trinitystats (Grabherr et al., 2011), BUSCO (Benchmarking Universal Single-Copy Orthologs; (Simão et al., 2015), and by calculating the fraction of reads mapping back to the transcriptome with a combination of bowtie2 (Langmead & Salzberg, 2012; Langmead et al., 2009) and flagstat (Li, 2011a; Li, 2011b; Li et al., 2009).

Sequence assembly and expression analysis

All eight transcriptomes were imported into Geneious (R10.2.6; BioMatters, Ltd.). We selected nine potential reference genes: actin5C (Act5C), α-tubulin (Tub), succinate dehydrogenase (SDH), elongation factor 2 (EF2), ribosomal protein L5 (RIBL5), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), annexin IX (Ann), armadillo (Arm) and heat shock protein 70 (Hsp70). We also selected four target genes: pacifastin-related peptide 4, pacifastin-related peptide 5, allatostatin (ast), and allatotropin (at). For consistency, we followed suggestions by Simonet et al. (2004) for naming pacifastin-related peptides and named the Schistocerca piceifrons sequences SPPP-4 and SPPP-5. We obtained nucleotide sequences for all 13 genes from other orthopteran insects from Genbank (https://www.ncbi.nlm.nih.gov/genbank/) and used Megablast or tblastx with default settings in Geneious to find related sequences in our four species. For each gene, sequences of all four species were aligned using MUSCLE version 3.8.24 (Edgar, 2004) in Geneious, with a maximum of nine iterations and default settings. Subsequently, sequences were manually curated, resulting in full-length coding sequences for each species. The identity of each gene was confirmed using the standard nucleotide Basic Local Alignment Search Tool (blastn) at the NCBI website, using the nr/nt database as reference (Altschul et al., 1990; Johnson et al., 2008). Sequence information can be found in Table 1. Additionally, mapping reads were counted as previously described (Wang et al., 2020), using bowtie2 (very sensitive end-to-end, disable no-mixed and no-discordant behavior) and SAMtools’ idxstats (Li et al., 2009) in Galaxy. Differential expression analysis was performed with DEseq2 (Love, Huber & Anders, 2014) within SARTools 1.7.1 (Varet et al., 2016), using R 3.5.3 (R Core Team, 2017). All settings were kept at their defaults; the expression in crowded-reared individuals was compared to that of conspecific isolated-reared individuals.

Primer design

Primers for reference genes (Table 2) were designed using Primer3 (Koressaar & Remm, 2007; Untergasser et al., 2012), based on conserved regions, identified using the nucleotide alignment in Geneious. In this way, we aimed to generate primers that would work in all four species. Standard settings were altered to an optimal melting temperature set at 60 °C, and the amplicon length at 150–200 bp. For the four target genes, primers were designed for just piceifrons, with identical settings in Primer3 as described above (Table 2). We designed three primer pairs for each gene. All sequences were ordered from Integrated DNA Technologies.

Real time quantitative PCR and statistics

For each qPCR reaction, 5 µl of cDNA was added to 10 µl of SsoAdvanced™ Universal SYBR^® Green Supermix (Bio-Rad #1725275) and 5 µl of primers at a final concentration of 500 nM. Every reaction was performed in duplicate or triplicate. To determine primer efficiency for reference genes, a dilution series of one sample was generated for each species by performing serial 10-fold dilutions ranging from dilutions of 1/1 to 1/10,000. For the target genes, a 5-fold dilution series from 1/1 to 1/3,125 was used instead for determining primer efficiency. Based on these data, the most efficient primer pair was used in all further analyses. All other reactions were run on 96-well plates, with 10 samples (five isolated-reared individuals and five crowded-reared individuals) run on the same plate for a particular gene following the sample maximization method (Hellemans et al., 2007). Reactions were run on a CFX connect real time system (Bio-Rad) using the following thermal cycling profile: 3 min at 95 °C, 40 cycles of (1) 15 s at 95 °C and (2) 45 s at 60 °C, and a melting curve from 65 °C to 95 °C. C_q values were exported from the Bio-Rad CFX manager using the default threshold.

For the analysis of reference gene stability, a total of five crowded-reared and five isolated-reared individuals per species were used for both head and thorax tissues. Gene stability was analyzed using three different programs: geNorm in qbase+ (Mestdagh et al., 2009; Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Orntoft, 2004), and BestKeeper (Pfaffl et al., 2004). Raw C_q values were used as input for BestKeeper and geNorm, while they were transformed to a linear scale for NormFinder. Gene efficiencies were assumed to be 100% for all tested genes (see Table 2 for actual gene efficiencies). To obtain a non-arbitrary ranking of reference genes, they were first ranked for each program separately, and subsequently the geometric mean of these scores was calculated as an accurate estimate of which reference genes were the most stable (Chen et al., 2011; Chen et al., 2013; Gong et al., 2016; Pereira-Fantini et al., 2016). The amount of reference genes needed to standardize gene expression was assessed using geNorm. Normalization factors were calculated in a stepwise manner, starting by taking the two most stable reference genes into account and adding another reference gene one by one. Subsequently, the variation between these normalization factors were compared. If the difference in variation is lower than 0.15, the addition of an additional reference gene was deemed unnecessary (Mestdagh et al., 2009; Vandesompele et al., 2002).

For the quantification of target gene expression, we used thorax tissues of five isolated-reared and five crowded-reared individuals of piceifrons. Relative expression of target genes (ast, at, SPPP-4, and SPPP-5) for thorax tissue in piceifrons was calculated with the ΔΔC_q-method (Livak & Schmittgen, 2001), using different sets of reference genes to normalize gene expression. P-values were calculated with a two-tailed student t-test in R (version 3.6.2) based on non-transformed ΔC_q-values. All raw qPCR data can be found in Data S1–S3, and the code for the student t-test is present in Data S1.

Selection of potential reference genes

A total of 9 potential reference genes were initially selected based on previous studies using qPCR for locust gene expression research (Chapuis et al., 2011; Van Hiel et al., 2009; Yang et al., 2014). All of these genes were members of different functional categories (Table 1), with the exception of the cytoskeleton genes, Tub and Act5C, which are both commonly used as reference genes in insect gene expression studies (Lü et al., 2018). We chose EF2 and RIBL5 over elongation factor 1α and other ribosomal genes, respectively, based on the low amount of variation for these genes in our transcriptome data. For 7 out of 9 genes, primer efficiencies for the selected primer pairs varied between 96 and 106% (Table 2). Two exceptions were Arm with an efficiency of 109.8% for piceifrons and 92.1% for nitens, and SDH for which we were unable to find primers with sufficient efficiency values in all four species. We attribute this to a surprisingly high amount of sequence differences among the four species for SDH, effectively restricting the sequence regions that were conserved enough for primer design. As a result, SDH was removed from all further analyses, and thus, only eight reference genes were used in further analyses. All included primer pairs showed melt curves with a single peak, suggesting the absence of aspecific amplification and contamination. C_q values varied between 17.5 and 28.5, with GAPDH and Hsp70 having the lowest values and Arm having the highest.

Comparing reference gene stability

We subsequently compared the stability of these eight potential reference genes by performing qPCR reactions in five isolated-reared and five crowded-reared individuals for head and thorax tissues in all four species (Fig. 1). In americana, cubense, and nitens, all eight genes showed similar expression levels under the tested rearing conditions. However, in piceifrons, several genes showed a trend towards differential expression between the two density conditions. Specifically, Ann, Act5C, Tub, and Hsp70 showed a trend towards lower expression in the isolated condition for both tissues, while Arm and EF2 showed this trend only in the thorax tissue. RIBL5 was the only gene showing no difference between the two density conditions. We were able to confirm most of these observed trends with our transcriptome data (Table 3). The only discrepancies between our transcriptome data and the qPCR data were found for Hsp70 in head tissue, GAPDH for thorax tissue, and Tub for both tissues. For thorax tissue, Act5C and Ann were even found to be significantly upregulated in crowded-reared individuals compared to isolated-reared individuals.

NormFinder (Andersen, Jensen & Orntoft, 2004) calculates both intra- and inter-group variations of gene expression, and combines these into a stability value, and thus, genes with a lower rank are considered to be more stably expressed. geNorm (Mestdagh et al., 2009; Vandesompele et al., 2002), the most popular algorithm, assumes that expression values of real reference genes are perfectly correlated with each other in all tested samples. Based on this correlation, genes get a stability value (M value), with genes with M < 1.5 considered to be stable and the most stable gene obtaining the lowest score. BestKeeper (Pfaffl et al., 2004), on the other hand, outputs the standard deviation for each reference gene, expected to be lower than 1, and a correlation coefficient for each gene with the so-called BestKeeper-index, which is essentially the geometric mean of all stable reference genes. Before comparing the stability values, we estimated the amount of necessary reference genes for each species and tissue using geNorm’s V-values. For americana, cubense, and nitens, geNorm suggested that using two reference genes was sufficient as its values for V_2∕3 were below 0.15 (Fig. 2). However, for piceifrons, the difference in variation between using two or three reference genes was higher than 0.15 for the head tissue and very close to 0.15 for the thorax tissue (Fig. 2). Thus, the use of three reference genes is suggested for this species. Overall, the three programs listed genes in a similar fashion, with only a few exceptions (Tables 4 and 5). Tub and GAPDH were the two least stable genes, and both were listed as the last in 3 out of the 8 tested species/tissue combinations. Act5C be the most stable reference gene overall, as it was generally listed in the top half of stability values by all three programs. Arm was also often listed in the top three of most stable genes, but was also listed as the seventh most stable gene in some cases (Normfinder, Bestkeeper in piceifrons head; Normfinder in nitens head).

For the americana head tissue, Act5C and EF2 were the two most stable genes overall, with Act5C rated the most stable gene by each program, and EF2 obtaining a second place for both geNorm and NormFinder. RIBL5 and Arm, which were also consistently ranked as stable reference genes, obtained a third and a fourth place. For the americana thorax tissue, RIBL5 and EF2 were our two reference genes of choice, as they obtained a first and a second place overall. Additionally, Arm and Act5C were ranked as third and fourth overall, but obtained very similar stability scores to EF2 and were ranked in the first half by all three programs. For the cubense head tissue, both Act5C and Hsp70 were ranked within the top 3 most stable reference genes by all three programs and as a result, obtained top overall rankings. For the cubense thorax tissue, Hsp70 obtained the best ranking, followed by Ann ranked second. Nonetheless, it should be noted that for this species and tissue combination, several reference genes obtained very similar scores for both NormFinder and geNorm, resulting in low rankings for some genes that were actually quite stable. For instance, Act5C was ranked as only fifth overall even though it performed very similar to Ann, which was ranked second. For the nitens head tissue, the rankings suggested by the three different programs disagreed strongly with each other, as Ann and Act5C were ranked first and second by geNorm and NormFinder, but fifth and fourth by BestKeeper, while Arm was ranked first by BestKeeper and sixth and seventh by the other two programs. Overall, Ann and Act5C were ranked first and second, and Hsp70, which was the only gene ranked within the first half by all three programs, was ranked as third. For the nitens thorax tissue, Ann and Arm were clearly the preferred reference genes as they were consistently ranked as the two most stable genes by all three programs. Finally, for both tissue of piceifrons, Act5C and Ann were ranked in the top 4 of most stable genes by all programs. RIBL5 and Arm were also ranked in the top 4 for head and thorax tissue respectively.

Validation of reference genes

To explore how the choice of reference genes would affect the outcome of qPCR experiments, we quantified and compared the expression levels of the four target genes (ast, at, SPPP-4, and SPPP-5) using different reference genes in the thorax tissue of piceifrons reared in two density conditions. As Tub was the least stable reference gene, it was not included for this analysis. Our data showed the choice of reference genes strongly affected the qPCR results (Fig. 3). Regardless of the reference genes used, ast showed a significantly higher expression level in isolated-reared individuals compared to crowded-reared individuals, and SPPP-5 never came close to reaching significant differences. When all of the reference genes (except Tub) were chosen, at and SPPP-4 showed trends toward decreased and increased expression, respectively, in crowded-reared individuals (p = 0.068 and p = 0.061, respectively). We subsequently removed RIBL5 from the analysis and thus only used reference genes of which the expression varied between both rearing conditions. As a result, the p-value for at dropped below 0.05, while the trend observed for SPPP-4 partially disappeared. Lastly, we included just Hsp70, RIBL5 and GAPDH as reference genes, which was our preferred set of reference genes based on the earlier data. As a result, SPPP-4 showed a significantly higher expression level in the crowded-reared condition, while no significant difference was observed for at anymore (p = 0.19). Three out of the four of these qPCR results agreed with our transcriptome data (Table 6).