Housekeeping Gene Expression in the Fetal and Neonatal Murine Thymus Following Coxsackievirus B4 Infection

The thymus fulfills the role of T-cell production and differentiation. Studying transcription factors and genes involved in T-cell differentiation and maturation during the fetal and neonatal periods is very important. Nevertheless, no studies to date have been interested in evaluating the expressions of housekeeping genes as internal controls to assess the varying expressions of different genes inside this tissue during that period or in the context of viral infection. Thus, we evaluated by real-time quantitative polymerase chain reaction (qPCR) the expression of the most common internal control genes in the thymus of Swiss albino mice during the fetal and neonatal period, and following in utero infection with Coxsackievirus B4. The stability of expression of these reference genes in different samples was investigated using the geNorm application. Results demonstrated that the expression stability varied greatly between genes. Oaz1 was found to have the highest stability in different stages of development, as well as following Coxsackievirus B4 infection. The current study clearly demonstrated that Oaz1, with very stable expression levels that outperformed other tested housekeeping genes, could be used as a reference gene in the thymus and thymic epithelial cells during development and following Coxsackievirus B4 infection.


Introduction
Coxsackievirus B (CV-B) group belongs to the human enterovirus B species, from the Enterovirus genus and Picornaviridae family. The virus particle is small, non-enveloped, and consists of a single-stranded RNA genome with positive polarity, encompassed in an icosahedral capsid. The genome comprises a single open-reading frame (ORF) flanked by two untranslated regions at both 5 and 3 ends, with a virus-encoded protein, VPg, linked to its 5 terminus and a polyadenylated tail at its 3 terminus. The ORF codes for a polyprotein that is proteolytically processed into 11 mature proteins, 4 of which are structural and 7 of which are non-structural (reviewed in Reference [1]). and streptomycin) was added and cells were washed by centrifugation at 500× g for 5 min. Cell pellets were resuspended and incubated for 10 min at room temperature with 20 µL of MagniSortTM Positive Selection Beads A. Up to 2.5 mL of cell separation buffer was added, and the tube was inserted into the 6-Tube Magnetic Separation Rack (Cell signal) and incubated for 5 min at room temperature. The supernatant containing CD45cells was collected and washed twice by centrifugation at 500× g for 5 min [12].

RNA Extraction
Total RNA was extracted using the NucleoSpin XS-kit (Macherey-Nagel, GmbH & Co, Düren, Germany) according to the manufacturer's instructions. Briefly, thymuses were lysed in 102 µL of RA1-TCEP mixture. Five microliters of carrier RNA working solution were added to the lysate. The sample was filtered through the NucleoSpin®Filter Column. One hundred microliters of 70% ethanol was added to the lysate and RNA bound to the NucleoSpin®RNA XS Column. The column was desalted using membrane-desalting buffer. DNA was digested using DNase reaction mixture. rDNase was inactivated with RA2 buffer and silica membrane washed twice with RA3 buffer. Pure RNA was eluted in 10 µL of nuclease-free water (Ambion). The RNA concentration and quality were assessed using a NanoDrop 2000 (UV-Vis Spectrophotometer; ThermoScientific, Waltham, MA, USA).

cDNA Synthesis
RNA was reverse-transcribed using the Transcriptor First Strand cDNA Synthesis Kit (Roche). Complementary DNA (cDNA) synthesis was first performed with approximately 500 ng of RNA in a final volume of 20 µL containing 20 U/µL of Transcriptor Reverse Transcriptase, 2.5 µM of anchored-oligo(dT)18 primer, 60 µM of random hexamer primer, 20 U of RNase inhibitor (all contained in the kit), and 1 mM of each deoxynucleoside triphosphate (dNTP, Promega).
According to the manufacturer's instructions, secondary structures were first denatured by heating the samples for 10 min at 65 • C in the presence of anchored-oligo(dT)18 primer, random hexamer primer, and nuclease-free water in 20 µL total volume. Tubes were then cooled on ice immediately, supplemented with buffer, protector RNase inhibitor, dNTP mix, and Transcriptor reverse transcriptase, and incubated for 10 min at 25 • C, followed by 30 min at 55 • C. The enzyme was inactivated by heating at 85 • C for 5 min, and the reaction was then stopped by placing tubes on ice. Two negative controls were performed in each reaction, one without RNA and the second without enzyme (RT minus control), to detect any genomic DNA contamination.

