Validation of reference genes for expression analysis in a murine trauma model combining traumatic brain injury and femoral fracture

Systemic and local posttraumatic responses are often monitored on mRNA expression level using quantitative real-time PCR (qRT-PCR), which requires normalisation to adjust for confounding sources of variability. Normalisation requests reference (housekeeping) genes stable throughout time and divergent experimental conditions in the tissue of interest, which are crucial for a reliable and reproducible gene expression analysis. Although previous animal studies analysed reference genes following isolated trauma, this multiple-trauma gene expression analysis provides a notable study analysing reference genes in primarily affected (i.e. bone/fracture callus and hypothalamus) and secondarily affected organs (i.e. white adipose tissue, liver, muscle and spleen), following experimental long bone fracture and traumatic brain injury. We considered tissue-specific and commonly used top-ranked reference candidates from different functional groups that were evaluated applying the established expression stability analysis tools NormFinder, GeNorm, BestKeeper and RefFinder. In conclusion, reference gene expression in primary organs is highly time point as well as tissue-specific, and therefore requires careful evaluation for qRT-PCR analysis. Furthermore, the general application of Ppia, particularly in combination with a second reference gene, is strongly recommended for the analysis of systemic effects in the case of indirect trauma affecting secondary organs through local and systemic pathophysiological responses.

Although pre-hospital and in-hospital care has profoundly advanced over the past decades 1,2 , major trauma remains a leading cause of death and disability, mainly in younger ages 3,4 . Treatment of patients with multiple injuries presents a medical challenge, as the complex posttraumatic response is injury-unique and includes the interaction between different organ systems [5][6][7] . Despite strenuous clinical and experimental efforts [8][9][10][11] , the underlying pathophysiological and immunological responses triggered by the combination of multiple injuries still remain to be elucidated.
Systemic and local posttraumatic gene response is commonly monitored on mRNA expression level [12][13][14][15] , using the cost-efficient, easy to adapt and widely established method of quantitative real-time reverse transcriptionpolymerase chain reaction (qRT-PCR). Despite its advantages, the results are influenced by confounding sources of variability during mRNA extraction, quantification and cDNA transcription. Thus, internal normalisation through a reference gene is required for accurate data analysis 16 . The ideal reference or housekeeping gene for normalisation is stable in the tissue of interest, throughout time and under all experimental conditions. However, previous studies demonstrated that commonly used reference genes such as Gapdh are indeed differentially regulated among distinct cell types and organ systems [17][18][19] . Therefore, the time-consuming identification of suitable housekeeping genes for the qRT-PCR data normalisation remains essential for each newly designed study, Bone/fracture callus analysis. The expression analysis in primarily traumatised organs is challenging, as the gene expression pattern fundamentally changes. In our fracture model, we aimed to compare the intact bone with fracture callus tissue, both consisting of highly different cell composition. As the callus in the applied model offers low sample volume, only five predominantly applied reference genes were tested.
At d3 after trauma in bone/callus, we observed a strong expression (low Ct) of all reference genes tested (Fig. 1A). While NormFinder ranked Hprt (sv = 0.225) first, Gapdh (sv = 0.245) second, and Actb (sv = 0.517) last (Fig. 1C), GeNorm identified the highest stability for B2m and Ppia (M = 0.509), but also Actb (M = 0.916) as the least stable. As none of the candidate reference genes reached the threshold for homogeneous sample panels (M < 0.5) during the GeNorm analysis, the standard threshold for heterogeneous samples (M < 1.5) was applied 27,28 . GeNorm's pairwise variation revealed the need for the combination of three reference genes (V 3/4 = 0.15) (Fig. 1B). The BestKeeper algorithm calculated a high coefficient of correlation for Hprt (r = 0.928) Table 1. qPCR primer set (murine) of all reference genes tested, contrasting the top-ranked generally applied (bold = applied in all tissues) and the tissue-specific (unbold = applied in selected tissues) candidates.

