Selection of suitable internal controls for gene expression normalization in rats with spinal cord injury

Wei Liu, Jie Yu, Yi-Fan Wang Qian-Qian Shan, Ya-Xian Wang

Abstract There is a lack of systematic research on the expression of internal control genes used for gene expression normalization in real-time reverse transcription polymerase chain reaction in spinal cord injury research. In this study, we used rat models of spinal cord hemisection to analyze the expression stability of 13 commonly applied reference genes: Actb, Ankrd27, CypA, Gapdh, Hprt1, Mrpl10, Pgk1, Rictor, Rn18s, Tbp, Ubc, Ubxn11, and Ywhaz. Our results show that the expression of Ankrd27, Ubc, and Tbp were stable after spinal cord injury, while Actb was the most unstable internal control gene. Ankrd27, Ubc, Tbp, and Actb were consequently used to investigate the effects of internal control genes with differing stabilities on the normalization of target gene expression. Target gene expression levels and changes over time were similar when Ankrd27, Ubc, and Tbp were used as internal controls but different when Actb was used as an internal control. We recommend that Ankrd27, Ubc, and Tbp are used as internal control genes for real-time reverse transcription polymerase chain reaction in spinal cord injury research. This study was approved by the Administration Committee of Experimental Animals, Jiangsu Province, China (approval No. 20180304-008) on March 4, 2018.

Key Words: geNorm analysis; reference genes; internal control genes; normalization; NormFinder analysis; reverse transcription-quantitative polymerase chain reaction; spinal cord injury; stability of gene expression

Introduction

Spinal cord injury impairs sensation and muscle functions and often leads to paraplegia or quadriplegia. Spinal cord injury has complex pathophysiology involving primary injury and subsequent secondary injury (Silva et al., 2014). Primary spinal cord injury refers to the initial mechanical injury and occurs within 2 hours, whereas secondary injury refers to a series of pathological processes, including apoptosis, hypoxia, edema, and inflammation, and can last weeks or even months (Anwar et al., 2016; Ahuja et al., 2017a). Secondary injury aggravates cell death and tissue damage and expands the zone and severity of tissue damage (Ahuja et al., 2017b; Zhang et al., 2020). These multifaceted injury events involve many types of cells, including neurons, glial cells (such as astrocytes, oligodendrocytes, and microglia), and inflammatory cells (such as macrophages, neutrophils, and T cells) (O’Shea et al., 2017; Huang et al., 2019). Recently, with the development of highthroughput analysis, gene profiles after spinal cord injury have been revealed, facilitating the understanding of the underlying molecular basis of and the discovery of potential targets for the treatment of spinal cord injury (Chen et al., 2013; Yu et al., 2019; Gong et al., 2020). However, for the precise determination of gene abundance, appropriate internal control genes should be used to normalize gene expression during quantification.

In quantitative real-time reverse transcription polymerase chain reaction (RT-PCR) experiments, internal controls are used to account for variations in gene expression determination caused by variations in experimental procedures and RNA sample quantity and quality (Vandesompele et al., 2002). Many reference genes have been identified as suitable internal controls in diverse species, tissues, and physiological and pathological processes (Benak et al., 2019; Herath et al., 2020; Liu et al., 2020; Xia et al., 2021). Nevertheless, there are few consensuses on proper internal controls for studies of spinal cord injury and repair. Therefore, we investigated the expression stability of 13 reference genes in normal or injured spinal cord with the aim to select optimal internal control genes following spinal cord injury. The 13 reference genes wereActb(beta-actin),Ankrd27(ankyrin repeat domain 27),CypA(cyclophilin A),Gapdh(glyceraldehyde-3-phosphate dehydrogenase),Hprt1(hypoxanthine phosphoribosyltransferase 1),Mrpl10(mitochondrial ribosomal protein L10),Pgk1(phosphoglycerate kinase 1),Rictor(RPTOR independent companion of MTOR, complex 2),Rn18s(18S ribosomal RNA),Tbp(TATA box binding protein),Ubc(ubiquitin C),Ubxn11(UBX domain protein 11), andYwhaz(tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta).