Quantitaive PCR with SYBR Green
qPCR was performed for seven genes using the iCycler iQ real-time detection system (Bio-Rad) using SYBR green detection. Primers (Eurogentec, Belgium) used for the PCR are listed in Table 1.
The cycling conditions were 1 cycle of denaturation at 95 • C for 10 min followed by 40 three segment cycles of amplification (95 • C for 30 s, 60 • C for 30 s, 72 • C for 25 s), where the fluorescence was automatically measured during PCR. A melting curve was done for each primer pair to check the presence of one gene-specific peak and the absence of primer dimers.
The endpoint used in the qPCR, Ct, was defined as the PCR cycle number that crossed an arbitrarily placed signal threshold. All samples were amplified in triplicate and the means were used for further analysis. In the same way, two controls were also included in each PCR reaction, a reverse-transcription negative control (without reverse transcriptase) to confirm the absence of genomic DNA, and a non-template negative control to check for primer dimers.

Viral Load
One microliter of cDNA was used to quantify CV-B4 RNA by qPCR with the same program as described for HKGs. The primer sequences used were: forward primer EV1 5 -CAAGCACTTCT GTTTCCCCGG-3 and reverse primer EV2 5 -ATTGTCACCATAAGCAGCCA-3 . Enterovirus 71 RNA control (Vircell, Kindly provided by Pr Didier Hober, Université Lille II) was used as a standard for a calibration curve containing five points ranging from 1.26 × 10 4 to 1.26 copies/mL.

Analysis of Expression Stability
The amount of internal reference gene relative to a calibrator (fold change between two Ct values) has been expressed as fold difference, which was calculated as [20] Fold difference = 2−∆Ct (1) For each gene, the coefficient of variation (CV) of the cycle threshold (Ct) value was calculated. The CV equals the SD divided by the mean of Ct (expressed as a percentage). The MFC was calculated to reflect the minor variation in expression of those candidate HKGs within the large dataset. The MFC was defined as a ratio of the maximum and minimum values observed within the dataset.
GeNorm is the best algorithm used to determine the most stable HKG from a set of tested candidate genes in a given sample subjected to different experimental conditions [21]. geNorm was originally developed by Vandesompele et al. [21] and is implemented in qbase+ software (version 3.1, Ghent, Belgium). The geNorm stability M value can be determined by this software, in which every HKG is tested against the other genes in a pairwise variation that serially excludes the least stable genes from the analysis.
A candidate HKG must have a small CV, MFC <2 [22], and the lowest M value [21]. The ranking of HKG was essentially based on the M value of geNorm (from lowest to highest), defined as the average pairwise variation for that gene with all the other tested control genes.

Statistical Analysis
Statistical analysis was performed with GraphPad Prism 5 software. The unpaired t-test was used for fold change analysis in CV-B4-versus mock-infected samples. A p value < 0.05 was considered statistically significant. * p value < 0.05, ** p value < 0.01, *** p value < 0.001, **** p value < 0.0001. Pearson's correlation test was used to analyze the correlation between the viral load and the fold change of a given HKG at a given time point. Pearson's r > 0: positive correlation; Pearson's r < 0: negative correlation.

Evaluation of RNA Purity, Integrity, and Contamination
High levels of RNA purity were obtained, as confirmed by A260/A280 ratios in the range of 1.8 to 2.0. RNA integrity was defined as an RNA Integrity Number (RIN) >8 by Agilent 2100 Bioanalyzer (Agilent technologies), and was revealed to be in the acceptable range for all samples. A single amplicon, represented by a generation of a single peak, for each candidate gene was confirmed by qPCR melting curve analysis. DNA contamination for all samples was ruled out by negative results (cycle threshold, Ct > 40) for RT (enzyme) minus control.

Housekeeping Gene Expression in the Developing Thymus
Seven HKGs commonly used in quantitative gene expression to normalize Ct values were used in this study, namely Hprt, Actb, Sdha, RpS29, Rpl4, Oaz1, and Gapdh, the transcripts of which were quantified using RT-qPCR.
Candidate reference genes showed Ct values between 15 and 27. The variability of selected HKGs in the thymus has been presented herein as standard deviations (SD) and ranges, expressed as an average fold change from the mean and as a maximum fold changes (MFC), respectively.
As shown in Figure 1, Oaz1 and Hprt were the most stable HKGs at the three tested time points. For Oaz1, the fold change was between 0.6 and 1.58 (average fold change = 1.033). For Hprt, the fold change was between 0.38 and 2.83 (average fold change = 1.087). Rpl4 varied moderately, 2.5-fold at Day 17G and Day 5, and 5.3-fold at Day 1. There was considerably greater variability for Actb (up to 6-fold at Day 1), Sdha (up to 9-fold at Day 1), RpS29 (up to 8-fold at Day 17G) and Gapdh (up to 7-fold at Day 5). Thus, we found Oaz1 and Hprt to be the best combination of reference genes for studies of transcript levels in the developing murine thymus.
Genes 2020, 11, x FOR PEER REVIEW 6 of 15 considered statistically significant. *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001. Pearson's correlation test was used to analyze the correlation between the viral load and the fold change of a given HKG at a given time point. Pearson's r > 0: positive correlation; Pearson's r < 0: negative correlation.