Scientific RepoRtS
| (2020) 10:15057 | https://doi.org/10.1038/s41598-020-71895-x www.nature.com/scientificreports/ and Gapdh (r = 0.866) and low standard deviations for Gapdh (SD = 0.53) and B2m (SD = 0.57). Our final ranking, as a summation of the single ranks, determined Gapdh, Hprt and B2m as the most suitable candidates for a combined bone and fracture callus analysis at three days after trauma (Fig. 1C). This combination was also verified by RefFinder (Fig. 1D). An additional inter-and intra-group variation analysis identified a strong dysregulation of Actb following trauma (Fig. 1E), which was confirmed by the separate evaluation of reference genes in the intact bone or in the fracture callus using RefFinder. In both tissues bone and callus, Gapdh, Hprt and B2m remained the most stable reference genes, while Actb showed the strongest variability (Fig. 1F,G). At d7 after trauma in bone/callus, we conducted the same evaluation as for d3, as the healing process provokes a great gene expression modulation in the fracture callus, as well as in intact bone ( Fig. 2A-D). Here NormFinder ranked Hprt first and Ppia second (Fig. 2C), whereas GeNorm calculated the highest stability for Actb and Ppia. The GeNorm analysis additionally revealed that the combination of two reference genes is sufficient for housekeeping in bone and callus tissue at d7 after trauma (Fig. 2B). BestKeeper calculated the highest coefficient of correlation for Ppia and the lowest SD for Hprt. Our final ranking and RefFinder confirmed Ppia and Hprt as the most suitable reference genes in bone and fracture callus seven days after trauma (Fig. 2C,D).

