Selection and stability validation of reference gene candidates for transcriptional analysis in Rousettus aegyptiacus

Bats are the only mammals capable of powered flight and their body temperature can reach up to 42 °C during flight. Additionally, bats display robust type I IFN interferon (IFN-I) responses and some species constitutively express IFN-α. Reference genes with stable expression under temperature oscillations and IFN-I release are therefore critical for normalization of quantitative reverse-transcription polymerase chain reaction (qRT-PCR) data in bats. The expression stability of reference genes in Rousettus aegyptiacus remains elusive, although this species is frequently used in the infection research. We selected ACTB, EEF1A1, GAPDH and PGK1 as candidate reference genes and evaluated their expression stability in various tissues and cells from this model bat species upon IFN-I treatment at 35 °C, 37 °C and 40 °C by qRT-PCR. We employed two statistical algorithms, BestKeeper and NormFinder, and found that EEF1A1 exhibited the highest expression stability under all tested conditions. ACTB and GAPDH displayed unstable expression upon temperature change and IFN-I treatment, respectively. By normalizing to EEF1A1, we uncovered that GAPDH expression was significantly induced by IFN-I in R. aegyptiacus. Our study identifies EEF1A1 as the most suitable reference gene for qRT-PCR studies upon temperature changes and IFN-I treatment and unveils the induction of GAPDH expression by IFN-I in R. aegyptiacus. These findings are pertinent to other bat species and may be relevant for non-volant mammals that show physiological fluctuations of core body temperature.


Scientific Reports
| (2021) 11:21662 | https://doi.org/10.1038/s41598-021-01260-z www.nature.com/scientificreports/ The knowledge about the exceptional immune system of bats has significantly advanced during the past decade 16 . However, experimental tools to systematically investigate bat immune responses, such as speciesspecific or cross-reactive antibodies, are largely missing [17][18][19] . Accordingly, investigations on host immunity heavily rely on gene transcription profiling by qRT-PCR. This method is sensitive, specific, highly reproducible and accurate 20,21 . A critical step in qRT-PCR setup is selection of several stable reference genes. The inclusion of such reference genes is crucial for gene expression normalization and subsequent data interpretation. The ideal reference genes should maintain stable expression levels across diverse tissues and cell types as well as under different experimental conditions 22 . Considering that bats display unique physiological features, notably oscillating metabolic rates and core body temperature depending on flying and roosting phases [23][24][25][26] , the expression stability of reference genes must be evaluated in context of these physiologically relevant conditions. Multiple reference genes, including glyceraldehyde-3-phosphate dehydrogenase (GAPDH), actin-beta (ACTB), small nuclear ribonucleoprotein Sm D3 (SNRPD3) and 18S ribosomal RNA (18S rRNA) have been employed in gene expression studies of P. alecto, E. fuscus, M. davidii and other bats 13,14,[27][28][29] . However, a comprehensive analysis of reference genes, particularly their expression stability under oscillating temperatures, has not been performed in bats, including the model bat R. aegyptiacus.
Here, we provide a first in-depth validation of four reference gene candidates, including ACTB, GAPDH, eukaryotic translation elongation factor 1 alpha 1 (EEF1A1) and phosphoglycerate kinase 1 (PGK1) for R. aegyptiacus. Our findings support EEF1A1 as an appropriate reference gene for normalization of qRT-PCR data in the Egyptian fruit bat and call for caution when using other candidate genes, i.e. GAPDH and ACTB, due to their instability under specific conditions.

Results
Performance of PCR primers targeting reference gene candidates. To evaluate the performance of the qRT-PCR assay, we first examined the specificity of the primer pairs with melting curve analysis, agarose gel electrophoresis and sequencing. Melting curve analysis revealed single peaks for all primer pairs (Fig. S1). Agarose gel electrophoresis further demonstrated single bands for all the PCR products with the predicted sizes, indicating high specificity of all primer pairs (Fig. S1). PCR products were sequenced and the specificity of the primers was confirmed (Fig. S2). To evaluate the amplification efficiency, standard curves with ten-fold dilution steps were generated (Fig. S3) and subsequently the linear dynamic range (LDR) and precision of each primer pair were assessed following the MIQE guideline 30 . Amplification efficiencies of all tested primers met the validation criteria, notably 100.51% for ACTB, 99.53% for EEF1A1, 106.52% for GAPDH and 106.44% for PGK1 ( Table 1). The correlation coefficient (R 2 ) of all candidates was above 0.99, suggesting excellent linearity of the standard curves. The LDR values of all primers were in the range of 3 to 300,000 copies and precision values varied from 0.31 (EEF1A1) to 1.27 (ACTB) ( Table 1). Thus, all the primers for reference gene candidates demonstrated satisfactory specificity and efficiency in qRT-PCR.
Expression profile of reference gene candidates in various tissues from R. aegyptiacus. To investigate the expression of reference gene candidates in tissues from R. aegyptiacus, qRT-PCR was performed with pooled cDNA from nose (nasal epithelium), trachea, lung, blood, spleen and duodenum. Threshold cycle (Ct) values were employed to determine the expression levels of the candidate reference genes. Their expression levels varied and EEF1A1 displayed the highest expression levels across different tissues as indicated by the lowest Ct values. The overall Ct values of EEF1A1, ACTB, GAPDH and PGK1 were 23 ± 1.5, 33 ± 3, 27 ± 2 and 27 ± 1.5, respectively (Fig. 1). EEF1A1 and PGK1 showed the lowest variability in Ct values, suggesting that these two genes display the most stable expression across diverse tissues from R. aegyptiacus.