Evaluation of RNA Purity, Integrity, and Contamination
High levels of RNA purity were obtained, as confirmed by A260/A280 ratios in the range of 1.8 to 2.0. RNA integrity was defined as an RNA Integrity Number (RIN) >8 by Agilent 2100 Bioanalyzer (Agilent technologies), and was revealed to be in the acceptable range for all samples. A single amplicon, represented by a generation of a single peak, for each candidate gene was confirmed by qPCR melting curve analysis. DNA contamination for all samples was ruled out by negative results (cycle threshold, Ct > 40) for RT (enzyme) minus control.

Housekeeping Gene Expression in the Developing Thymus
Seven HKGs commonly used in quantitative gene expression to normalize Ct values were used in this study, namely Hprt, Actb, Sdha, RpS29, Rpl4, Oaz1, and Gapdh, the transcripts of which were quantified using RT-qPCR.
Candidate reference genes showed Ct values between 15 and 27. The variability of selected HKGs in the thymus has been presented herein as standard deviations (SD) and ranges, expressed as an average fold change from the mean and as a maximum fold changes (MFC), respectively.
As shown in Figure 1, Oaz1 and Hprt were the most stable HKGs at the three tested time points. For Oaz1, the fold change was between 0.6 and 1.58 (average fold change = 1.033). For Hprt, the fold change was between 0.38 and 2.83 (average fold change = 1.087). Rpl4 varied moderately, 2.5-fold at Day 17G and Day 5, and 5.3-fold at Day 1. There was considerably greater variability for Actb (up to 6-fold at Day 1), Sdha (up to 9-fold at Day 1), RpS29 (up to 8-fold at Day 17G) and Gapdh (up to 7-fold at Day 5). Thus, we found Oaz1 and Hprt to be the best combination of reference genes for studies of transcript levels in the developing murine thymus.  Figure 1. Fold change in HKG expression in the developing mouse thymus. mRNA was extracted and converted to cDNA. The cDNA was subjected to real-time quantitative PCR using gene-specific primers for seven selected HKGs. Expression variability of each HKG in thymuses (n = 7-12), sampled from fetuses at Day 17 of gestation (Day 17G) (A), and newborns at Day 1 (B) and Day 5 (C) from birth, is represented by log-scale of fold change (lines represent a median). The fold change represented was calculated using the equation presented by Livak and Schmittgen, 2001 [20]: Fold difference = 2−∆Ct.