Materials and Methods

Animals

The animal experiments were ethically approved by the Administration Committee of Experimental Animals, Jiangsu Province, China (approval No. 20180304-008) on March 4, 2018. Surgical procedures were performed following the Institutional Animal Care Guidelines of Nantong University, Nantong, China. All experiments were designed and reported according to the Animal Research: Reporting ofIn VivoExperiments (ARRIVE) guidelines.

A total of 45 adult female Sprague-Dawley rats (specificpathogen-free; 2 months old; 180-220 g) were purchased from the Animal Center of Nantong University (license no. SCXK [Su] 2019-0001 and SYXK [Su] 2017-0046). The rats were randomized into experiment and sham groups. Rats were anesthetized by intraperitoneal injection of mixed narcotics (MilliporeSigma, Burlington, MA, USA; Merck KGaA, Darmstadt, Germany) containing 85 mg/kg trichloroacetaldehyde monohydrate, 42 mg/kg magnesium sulfate, and 17 mg/kg sodium pentobarbital (Wang et al., 2018, 2019). A longitudinal incision was made in the dorsal midline skin, and the spinal cord was exposed by a laminectomy at vertebral segments T9-10. For the experiment group, a hemisection was performed at vertebral segment T9 with an ophthalmic iris knife (Beaver-Visitec International, Waltham, MA, USA). For the sham group, the hemisection was not performed. All rats were housed in temperature- and humidity-controlled home cages with free access to water and food after surgery.

Immunostaining observation

Rats were anesthetized with mixed narcotics at 1, 4, 7, 14, and 21 days after spinal cord injury and were transcardially perfused with 4% paraformaldehyde. Spinal cord segments at the caudal region of the injured site were harvested, dehydrated in 30% sucrose, embedded in Tissue-Tek O.C.T. compound (Sakura Finetek USA, Inc., Torrance, CA, USA), and coronally cryo-sectioned with a cryostat microtome (Leica, Wetzlar, Germany). The 20-μm spinal cord slices were mounted on microscope slides and incubated in primary antibody overnight at 4°C. The primary antibodies were rabbit anti-glial fibrillary acidic protein (GFAP; 1:400, Cat# Z0334, RRID:AB_10013382; Agilent Technologies, Santa Clara, CA, USA) or rabbit anti-ionized calcium binding adapter molecule 1 (IBA1; 1:200, Cat# 019-19741, RRID:AB_839504; FUJIFILM Wako Shibayagi Corporation, Osaka, Japan). This was followed by incubation with the secondary antibody Alexa Fluor 488-conjugated anti-rabbit IgG (1:400, Cat# ab150077; Abcam, Cambridge, UK) for 2 hours at room temperature. Spinal cord slices were photographed under fluorescence microscopy (Axio Imager M2, Carl Zeiss Meditec AG, Jena, Germany).

Quantitative real-time RT-PCR

At 1, 4, 7, and 14 days after surgery, 5-mm-long rat spinal cord segments at the caudal region were collected after laminectomy or hemisection. Total RNA samples were prepared from spinal cord segments with TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was confirmed by agarose gel electrophoresis. RNA samples (1 μg) were reverse transcribed to complementary DNA using HiScript III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd., Nanjing, China). Sequence-specific primers were designed according to previous publications (Peinnequin et al., 2004; Bonefeld et al., 2008; Martínez-Beamonte et al., 2011; Bangaru et al., 2012; Gambarotta et al., 2014) or with Primer Express software (v3.0.1; Thermo Fisher Scientific) and synthesized by Sangon Biotech (Shanghai) Co., Ltd. (Shanghai, China). The primer sequences of the candidate reference genes tested are inTable 1. Primer sequences ofNONRATT020851.2were (sense) 5′-TTG GCA ATG GCA GTG TAT GTC-3′ and (antisense) 5′-TCA AGA GCA ACG GAA CAA GAA C-3′, and primer sequences ofNONRATT007410.2were (sense) 5′-GCA CGG TCA GAG GTC AGT T-3′ and (antisense) 5′-GAT GGC GGC TCA GTG GTA A-3′. Quantitative RT-PCR was performed with ChamQ Universal SYBR qPCR Master Mix (Vazyme) and specific primers in a StepOne Real-Time PCR System (Thermo Fisher Scientific). The thermocycler program was one initial denaturation cycle at 95°C for 5 minutes followed by 40 cycles of 95°C for 15 seconds and 60°C for 30 seconds. Each sample was tested in triplicate and analyzed using the 2-∆∆CTmethod (Livak and Schmittgen, 2001) to determine the relative gene expression.