Expression of candidate reference genes in primary bat fibroblasts upon IFN-I stimulation and incubation at various temperatures.
To investigate the expression stability of candidate reference genes, bat primary fibroblasts were incubated at 35 °C, 37 °C or 40 °C in the presence or absence of universal type I interferon (uIFN) for 4 h. We employed 35 °C, 37 °C and 40 °C to mimic the physiological daily oscillation of body temperature in R. aegyptiacus 9 . The expression levels of all candidate reference genes under these conditions were assessed by qRT-PCR. ACTB showed a broad variation in Ct values, ranging from 29.7 to 39.01 at 40 °C. The Ct values of ACTB were lower at 35 °C, 31.5 ± 0.7 compared to 33.5 ± 1 at 37 °C, suggesting an unstable expression of ACTB upon temperature changes (Fig. 2). The expression levels of EEF1A1, GAPDH or PGK1 were comparable at 35 °C, 37 °C and 40 °C, demonstrating stable expression of these potential reference genes under temperature oscillations. Since uIFN activates IFN-I pathway in bats 31 , we employed this cytokine to further investigate the expression of the candidate reference genes upon IFN-I stimulation. The mean Ct values of GAPDH decreased upon IFN-I stimulation at 35 °C, 37 °C and 40 °C, indicating increased expression of  (Table 2). NormFinder was employed as the second algorithm to determine gene stability, since it allows evaluation of the overall stability as well as the individual stability for each condition 35 . The most stable reference genes display stability values close to 0 according to this algorithm. Based on this method, we calculated the overall and individual stability values. The total stability values for ACTB, EEF1A1, GAPDH and PGK1 were 0.097, 0.013, 0.05 and 0.015, respectively, suggesting EEF1A1 as the most stable reference gene under all tested conditions. EEF1A1 also displayed the highest expression stability under high temperature or upon IFN-I treatment (Table 3). Altogether, both algorithms suggest that EEF1A1 is the reference gene with the highest expression stability under the selected conditions, followed by PGK1, GAPDH, whereas ACTB has the lowest stability.  ACTB expression was also stable at 35 °C and 37 °C in the presence of IFN-I. However, the expression level was significantly reduced following a 2 h IFN-I stimulation at 40 °C (Fig. 3C). Intriguingly, the expression of GAPDH was significantly increased in cells treated with IFN-I at all selected temperatures, suggesting that the induction of GAPDH by IFN-I is temperature independent in R. aegyptiacus ( Fig. 3A-C). To investigate whether GAPDH induction by IFN-I is specific to R. aegyptiacus, we stimulated human fibroblasts with IFN-I at 37 °C and 40 °C (Fig. 3D, E). Both human GAPDH and ACTB have been previously used as reference genes in several qRT-PCR studies 36,37 . Indeed, both genes displayed good expression stability at 37 °C, yet the variation of Ct values at 40 °C argued against stability of human GAPDH and ACTB at 40 °C. The relative expression of human GAPDH normalized against ACTB remained unchanged upon IFN-I treatment either at 37 °C or 40 °C, revealing that human GAPDH expression in fibroblasts is not modulated by IFN-I (Fig. 3D, E). Overall, we conclude that IFN-I triggers GAPDH expression specifically in R. aegyptiacus.