Figure 1.
Reference gene expression stability analysis in bone/callus tissue at 3 days following surgery. (A) Distribution of threshold cycle (Ct) values for each group (Fx = fracture; TBI = traumatic brain injury; TBI + Fx = combined injury) of all candidate reference genes tested. (B) Pairwise variation (V n/n+1 ) analysis by GeNorm, to determine the optimal number of reference genes for normalisation. (C) Stability analysis by NormFinder (stability value: sv), GeNorm (M value) and BestKeeper (coefficient of correlation: r, standard deviation: SD), with a sum displayed in the final rank (green = most stable, red = least stable). (D) Overall expression stability confirmation of RefFinder, which combines NormFinder, GeNorm, BestKeeper and the comparative delta-Ct method. (E) inter-and intra-group variation calculated with NormFinder. (F/G) expression stability confirmation of RefFinder for intact bone (F) or callus tissue (G).  Fig. S1 online). Like the fracture callus, the hypothalamus offered a low sample volume, which again allowed the testing of only the five predominantly applied reference genes. At d3 after trauma in the hypothalamus, we observed a strong expression (low Ct) of all reference genes tested (Fig. 3A). NormFinder and GeNorm identified Actb (sv = 0.125, M = 0.243) and Hprt (sv = 0.148, M = 0.243) as the most stable reference genes (Fig. 3C). The GeNorm algorithm revealed that the combination of two genes (V 2/3 = 0.111) is sufficient as a reference in the hypothalamus at d3 after trauma (Fig. 3B). BestKeeper calculated the highest coefficient of correlation for Actb (r = 0.780) and the lowest SD for Hprt (SD = 0.15). The results of NormFinder, GeNorm, BestKeeper and the final rank identified Actb and Hprt as the most suitable combination of housekeepers (Fig. 3C), which was confirmed by RefFinder (Fig. 3D). All algorithms defined B2m as the least suitable reference gene, which exclusively missed the GeNorm threshold. An additional inter-and intra-group variation analysis revealed a strong regulation of B2m in the injured hypothalamus in contrast to the intact hypothalamus (Fig. 3E). Although B2m is highly regulated in the injured hypothalamus (Fig. 3G), it remained the most stable housekeeper in the intact hypothalamus ( Fig. 3F).
At d7 after trauma in the hypothalamus, we observed a strong expression (low Ct) of all reference genes tested (Fig. 4A). NormFinder determined the best stability for Gapdh, while GeNorm identified Hprt and Ppia as the most adequate reference gene combination (Fig. 4C). During the GeNorm analysis Actb and B2m failed the suitability threshold. GeNorm revealed that the combination of two reference genes is sufficient in the hypothalamus at d7 (Fig. 4B). BestKeeper calculated the lowest SD for Gapdh and Hprt. Our final ranking and RefFinder confirmed Gapdh and Hprt as the most suitable reference genes in hypothalamus at d7 after trauma (Fig. 4C,D).
White adipose tissue (WAT) analysis. As metabolism is strongly regulated after trauma, we investigated the stability of recommended reference genes in WAT at d3. Here gene expression was highest (low Ct) for Ppia and lowest for Tbp (Fig. 5A). NormFinder ranked Psmb2 (sv = 0.084) the best target gene and Ppia (sv = 0.105) second. Highest variation was observed for Gapdh (sv = 0.286) (Fig. 5C). GeNorm revealed the best M value for Ppia and Tbp (M = 0.234), as well as a highly variable regulation of Gapdh (M = 0.453). During the GeNorm analysis all investigated housekeepers reached the threshold of suitability of M < 0.5, while showing a sufficiency for the application of two reference genes (V 2/3 < 0.15) in WAT d3 after trauma (Fig. 5B). The BestKeeper algorithm determined the highest correlation for Psmb2 (r = 0.988) and Ppia (r = 0.983). In contrast, Gapdh was www.nature.com/scientificreports/ the least eligible housekeeping gene tested with a low coefficient of correlation and high SD (Fig. 5C). Overall, our final ranking identified Ppia and Psmb2 as the most stable reference genes in WAT three days after trauma (Fig. 5C). RefFinder verified Ppia and Psmb2 as the best combination, with Gapdh showing again the strongest differential regulation (Fig. 5D).
To prove that the suggested combination of Ppia and Psmb2 from d3 remains suitable even at seven days after trauma, we conducted the same evaluation for d7 ( Fig. 5E-H). Ppia and Psmb2 showed the best outcome among all algorithms applied, while Gapdh remained the most variable. In GeNorm all tested genes reached the threshold, and the algorithm revealed that two housekeeping genes are sufficient for a subsequent gene expression analysis (Fig. 5F). Our final rank and RefFinder confirmed Ppia and Psmb2 as the best reference genes in white adipose tissue at seven days after trauma (Fig. 5H).
Liver analysis. Another secondary organ analysed was the liver, as altered liver metabolism is a core element of posttraumatic response 23,24 . In liver tissue at d3 after trauma high gene expression (low Ct) of Ppia and B2m was observed, while Tbp showed the lowest expression (Fig. 6A). NormFinder revealed the highest stability for Ppia (sv = 0.157) and Hprt (sv = 0.199), as well as great variability of Actb (sv = 0.395). GeNorm calculated the www.nature.com/scientificreports/ best stability for Hprt and B2m (M = 0.261), followed by Ppia (Fig. 6C). Among all housekeeping genes tested, Actb (M = 0.675), Gapdh (M = 0.588) and Hmbs (M = 0.544) did not reach the threshold of suitability during the GeNorm analysis. GeNorm also revealed that the combination of two genes (V 2/3 = 0.134) is sufficient for housekeeping in the liver tissue at d3 (Fig. 6B). BestKeeper identified the best coefficient of correlation for Ppia (r = 0.996) and Hprt (r = 0.981) with a low SD. Again, our final rank as well as RefFinder confirmed Ppia and Hprt as the best reference gene combination, while B2m was ranked third at d3 after trauma in the liver (Fig. 6C,D). At d7 after trauma in liver tissue ( Fig. 6E-H), NormFinder revealed a low stability value for all investigated genes, ranking B2m first, Hmbs second and Ppia third (Fig. 6G). GeNorm identified Ppia and Hmbs as the best housekeeping gene combination. Although BestKeeper calculated the highest coefficient of correlation for B2m and Hmbs, Ppia still offered a high coefficient of correlation with the lowest SD. Our final ranking and RefFinder both confirmed B2m, Ppia and Hmbs as the best reference gene combination in liver at d7 (Fig. 6G,H).
Besides the secondary organs WAT and liver, our candidate reference genes were additionally evaluated in the metabolic-and immuno-active tissues muscle and spleen (see Supplementary Fig. S2 online). In both tissues, Ppia was once again identified as the most stable reference gene.