Table 1 ｜Primer sequences of candidate reference genes

Bioinformatic analysis

Raw Ct values of candidate internal control genes were analyzed using geNorm and NormFinder bioinformatic software. geNorm analysis (v3.4; https://genorm.cmgg.be/) was applied using the Microsoft Excel 2007 Visual Basic plug-in (Microsoft Corporation, Redmond, WA, USA) to calculate the gene expression stability value (M-value), which is negatively correlated with gene expression stability. Additionally, we used geNorm software to calculate the pairwise variation (Vn/n+1), a normalization factor. Vn/n+1less than 0.15 indicates that no further (n + 1) internal control genes are needed (Vandesompele et al., 2002). NormFinder analysis (v19; https://moma.dk/normfinder-software) was applied using the Microsoft Excel 2007 Visual Basic plug-in (Microsoft Corporation) to calculate intra-group and intergroup variations and estimate expression stability values (Andersen et al., 2004).

Statistical analysis

Statistical analyses were performed using GraphPad Prism 6.0 (GraphPad Software, San Diego, CA, USA). Differences between the experiment and control groups were analyzed by one-way analysis of variance followed by the Holm-Šídák multiple comparisons test. Significant differences were set atP＜ 0.05.

Results

Spinal cord hemisection elicits significant cellular changes

Rat spinal cord segments at the caudal region of the injured site were collected at 1, 4, 7, 14, and 21 days after spinal cord hemisection at vertebral segment T9 and immunostained to observe injury-induced morphological changes. The expressions of two main glial cell types, astrocytes and microglia, were visualized by immunofluorescence assay, as glial cells make up a large proportion of cells in mammals and play essential roles after spinal cord injury (Bahney and von Bartheld, 2018).

Few astrocytes were present at 1 day after spinal cord injury by immunostaining with the astrocyte marker GFAP (Zhang et al., 2017). GFAP immunopositivity was elevated at later time points, reaching its peak at 7 days after injury and gradually reduced toward the basal value at 14 and 21 days after spinal cord injury (Figure 1A). Immunostaining from the astrocyte marker GFAP was mainly localized in the gray matter, whereas immunostaining from the microglia marker IBA1 (Waller et al., 2019) was observed in the gray and white matter. The immunopositivity of IBA1 was elevated at 4 days after spinal cord injury and recovered to the basal level at 21 days after injury (Figure 1B), similar to GFAP results. These observations showed that there are robust morphological and molecular changes after spinal cord hemisection. These remarkable cellular changes indicate that many genes, including reference genes, might be differentially expressed after spinal cord injury.

Changes in candidate reference gene expression in spinal cord

To select potential internal controls for the normalization of target genes, the commonly used reference genesActb,Ankrd27,CypA,Gapdh,Hprt1,Mrpl10,Pgk1,Rictor,Rn18s,Tbp,Ubc,Ubxn11, andYwhazwere analyzed by RT-PCR. All these reference genes had linear correlation coefficients (R2) close to 1 and amplification efficiencies between 90% to 110%, which suggests high efficiency and quality of the synthesized primers of these reference genes (Table 2).

Table 2 ｜The amplification efficiencies of candidate reference genes

Next, the raw Ct values of these reference genes were determined (Figure 2). The raw Ct values ofRn18swere less than 7, which are the lowest values among the tested reference genes and are consistent with the high abundance of theRn18sribosomal RNA gene.Ubxn11had the highest raw Ct values (～28), and other reference genes had raw Ct values ranging from 16 to 25. Notably, some reference genes, such asRn18sandActb, showed curved trends in their raw Ct values, especially in samples from spinal cord after hemisection. Conversely, the raw Ct values of some reference genes, such asAnkrd27,Gapdh,Mrpl10,Tbp,Ubc, andYwhaz, did not exhibit obvious changes.

Determination of expression stability of candidate reference genes

The raw Ct values of the 13 reference genes were analyzed using geNorm to determine gene expression stability, and the M-value, which negatively correlates with gene stability, was calculated.Actbhad the highest M-value, which suggests that the stability of the expression ofActbwas relatively low. Conversely,Ankrd27andUbchad the lowest M-values and the highest gene expression stabilities (Figure 3A). Calculation of pairwise variation (Vn/n+1), a normalization factor that advises the ideal number of internal control genes, showed that the pairwise variation of V2/3was low (Figure 3B). Bioinformatic analysis using NormFinder indicated thatTbp,Ankrd27, andUbcwere the most stable genes, whileRn18sandActbwere the least stable genes in spinal cord segments following rat spinal cord injury (Figure 4). Therefore, the reference genesTbp,Ankrd27, andUbccan be used as internal controls for the determination of gene expressions in injured spinal cord segments.

Application of stable or unstable candidate reference genes as internal control genes

The candidate reference genes were used to determine expression patterns of the genes NONRATT020851.2 andNONRATT007410.2, which were randomly selected from previously obtained sequencing data of rat spinal cord segments at 1, 4, and 7 days after injury (stored in the National Omics Data Encyclopedia [NODE], https://www.biosino.org/node/, accession number OEP000369) (Yu et al., 2019). Sequencing data showed that the expression levels ofNONRATT020851.2were elevated at 1, 4, and 7 days after spinal cord injury compared with the control (P＜ 0.05), reaching a peak at 4 days after injury (Figure 5A). The expression patterns ofNONRATT020851.2were different when usingActb, a reference gene with low stability, as the internal control for normalization. The relative abundance ofNONRATT020851.2decreased at 1 day after injury and increased to peak levels at 7 days after injury when using the unstable reference gene Actb to calculate ∆∆Ct (Figure 5B). Conversely, the use of the stable reference genesTbporAnkrd27+Ubcas internal controls achieved similar trends ofNONRATT020851.2expression as sequencing data (Figure 5CandD). Similarly, the expression patterns of NONRATT007410.2 were more comparable to sequencing observations when usingTbporAnkrd27+Ubcas internal controls instead ofActb(Figure 5E-H).

Figure 1｜ Immunostaining of rat spinal cord segments at the caudal region at 1, 4, 7, 14, and 21 days following spinal cord injury.

Figure 2｜ Raw Ct values of candidate reference genes at 1, 4, 7, and 14 days following spinal cord injury.

Figure 3｜geNorm-determined expression stability of candidate reference genes.

Figure 4｜NormFinder-determined gene expression stability of candidate reference genes.

Figure 5｜Relative expression of target genes NONRATT020851.2 and NONRATT007410.2 at 1, 4, and 7 days following spinal cord injury.

Discussion

qPCR is the gold standard for the determination of gene abundance and the validation of large-scale transcript data, such as from microarrays or sequencing (Derveaux et al., 2010). qPCR has advantages, such as high sensitivity and simplicity, and disadvantages, including stringent sample quality requirements, operator variability, and the inefficiency of reverse transcription (Bustin and Nolan, 2004). The use of internal controls largely reduces variations caused by technological defects or operational errors and increases the reliability and accuracy of RT-PCR results (Kadegowda et al., 2009).

Reference genes are essential genes with relatively constant abundance in all cells for various biological activities (Eisenberg and Levanon, 2013). Therefore, reference genes are commonly used as internal controls for calibration of RTPCR. However, emerging studies have demonstrated that reference gene expression is tissue- and condition-dependent. For instance,ActbandGapdh, two widely used reference genes, had dissimilar histological expression profiles in the liver, the vascular system, and the digestive tract (Lin and Redies, 2012). Our previous sequencing analyses also revealed that the expression level ofActbwas upregulated in rat sciatic nerve segments after peripheral nerve injury, which suggests that the expression of reference genes may be affected by pathophysiological changes (Wang et al., 2018). Recently, Košuth et al. (2019) analyzed the expression stability of several reference genes in the rat spinal cord using two different injury models, minimal and contusion spinal cord injury. They foundGapdhandHprt1were suitable for contusion spinal cord injury experiments, whileGapdhwas suitable for minimal spinal cord injury experiments. The expression stability of reference genes depended on the severity of the injury. Therefore, it is of great importance to select reference genes as internal controls for normalization that are stably expressed in the target cells, tissues, and biological conditions of interest.

Previously, we examined the stability of many frequently used reference genes in different tissues after rat peripheral nerve injury. We discovered thatMrpl10,Tbp, andHprt1are suitable internal controls for sciatic nerve segments, dorsal root ganglion samples, and gastrocnemius muscles, respectively (Wang et al., 2017, 2019). In the current study, we focused on spinal cord injury—injury to the central nervous system that induces even more serious consequences—and measured the expression of reference genes in spinal cord segments after central nervous system injury. For all tested reference genes, geNorm-calculated gene expression stability M-values were less than the default limit of 1.5 and pairwise variation V-values were less than the default limit of 0.15 (Bangaru et al., 2012). Therefore, the tested reference genes met the basic requirements for use as internal controls. Raw Ct values and geNorm and NormFinder outcomes showed thatTbp,Ankrd27, andUbcwere the most stably expressed reference genes in spinal cords injured by hemisection. Moreover, a comparison study of the gene expression ofNONRATT020851.2andNONRATT007410.2using the stable reference genesTbp,Ankrd27, andUbcand an unstable reference geneActbfurther validated and supported the expression stability ofTbp,Ankrd27, andUbc. The outcome of our current study is suitable for the spinal cord hemisection model, but we still need to find suitable internal control genes for use in other injury models, such as contusion or full transection. We used female rats in our study, as their use is preferred in spinal cord injury research (Onifer et al., 2007; Stewart et al., 2020). Many sex-related genes contribute to morphological, physiological, and behavioral differences between males and females (Shi et al., 2019), so it is worth exploring in the future whether sex differences affect the selection of internal control genes.

In sum, the identification of suitable internal controls provides guidance for RT-PCR experiments and facilitates the exposure and understanding of the underlying genetic infrastructure of spinal cord injury and regeneration.

Author contributions:Study conception and design and manuscript draft: YXW; experiment implementation and data analyses: WL, JY, YFW, QQS; reagents/materials/analysis tools support: QQS, YXW. All authors read and approved the final manuscript.

Conflicts of interest:The authors declare that they have no conflict of interest.

Financial support:This study was supported by the National Natural Science Foundation of China, No. 81901257 (to YXW); the Natural Science Foundation of Jiangsu Province of China, No. BK20180951 (to YXW); Postgraduate Research and Practice Innovation Program of Jiangsu Province of China, No. KYCX20_2818 (to WL); and Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD, to Key Laboratory of Neuroregeneration of Jiangsu and Ministry of Education). The funding sources had no role in study conception and design, data analysis or interpretation, paper writing or deciding to submit this paper for publication.

Institutional review board statement:This study was approved by the Administration Committee of Experimental Animals, Jiangsu Province, China (approval No. 20180304-008) on March 4, 2018.

Copyright license agreement:The Copyright License Agreement has been signed by all authors before publication.

Data sharing statement:Datasets analyzed during the current study are available from the corresponding author on reasonable request.

Plagiarism check:Checked twice by iThenticate.

Peer review:Externally peer reviewed.

Open access statement:This is an open access journal, and articles are distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 License, which allows others to remix, tweak, and build upon the work non-commercially, as long as appropriate credit is given and the new creations are licensed under the identical terms.