Discussion
Accurate gene transcription measurements require selection of reference genes that maintain high stability under various experimental conditions. In this study, we selected ACTB, GAPDH, EEF1A1 and PGK1 as reference gene candidates due to their wide applications in other species 14 Our finding that the expression of GAPDH, one of the key enzymes in glycolysis, is induced by IFN-I in R. aegyptiacus may have implications for the immunometabolism of bats. It is well established that metabolic reprogramming in cells controls their immune responses 53,54 . Glycolysis is commonly utilized by various immune cells to enable prompt responses to infections. In bone marrow derived macrophages, aerobic glycolysis promotes IL-1β production upon LPS stimulation or Bordetella pertussis infection 55 . In CD4 + and CD8 + T-cells 56-60 , as well   www.nature.com/scientificreports/ as in NK cells 61,62 , glycolysis is required for their effector functions, such as IFN-γ production and cytotoxicity. In human plasmacytoid dendritic cells (pDCs) and monocyte-derived DCs (moDC), glycolysis promotes IFN-I production upon TLR9 or RIG-I activation, respectively 63 . Further, IFN-I induces a metabolic shift towards glycolysis, contributing to the antiviral activity in fibroblasts and antigen presentation in DCs 64,65 . On the other hand, lactate, the end metabolite of glycolysis, directly binds to mitochondrial antiviral-signaling protein (MAVS), and consequently inhibits its activation and IFN-I production 66 . Thus, outcomes of such metabolic shift vary in diverse cells or under different stimulations. The nectarivore bat Glossophaga soricina employs high rates of glycolysis to generate ATP during flight 67 . As a fruit bat, R. aegyptiacus could also directly utilize dietary sugars to fuel both roosting and flight metabolism 68 . Whether induction of GAPDH expression in this species impacts on metabolic reprogramming towards glycolysis, needs to be clarified. Moreover, whether and how such metabolic shifts affect IFN-I signalling or other immune pathways in R. aegyptiacus remains to be uncovered. In addition to its roles in glycolysis, GAPDH also modulates cell death, RNA export and cytoskeleton dynamics 69 . Upon serum deprivation and DNA damage, GAPDH translocates to the mitochondria and interacts with voltage-dependent anion channel (VDAC), leading to apoptosis 70 . It can also bind to the 3' untranslated region of TNFα mRNA and represses TNFα expression in human monocytes and macrophages 71 . Hence, GAPDH upregulation by IFN-I may contribute to apoptosis induction and TNFα repression in R. aegyptiacus, which could represent novel mechanisms for the prevention of excessive inflammation during viral infections.
Overall, our study provides an extensive analysis of reference genes and identifies EEF1A1 as the most stable reference gene in R. aegyptiacus under temperature changes and IFN-I stimulation, which allows us to perform accurate gene transcription studies in this species. Our findings also open new investigation avenues by showing that GAPDH is regulated by IFN-I which has a broad relevance in context of immunometabolism.

Material and methods
Selection of reference gene candidates and design of primer pairs. Reference gene candidates (ACTB, EEF1A1, GAPDH and PGK1) for R. aegyptiacus gene expression studies were selected based on their utilization in other bat species 13,14,17 . Primers against these reference genes were designed using the PrimerQuest tool (Integrated DNA Technologies, Inc.). The criteria for primer design were as follows: primer lengths around 17-30 bp, GC content of 40-55%, optimal melting temperature at 62 °C, and amplicon lengths within a range of 100-250 bp. Derived primer pairs were evaluated using the OligoAnalyzer tool to exclude primers with hairpin structures and homo-and/or heterodimer formation (Integrated DNA Technologies, Inc.). Primer sequences were also blasted using the NCBI BLAST tool to ensure their specificity for R. aegyptiacus. The primer pairs meeting all criteria were selected for further experiments. The characteristics of these primers are shown in Table 4.   -GCC TTG GTC GTG GAT AAT G  R-GGG ATA CTT CAG GGT CAG GATA  193   EEF1A1  eukaryotic translation elongation factor  1 alpha 1  107509282  F-GTA TGC CTG GGT CTT GGA TAAA  R-GCC TGT GAT GTG CCT GTA A  162   GAPDH  glyceraldehyde-3-phosphate dehydrogenase  107519804  F-CAA GTT CAA AGG CAC AGT CAAG  R-TAT TCA GCA CCA GCA TCA CC  120   PGK1  phosphoglycerate kinase 1  107503843  F-GAT TAC CTT GCC TGT TGA CTTTG  R-GAC AGC CTC AGC ATA CTT CTT  Establishment of standard curves and examination of amplification efficiency via qRT-PCR. Standard curves of all reference genes in the qRT-PCR reaction were generated with copy numbers from 300,000 to 3 copies in ten-fold dilution steps. To achieve accurate copy numbers, amplicon sizes of each reference gene were used to calculate the specific weight of each amplicon 78 , resulting in 2.12 × 10 -19 g for ACTB, 1.78 × 10 -19 g for EEF1A1, 1.32 × 10 -19 g for GAPDH and 1.62 × 10 -19 g for PGK1. The amplification efficiency of all primer pairs was subsequently determined with the slope of the standard curve according to the equation 10 -1/slope -1. The amplification efficiency of favourable primers ranges between 90-110%. Linear dynamic range (LDR) is described as the highest to the lowest quantifiable copy numbers from standard curves. LDR should cover at least 3 orders of magnitude, ideally 5-6 orders. Precision refers to intra-assay variation and is defined as standard deviation (SD) of technical replicates 30 .

Cells, tissues and stimulation experiments.
Amplicon purification and sequencing. Amplicons of all candidate reference genes were visualized in 1.5% agarose gels and bands were cut and purified using the QIAquick gel extraction kit (#28506, Qiagen). Purified amplicons were subsequently sequenced using the Eurofins tubeseq platform. with Holm-Šidák's post-hoc test was performed. A P value of < 0.05 was considered to be significant.