Target gene analysis.
To illustrate the effect of the reference gene analysis conducted, we normalised trauma relevant target genes ( Table 2) [31][32][33] to strongly regulated as well as to the most suitable reference genes identified (Fig. 7A-D). In bone/callus (d3), the normalisation of the target gene Atf4 to Actb resulted in a large variance with no informative value (Fig. 7A), whereas the use of B2m, Gapdh and Hprt disclosed a significant lower expression of Atf4 in the fracture callus (groups: Fx and TBI + Fx) compared to the intact bone (groups: control and TBI). As B2m expression is strongly modulated in the injured hypothalamus, an artificial Creb1 regulation was observed when applying B2m for housekeeping in the hypothalamus at d3. However, using the recommended reference gene combination of Hprt and Actb, a constant Creb1 expression across all groups was revealed (Fig. 7B). In WAT (d3), Creb1 seemed to be significantly regulated following fracture, when normalised to Gapdh. However, when applying the most suitable reference gene combination Ppia and Psmb2, the regulation was no longer valid (Fig. 7C). In liver (d3), ApoE showed no distinct regulation, after being normalised to the most unstable candidate Actb. When using the most suitable reference gene combination of Ppia and Hprt however, a significant regulation was revealed (Fig. 7D).
In summary, these examples strongly underline that the application of suitable reference genes reduces false positive results (Fig. 7B,C) and strengthens relevant regulations in gene expression (Fig. 7A,D).

Discussion
The choice of adequate reference genes for qRT-PCR normalisation is essential to detect target gene regulation, especially if the expression alteration is smaller than two fold 34,35 . Our analysis provides recommendations for the use of reference genes in a murine multiple-trauma model (Table 3). In primarily traumatised organs, housekeeping appears to be much more complicated than in secondary organs, particularly when comparing injured and intact tissues. Therefore, an individual evaluation for each condition and time point is required, which was also concluded by previous studies focusing on isolated traumatic brain injury 20,36 . For the comparison of gene expression in the intact bone and fracture callus we suggest the use of Gapdh, Hprt and B2m at three days (d3) as well as Ppia and Hprt at seven days (d7) following trauma. Although Gapdh 37,38 and Actb 39 were often used for housekeeping in bone studies, we recommend at least two reference genes in order to allow adequate normalisation for such a complex and multicellular organ under the challenging circumstances of trauma 34,40 . As in the hypothalamus the results of the reference gene analysis differ between d3 and d7 after trauma, no general guidance for this primary organ can be provided. Therefore, analogue to the bone/callus, we suggest an individual evaluation of the time point analysed in primarily traumatised organs to identify the most suitable housekeeping. At both time points Hprt, which has previously been recommended for isolated TBI and ischemic  22,41 , showed a high expression stability in the hypothalamus. Hprt, a hypoxanthine-phosphoribosyltransferase (Table 1), is essential for the synthesis of purine nucleotides 42 . The mRNA expression of Hprt seems to stay unaffected by the induction of TBI. In contrast, the expression of B2m was strongly regulated following TBI. Beta-2-microglobulin, part of the major histocompatibility complex (MHC), is directly involved in the immune response of the injured brain and therefore inappropriate for housekeeping 41,43 . As gene expression following TBI strongly depends on the trauma incidence and time point of healing, in line with previous studies 20, 36 , we strongly recommend an individual assessment of the reference genes for each study condition and time point analysed prior to target gene normalisation. Nonetheless, the application of different reference genes at continuative time points allows the analysis of time-dependent target gene expression, if the data is continuously normalised to the same control group. In studies primarily focussing on time-dependent changes of gene expression in primary organs following a specific trauma, the method of gene expression analysis needs to be adapted. Therefore, the evaluation of suitable reference genes is suggestively performed analysing the CT values of a specific trauma throughout all time points. While the reference gene evaluation in primarily traumatised organs appears to be very complex, in secondary organs, i.e. white adipose tissue (WAT), liver, muscle and spleen, Ppia expression remained continuously stable following isolated and combined trauma. The enzyme Ppia, also called cyclophilin A, catalyses the cis-trans isomerisation of proline imidic peptide bonds thereby accelerating protein folding 44 . In line with our results, the reference gene Ppia previously showed high stability in adipose tissue and liver even under challenging metabolic conditions 45,46 .
In WAT at d3 and d7 following trauma, our analysis suggests the additional use of the reference gene Psmb2 along with Ppia. The encoded protein of Psmb2 is part of the proteasome complex that controls protein degradation 47 . The application of Psmb2 has been previously approved for housekeeping in brown adipose tissue 48 . In liver at d3 and d7, our results suggest the additional use of B2m, which was already approved as reference gene by Gong et al. 45 . Since reference gene evaluation in secondary organs is indispensable but unfortunately uncommon 49-51 , we suggest the use of Ppia in combination with a second reference gene as a simple guideline. However, it has to be pointed out that our study is limited to specific types of trauma induced in young female C57Bl/6J mice. Despite the high degree of expression stability of Ppia, further investigations are recommended to verify Ppia as a suitable reference gene in secondary organs of other trauma mouse models.
In summary, our work raises awareness for the precise evaluation of the reference genes used in primary organs and points to Hprt as a stable housekeeper. Additionally, we provide a basic guideline for the gene expression analysis in secondary organs of a murine combined trauma model: Our study revealed Ppia as the most suitable reference gene in white adipose tissue, liver, muscle and spleen, which in combination with a second reference gene ensures the reliability of the normalised data. This study provides the basis for reliable and reproducible gene expression analysis for future biomedical research, to elucidate the pathophysiological mechanisms underlying post traumatic alterations with the aim to enhance specific therapy development in the field of regenerative trauma surgery.