Housekeeping Genes Following CV-B4 Infection
The presence of CV-B4 E2 RNA in the thymus was checked at different p.i. times by RT-qPCR using specific primers of the highly conserved 5 untranslated region of enterovirus genome, and Fold changes in HKG expression in CV-B4-infected compared to mock-infected thymuses, at different p.i. times, are represented in Figure 2.
The presence of CV-B4 E2 RNA in the thymus was checked at different p.i. times by RT-qPCR using specific primers of the highly conserved 5′ untranslated region of enterovirus genome, and infection was demonstrated in different samples. All samples retained for analysis of the effect of the infection on HKG expression were positive for CV-B4 RNA.
Oaz1 expression remained fairly stable from Day 17G to Day 5, even in infected thymuses, with fold changes between 0.5 and 1.7 following infection at Day 10G or 17G, as in mock-infected thymuses. Hprt seemed to be a little less stable, with fold changes between 0.37 and 2.92. A significant decrease in Actb was noted at Day 1 following CV-B4 inoculation at Day 17G (fold change between 0.45 and 2.2 versus 0.58 and 6.17 in mock-infected thymuses, p = 0.022). For Sdha, a significant variation was noted at Day 1 (fold change from 0.21 to 6.73 versus 0.18 to 17.75 in mock, p = 0.015) and Day 5 (fold change from 0.29 to 4.06 versus 0.54 to 10.72 in mock, p = 0.041) following CV-B4 inoculation at Day 10G. For Rpl4, a significant variation in expression was observed at Day 1 and Day 5 following CV-B4 inoculation at Day 17G (fold change from 0.369 to 2.71 versus 0.29 to 5.55 in mock, p = 0.012 and from 0.63 to 1.55 versus 0.42 to 2.58 in mock, p = 0.0098, respectively). A significant variation in Gapdh expression following CV-B4 inoculation at Day 17G was equally observed at Day Variances in CV-B4-versus mock-infected samples were compared using the unpaired t-test. A p value < 0.05 was considered statistically significant. *p value < 0.05, ** p value < 0.01, *** p value < 0.001, **** p value < 0.0001.
Oaz1 expression remained fairly stable from Day 17G to Day 5, even in infected thymuses, with fold changes between 0.5 and 1.7 following infection at Day 10G or 17G, as in mock-infected thymuses. Hprt seemed to be a little less stable, with fold changes between 0.37 and 2.92. A significant decrease in Actb was noted at Day 1 following CV-B4 inoculation at Day 17G (fold change between 0.45 and 2.2 versus 0.58 and 6.17 in mock-infected thymuses, p = 0.022). For Sdha, a significant variation was noted at Day 1 (fold change from 0.21 to 6.73 versus 0.18 to 17.75 in mock, p = 0.015) and Day 5 (fold change from 0.29 to 4.06 versus 0.54 to 10.72 in mock, p = 0.041) following CV-B4 inoculation at Day 10G. For Rpl4, a significant variation in expression was observed at Day 1 and Day 5 following CV-B4 inoculation at Day 17G (fold change from 0.369 to 2.71 versus 0.29 to 5.55 in mock, p = 0.012 and from 0.63 to 1.55 versus 0.42 to 2.58 in mock, p = 0.0098, respectively). A significant variation in Gapdh expression following CV-B4 inoculation at Day 17G was equally observed at Day 1 and Day 5 (fold change from 0.21 to 16.56 versus 0.42 to 1.67 in mock, p < 0.0001 and from 0.14 to 2.51 versus 0.14 to 7.93 in mock, p = 0.0036, respectively), together with a significant variation at Day 5 following CV-B4 inoculation at Day 10G (fold change from 0.51 to 2.36 versus 0.14 to 7.93 in mock, p = 0.0008). RpS29 gene was the most variable gene, with variation at Day 1 (fold change between 0.06 and 22.45 versus 0.5 and 3.48 in mock, p = 0.0001) consequent to virus inoculation at Day 17G, and significant variations at Day 5 following CV-B4 inoculation at either Day 10G or 17G (fold change from 0.42 to 2.93 and from 0.25 to 2.28 versus 0.13 to 5.39 in mock, p = 0.017 and p = 0.0026, respectively).
To better assess the effect of CV-B4 infection on HKG expression, the intensity of the infection was evaluated through measurement of the viral RNA copy number in the different samples. As illustrated in Figure 3A, viral loads could be segregated into two levels: low loads (around 8000 viral RNA copies/100ng of total RNA at Day 17G, following CV-B4 inoculation at day 10G, as well as at Days 1 and 5, following inoculation at Day 17G) and high loads (around 50,000 viral RNA copies/100 ng of total RNA at Days 1 and 5, following inoculation at Day 10G). Thus, at a given time point (mainly Day 1 or Day 5, because of the above-mentioned effects of the infection at those time points, as well as the marked differences in viral loads), we could appreciate an eventual correlation between the viral load and fold changes of the different HKG. As illustrated in Figure 3B, there were only few inverse correlations (Pearson's r < 0) between the viral load and the fold change of a given HKG. In thymuses sampled the day of birth, this was the case for Actb and Sdha following CV-B4 inoculation at Day 10G (r = −0.8785 and r = −0.7226 respectively), and for Gapdh following CV-B4 inoculation at Day 17G (r = −0.8215). At day 5 from birth, only a negative correlation between the viral load and the fold change of Actb expression was observed in thymuses from mice inoculated at Day 10G (r = −0.7366). There was no correlation when, at a given time point, analyzing together samples from mice inoculated at Day 10G (high viral loads) and those inoculated at Day 17G (low viral loads). Any evident correlation could be found in thymuses sampled at Day 17G (data not illustrated).

Oaz1 as a Candidate Gene
A candidate HKG is defined as a gene with minimal variability in a tissue under different experimental conditions. It must have a small coefficient of variation (CV) and a MFC < 2.
All these parameters were calculated (Table 2) in all samples together, regardless of infection conditions and developmental stages. Tested HKGs were ranked from the most to the least stable according to different parameters (Mean Ct, SD, CV, MFC, and M value).

Oaz1 as a Candidate Gene
A candidate HKG is defined as a gene with minimal variability in a tissue under different experimental conditions. It must have a small coefficient of variation (CV) and a MFC < 2.
All these parameters were calculated (Table 2) in all samples together, regardless of infection conditions and developmental stages. Tested HKGs were ranked from the most to the least stable according to different parameters (Mean Ct, SD, CV, MFC, and M value). Taken together, all values within the dataset showed that Oaz1 was the most stable among all seven tested HKGs, with the lowest CV, MFC, and M value (M = 1.665, calculated by geNorm).
Similarly, for TEC analysis (Table 3), Oaz1 and Hprt had the lowest CV, MFC, and M values (1.308 for both genes), followed by Gapdh and Sdha which varied moderately with M value, equal to 1.693 and 1.839 respectively. For Rpl4, Actb, and RpS29, there was an extreme variation in expression, with the biggest CVs and MFCs and with an M value > 2. Table 3. Ranking of the seven tested housekeeping genes in isolated TECs, from the most to the least stable.

Discussion
Real-time qPCR is a rapid, sensitive, and reliable method for the detection and quantification of mRNA transcription level. The principle of real-time PCR quantification is based on two concepts with different bases: absolute and relative quantification [7]. Absolute quantification determines the input copy number based on standard dilutions of known concentration. In relative or comparative quantification, the concentration of the target gene is expressed in comparison to a reference gene in the same sample.
HKGs are used as internal standards to normalize mRNA levels between different samples for a reliable comparison of mRNA transcription level. Although many studies have used those classic constitutive genes as presumed stable references, their expression may vary considerably depending on tissues, developmental stages, and specific experimental conditions [23]. It is therefore important to use as an internal control an HKG of which the expression remains stable during the studied growth and experimental conditions. Until now, no validated reference genes have been reported for the normalization of transcript levels in the thymus during development or following CV-B infection. It is obvious that selection of a reference gene of high quality is of crucial importance for normalization, and must be identified with increasingly stringent standards. Thus, our aim was to identify an HKG with minimal variability under our different experimental conditions.
For this purpose, a set of seven genes frequently used and validated (during development or under different experimental conditions in study conducted on humans or mice) as controls were chosen in our study, namely, Hprt, Actb, Sdha, RpS29, Rpl4, Oaz1, and Gapdh. We purposely chose genes with different functions (Table 1) in order to avoid genes belonging to the same biological pathway and that might be co-regulated.
First of all, gene expression was evaluated in thymuses (minimum n = 7) sampled from Swiss albino mice fetuses at Day 17G and neonates at Day 1 and Day 5. cDNA was obtained by reverse transcription of extracted RNA, and then quantified by qPCR.
Given the small size of the thymus during the fetal and neonatal stages, in order to obtain a pure RNA extract, total RNA was extracted by the NucleoSpin XS-kit since it seemed to be the most effective method for the extraction of high-quality RNA [24]. RNA integrity and purity was verified by RIN and the ratio A260 nm/A280 nm absorbance. The absence of DNA contamination was checked by a single peak qPCR analysis, Ct > 40 for RT minus control. Obtaining a unique qPCR product on gel electrophoresis confirmed the absence of primer dimers.
To determine the most stable gene, a combination between CV, MFC, and M value (calculated by geNorm) parameters was used. Our study demonstrated that the expression stability varied greatly between genes, and the expressions of six out of the seven investigated genes were found to be highly variable.
Variation in the expression of some HKGs during cell proliferation and differentiation has already been described. Furthermore, Piechaczyk et al. [25] found that Gapdh mRNA abundance varies from rat tissue to another. Similarly, an elevated expression of Hprt mRNA in the central nervous system, compared to other tissues, has been reported [26].
The combination of RpS29, Rpl4, and Oaz1 was used in a study conducted by Taves et al. [27] to evaluate genes implicated in the glucocorticoid metabolic pathway in lymphoid organs, without addressing the issue of expression stability. In contrast to our study, the evaluation of gene expression using those three HKGs was in the thymus of neonatal (6 days old) and adult C57BL/6 mice.
Among one subset of 535 genes expressed in 11 major fetal and adult human tissues, only 47 transcripts were expressed at the same level in fetal and adult tissues [28]. An increase in Gapdh expression was also noted during fetal rat brain development [29].
Our current results showed that Oaz1 was the most stable among all seven tested candidate genes in fetal and neonatal Swiss albino mice thymus.
Recently, Medrano and associates [13] demonstrated that, in the thymus, among six tested genes, Tbp (TATA box binding protein) gene was the most stable gene during the period, versus Hprt1 during the post-natal period (Day 21). Contrary to our current study, that investigation was performed on inbred mice (C57BL/6J strain) and only in the context of development, without addressing the infectious context.
Compared to the thymus, HKG expression in isolated TECs did not present a big difference. Oaz1, Hprt, then Gapdh had the smallest variations and RpS29 was the least stable gene. Similarly, in a study conducted by St-Pierre et al. [30], whole-transcriptome sequencing revealed that Oaz1, Hprt, and Gapdh presented the same transcript abundance in cTECs and mTECs isolated from the thymus of 7 day old C57BL/6 mice.
Kermani and colleagues [31] evaluated the expression of the growth-hormone-and insulin-related antigens in Balb/c thymus during the embryonic and neonatal period. Although they chose Hprt as an internal gene in their study, the question of instability of that reference gene was not addressed. In the current study, we worked with a different mouse strain, namely Swiss albino mice, and in an infectious context.
Despite the scarcity of investigations, instability in HKG expression following an infection has been previously demonstrated. Furthermore, an increase in ribosomal RNA (rRNA) expression was reported following cytomegalovirus (CMV) infection of human fibroblast cells [32]. An extreme variation in Actb and Gapdh was found during Herpes Simplex Virus-1 (HSV-1) infection of human embryonic lung fibroblasts [33]. In a study conducted by Watson and co-workers [34], among 10 tested reference genes, only 3 (Peptidylprolyl isomerase A, Gapdh, and Sdha) showed any variation in cell cultures infected by human immunodeficiency virus type 1 (HIV-1), HSV-1, CMV, and varicella zoster virus (VZV). In the current study, CV-B4 infection had a significant effect on the expression of Gapdh and RpS29, then Actb and Rp14, then Sdha. These genes were affected by either increase and/or decrease of their expression, observed at Day 1 and/or Day 5 from birth, and mainly following CV-B4 inoculation at Day 17G but also following infection at Day 10G. Plotting fold changes with the corresponding viral loads in an attempt to better assess the effect of the infection revealed only a few negative (inverse) correlations among samples showing a significant effect of the virus (one with Sdha and the other with Gapdh). A negative correlation was also observed between the viral load following inoculation at Day 10G and the corresponding fold change of Actb (at both Day 1 and Day 5 from birth), although a significant effect of the virus on Actb expression could be observed only following inoculation at Day 17G. A more important variation was observed following mouse inoculation at Day 10G. This pronounced effect may logically be attributed to the higher viral load in thymuses sampled following an inoculation at Day 10G than in those sampled following inoculation at Day 17G.
During development, a gene's expression is subject to multiple factors that cause it to vary. In the context of viral infection, intracellular replication and protein expression of CV-B may interfere with cellular physiology. In general, the virus diverts cellular machinery in its favor. Ribosomal genes are among those that may be affected. Rps20 gene expression was noted in a previous study following CV-B3 infection of rat cardiac myocytes [35]. Nevertheless, no previous study has demonstrated an effect of CV-B on such genes. As described in the introduction, CV-B preferentially infects the pancreas, heart, and central nervous system, and thus, studies have generally focused on the effect of the virus on the expression of specific genes of those tissues or of genes involved in the immune response.
Other conditions and treatments can also affect HKG expression, which necessitates their validation as internal reference genes before use. Indeed, HKG expression has been reported to vary considerably in leukemic patients' tumor samples [36] and in glioblastoma samples [37]. Additionally, the expression of Gapdh mRNA increases in case of insulin stimulation [38]. Furthermore, Burger and co-workers [39] demonstrated that chemotherapeutic drugs inhibit rRNA synthesis at different levels.
According to our knowledge, this work is the first validation of HKG in the thymus and isolated TECs of Swiss albino mice following in utero infection with Coxsackievirus B4, and provides useful data supporting further gene analysis in thymus mice during the fetal and neonatal period.

Conclusions
In conclusion, based on both expression stability and expression level, our data suggest that Oaz1 is a good reference gene for evaluation and normalization of gene expression in the thymus during ontogeny and after CV-B4 infection. To our best knowledge, this was the first study to assess the expression stability of candidate reference genes in the thymus during fetal development (Day 17G), in neonatal life (Day 1 and 5), and after in utero viral infection.