Methods
Animals. Female C57Bl/6J, referred to as wildtype (WT) mice, were purchased from Charles River and maintained following the regulations of animal protection and care at the Research facilities for experimental medicine (FEM, Charité-Universitätsmedizin Berlin). Upon arrival, animals were kept under standard controlled conditions: temperature 20 ± 2 °C, 12 h light/dark cycle, water ad libitum and standard diet. All experimental procedures were approved by the local legal representative animal rights protection authorities (Landesamt für Gesundheit und Soziales Berlin, G0009/12) and performed according to the policies and principles established by the national institute of health guide for care and use of laboratory animals as well as the Welfare Act (Federal Law Gazett I, p.1094).  48 12-week old female mice (body weight: 20-25 g) were randomly assigned to the following groups (n = 6 for each group and time point): control, isolated fracture (Fx) and traumatic brain injury (TBI) as well as combined injury (TBI + Fx). Fracture was induced conducting a standardised mid-diaphyseal femoral osteotomy of the left limb. Therefore, an anterolateral skin incision was performed following an imaginal line from knee to hip joint for approximately 2 cm. Once the fascia lata was opened, the Musculus vastus lateralis and Musculus biceps femoris were bluntly dissected until the femoral bone was exposed. Now the external fixator (RISystem) was placed to drill the orthogonal pin holes by hand perpendicular to the longitudinal axis of the femur along the cortical surface. With the first pin hole just proximal of the distal metaphysis, the external fixator was screwed strictly parallel to the femur, allowing a rigid fixation. With the external fixator in place, a 0.70 mm osteotomy of was performed, applying a Gigli wire saw (RISystem) between the two middle pins. Throughout the whole procedure, the sciatic nerve was spared and its function postoperatively tested. TBI was induced applying a standardised controlled cortical impact model (CCI) to the parietotemporal cortex. Therefore, the neurocranium was positioned and fixed in a stereotactic device, allowing a skin incision sagittal and temporal on the left side of the skull. After the Musculus temporalis was mobilised and the cranial bone exposed, a 7 × 7 mm temporoparietal craniotomy was microsurgically performed, while ensuring an entirely intact dura mater. Now the actual TBI was caused, applying a stereotactic computer controlled pneumatically driven bolt (3 mm flat, 45-degree angle, impact velocity of 3.5 m/s, 0.25 mm penetration depth, contact duration of 0.15 s). After the impact, the preserved piece of skull was repositioned and adapted using dental cement to suppress cranial bone healing. In the combined injury group, TBI was induced first, which was followed by fracture. All wounds were closed with Ethilon 6.0 suture. Each trauma surgery consisted of perioperative anesthesia and antibiotic prophylaxis as well as postoperative analgesia and monitoring 52 .
Three (d3) and seven days (d7) after surgery directly (primary organs), as well as indirectly trauma-affected organs (secondary organs) were obtained, snap frozen in liquid nitrogen and stored at −80 °C. Primary organs include the intact femur/fracture callus and the hypothalamus, whereas secondary organs include white adipose tissue (WAT), liver, spleen and muscle. The fracture callus was gathered between the middle pins of the external fixator (groups: Fx, TBI + Fx), with an extraction of the intact femur accordingly (groups: control, TBI). The hypothalamus was dissected and sampled macroscopically after decapitation. For white adipose tissue a gonadal fat pad, for muscle the intact Musculus quadriceps femoris and for liver/spleen a 0.1 cm 3 part of the organ was obtained for RNA extraction.
RNA extraction and cDNA synthesis. Frozen tissues were homogenised in TRIzol utilising an UltraTurrax (Sigma Aldrich) and isolated with the RNeasy mini Kit (Qiagen) combined with DNase I (Thermo Fisher) treatment. RNA concentrations were determined and purity monitored (A260:A280 ratio of 1.9-2.1) using a NanoPhotometer P360 (Implen GmbH). A tissue-dependent volume of 0.5-5 µg RNA was reverse transcribed to complementary DNA (cDNA) applying the RevertAid First Strand cDNA Synthesis Kit (ThermoFisher) and finally stored at −20 °C. qRT-PCR. Quantitative real-time PCR (qRT-PCR) was performed within a 384 well-plate in a 7900HT Fast Real-Time PCR System employing the software SDS v2.4. (Applied Biosystems). Power SYBR Green PCR Master Mix (Sigma Aldrich) with 2 ng cDNA per well were run at 60° annealing temperature. For each run SDS v2.4. calculated the threshold cycle (Ct), as well as the melting curve, which confirmed PCR product specificity.
The murine qRT-PCR primers of the candidate reference (Table 1) and target genes ( Table 2) were individually designed using the GRCm38.p6 C57BL/6J reference genome and employing the software Primer3 (https :// bioin fo.ut.ee/prime r3-0.4.0/), spanning at least two exons with a large intron in between to avoid amplifications from genomic DNA. The primers were purchased from Eurofins Genomics GmbH. While the reference genes Actb, B2m, Gapdh, Hprt and Ppia were analysed in all tissue, Hmbs, Psmb2, Rplp0, Srfs4 and Tbp were applied tissue-specific. For the finally analysed target genes, the expression was normalised per individual reference genes by the delta-delta Ct method 53 .
Expression stability analysis. The candidate reference genes (Table 1) were examined for their expression stability, using the following Microsoft Excel-based algorithms: NormFinder 54 , GeNorm 28 and BestKeeper 55 . Once all individually calculated ranks were summed to identify the "final rank", the result was confirmed applying the web tool RefFinder (https ://www.heart cure.com.au/reffi nder/).

NormFinder analysis.
NormFinder incorporates the inter-and intra-group variances to determine the gene expression stability value (sv) 54 . The reference gene scoring the lowest sv is considered the most stable expressed. Finally, the algorithm ranks all assessed candidate reference genes based on their sv. GeNorm analysis. The GeNorm algorithm defines the most suitable reference gene pair whose expression ratio is least affected by external factors of varying samples and experimental conditions 28 . Therefore, GeNorm calculates the stability value (M), which is based on the pairwise variation of one gene compared to all other assessed candidate reference genes, with a stepwise exclusion of the highest M value. As a result, the candidates with the lowest M value are considered the reference genes with the highest gene expression stability. The GeNorm file suggests a threshold of 1.5 for suitable reference genes, however, when dealing with homogeneous sample panels the threshold can be set to M < 0.5 27 . Additionally, GeNorm determines the minimal number of reference genes necessary to allow a sufficient gene expression normalisation. The accuracy of normalisation will