Analysis of transcriptome sequencing of sciatic nerves in Sprague-Dawley rats of different ages.

An aging-induced decrease in Schwann cell viability can affect regeneration following peripheral nerve injury in mammals. It is therefore necessary to investigate possible age-related changes in gene expression that may affect the biological function of peripheral nerves. Ten 1-week-old and ten 12-month-old healthy male Sprague-Dawley rats were divided into young (1 week old) and adult (12 months old) groups according to their ages. mRNA expression in the sciatic nerve was compared between young and adult rats using next-generation sequencing (NGS) and bioinformatics (n = 4/group). The 18 groups of differentially expressed mRNA (DEmRNAs) were also tested by quantitative reverse transcription polymerase chain reaction (n = 6/group). Results revealed that (1) compared with young rats, adult rats had 3608 groups of DEmRNAs. Of these, 2684 were groups of upregulated genes, and 924 were groups of downregulated genes. Their functions mainly involved cell viability, proliferation, differentiation, regeneration, and myelination. (2) The gene with the most obvious increase of all DEmRNAs in adult rats was Thrsp (log2FC = 9.01, P < 0.05), and the gene with the most obvious reduction was Col2a1 (log2FC = -8.89, P < 0.05). (3) Gene Ontology analysis showed that DEmRNAs were mainly concentrated in oligosaccharide binding, nucleotide-binding oligomerization domain containing one signaling pathway, and peptide-transporting ATPase activity. (4) Analysis using the Kyoto Encyclopedia of Genes and Genomes showed that, with increased age, DEmRNAs were mainly enriched in steroid biosynthesis, Staphylococcus aureus infection, and graft-versus-host disease. (5) Spearman's correlation coefficient method for evaluating NGS accuracy showed that the NGS results and quantitative reverse transcription polymerase chain reaction results were positively correlated (rs = 0.74, P < 0.05). These findings confirm a difference in sciatic nerve gene expression between adult and young rats, suggesting that, in peripheral nerves, cells and the microenvironment change with age, thus influencing the function and repair of peripheral nerves.

Introduction

The treatment of limb paralysis and dysfunction caused by peripheral nerve injury is a major challenge in medical science. In the USA alone, the number of peripheral nerve injury surgeries performed each year has reached 50,000 cases, and the economic burden is over 7 billion US dollars. Therefore, the question of how to improve treatment and effectively repair peripheral nerve injury has become a major public health issue that needs to be addressed by the global healthcare community (Wang et al., 2010; Yang et al., 2011; Chen et al., 2018). Thus, studies of peripheral nerve regeneration have significance in both scientific research and clinical practice (Zheng et al., 2014; He et al., 2015; Hung et al., 2015; Qiu et al., 2015; Pan et al., 2017; Zou et al., 2018).

The reduction in nerve repair ability in older animals is a research hotspot in the field of peripheral nerve injury repair (Verdu et al., 2000; Amer et al., 2014). For a long time, the majority of researchers have considered that a decline in neuronal regeneration ability is the major reason that axonal regeneration is slow in old animals (Graciarena et al., 2014; Moldovan et al., 2016; Lim et al., 2017). For example, Zou et al. (2013) showed that intrinsic changes in older neurons influence neuronal regeneration ability in a Caenorhabditis elegans model. They found that in older neurons, let-7 down-regulating LIN-41, which contributes to a developmental decline in axonal regeneration. However, in the younger neurons, LIN-41 inhibits let-7 expression via Argonatute ALG-1 to ensure that axonal regeneration could only be inhibited in older neurons. Some other scholars focused on cells other than neurons. For example, Painter et al. (2014) demonstrated that the peripheral nerve regeneration ability of 24-month-old mice is worse than that of 2-month-old mice; although in vitro gene detection revealed that only 14 genes in the dorsal root ganglion cells of old mice are significantly different from those in young mice. In addition, genes that are directly associated with regeneration ability (such as ATF3, GAP and prr1a) did not markedly differ between old and young mice. Therefore, some researchers have suggested that the major factor influencing peripheral nerve regeneration in mammals such as mice, rats, and humans is not the neurons themselves, but instead the decline of Schwann cell viability in peripheral nerves (Michio et al., 2014; Couve et al., 2017; Xu et al., 2017).

This study therefore aimed to compare mRNA expression in the sciatic nerves of Sprague-Dawley rats at different ages using next-generation, high-throughput whole RNA sequencing and bioinformatics. Using these techniques, the aim was to discover age-related molecules that influence the biological functions of peripheral nerves, and to investigate any underlying mechanisms, to provide theoretical bases for improving the treatment of aging peripheral nerves after injury (Canta et al., 2016; Sakita et al., 2016; Zhou. et al., 2017).

Materials and Methods

Animals

Ten 1-week-old and ten 12-month-old healthy male Sprague-Dawley rats, bought from the Animal Center of Medical School of Sun Yat-Sen University of China, were randomly selected as experimental animals. This study was approved by the Experimental Animal Administration Committee of Sun Yat-Sen University (approval number: (2013)A-055). Efforts were taken to minimize animal suffering during the experiment.

Collection of sciatic nerves

Four rats in the young or adult group were randomly selected for the next-generation RNA high-throughput sequencing detection of sciatic nerves. Rats in the adult group were deeply anesthetized using an intraperitoneal injection of 10% chloral hydrate (Sinopharm Chemical Reagent Co., Ltd., Shanghai, China; 0.3 mL/100 g body weight). Rats in the young group were sacrificed by cervical dislocation. The specific operational procedures have been previously published (Zhu et al., 2015, 2017). Under aseptic conditions, the skin of the left leg was cut parallel to the femur, and the sciatic nerve was exposed by splitting the superficial gluteus muscle. The bilateral sciatic nerves (≥ 15 mm in length) were excised from the lower edge of piriform muscle from the rats and then external fat and connective tissue were removed.

Library preparation and RNA sequencing

Comparison of the senescence pathway genes was performed by next-generation sequencing (NGS). All sequencing programs were performed by Shanghai Biotechnology Corporation (Shanghai, China). Briefly, RNA extraction was performed using the mir Vana™ miRNA isolation kit (Cat# AM1561, Ambion, Austin, TX, USA) following the standard operating procedures provided by the manufacturer. The quality of the extracted RNA was detected using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Total RNA was then purified using the RNAClean XP kit (Cat A63987, Beckman Coulter, Inc., CA, USA) and the RNase-Free DNase Set (Cat#79254, Qiagen, Germany). The quality of the initial total RNA sample for the sequencing experiment was detected using a NanoDrop ND-2000 spectrophotometer and an Agilent Bioanalyzer 2100. Total RNA that passed the quality control was used in subsequent sequencing experiments (Pertea et al., 2015).

For transcriptome sequencing, we used the number of reads mapping to the gene region to estimate the gene expression level (Pertea et al., 2016). However, in addition to being proportional to the gene expression level, the number of reads was also associated with the length of genes and the sequencing data volume. To compare gene expression levels between different genes and different samples, reads were converted into fragments per kilobase of exon model per million mapped reads (FPKM) for normalization of gene expression levels (Mortazavi et al., 2008; Yu et al., 2017). First, the number of fragments of each gene after Hisat2 comparison was counted using Stringtie v1.3.0 (The Center for Computational Biology at Johns Hopkins University, Baltimore, MD, USA). Next, normalization was performed using the trimmed mean of M values method (Robinson et al., 2010). Finally, the FPKM value of each gene was calculated using the Perl program.

Quantitative real time polymerase chain reaction (RT-PCR)

Quantitative analysis of targeted mRNA expression was performed using RT-PCR and the relative standard curve method. Six rats in each group were randomly selected. The bilateral sciatic nerves were collected using the previously described method, and total RNA was extracted. The RNA concentration and purity were detected using a Nanodrop 2000. RNA samples with high concentrations were diluted to a final concentration of 200 ng/μL. RNA (1 μg) was used for RNA reverse transcription using the RevertAid first strand cDNA synthesis kit (Thermo Fisher, Waltham, MA, USA). An appropriate amount of cDNA was amplified using FastStart Universal SYBR Green Master Mix (Roche, Basel, Switzerland) in a fluorescence quantitative PCR machine (Stepone Plus, Thermo Fisher). The specific procedures were performed according to the manufacturer's product manual. Each sample was run in triplicate. GAPDH expression was used to normalize mRNA expression. Gene information and primers are shown in . The specificity of RT-PCR was confirmed via routine agarose gel electrophoresis and melting curve analysis. The 2-ΔΔCt method was used to calculate relative gene expression (Huang et al., 2017).

Bioinformatics and statistical analyses

Analyses of differential genes among samples were performed using edgeR (Robinson et al., 2010). After the P-value was obtained, multiple testing was performed for correction. The threshold value of the P-value was determined by controlling for false discovery rate (Benjamini et al., 1995, 2001). The corrected P-value was the Q-value. In addition, the differential expression fold was calculated based on the FPKM value, which was the fold change. The screening condition of differential genes was limited to Q-value ≤ 0.05 and fold change ≥ 2. Both gene ontology (GO) analysis (http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (http://www.genome.jp/kegg/) were performed to determine the roles of the differentially expressed mRNAs. Fisher's exact test was used to test the significance of GO terms and pathway identifiers enriched in the differentially expressed gene list. RNA sequencing and bioinformatics analyses were performed by Shanghai Biotechnology Corporation (Shanghai, China).

The gene data of RT-PCR conform to normal distribution and were presented as the mean ± SD. SPSS 13.0 software (SPSS, Chicago, IL, USA) was used to analyze the accuracy of NGS. Quantitative RT-PCR results were compared with sequencing results using Spearman's correlation coefficient method. P < 0.05 was considered statistically significant.

Results

Analyses of mRNA expression levels in the sciatic nerve of young and adult rats

Transcriptome gene expression of sciatic nerves in young and adult rats was analyzed using NGS. A scatter map and a heat map were used to present the distribution of differentially expressed mRNAs (DEmRNAs; ). A total of 30,957 mRNA groups were detected using RNA sequencing (RNA-seq), and the DEmRNA screening was performed using fold change ≥ 2.0 and P < 0.05 as the criteria. Compared with young rats, adult rats had 3608 groups of DEmRNAs. Of these, 2684 groups were upregulated genes, and 924 groups were downregulated genes. The gene with the greatest mRNA increase in adult rats was Thrsp (log2FC = 9.01, P < 0.05), and the gene with the greatest reduction was Col2a1 (log2FC = –8.89, P < 0.05). The 10 groups of genes with the most obvious DEmRNAs are shown in . Additionally, 208 mRNAs were only expressed in adult rats and 50 mRNAs were only expressed in young rats. The gene with the greatest mRNA increase in adult rats was Rn60_X_0752.3 (FPKM = 227.23, P < 0.05), and the gene with the most obvious increase in young rats was Ddx3y (FPKM = 32.69, P < 0.05). The top 10 most highly expressed mRNAs which were not expressed by rats in the opposite group are shown in .

Scatter map and heat map of mRNA expression in the sciatic nerve in young (1 week old, n = 4) and adult (12 months old, n = 4) Sprague-Dawley rats.

(A) Scatter map of mRNA expression. Red indicates upregulated genes in adult rats compared with young rats, while blue indicates downregulated genes. (B) Heat map of mRNA expression. Green in the heat map indicates a fold change ≤ 2, and red indicates a fold change > 2.

GO and KEGG analyses of DEmRNAs

Clustering analysis of DEmRNAs can help to identify the function of unknown transcripts, or the unknown function of known transcripts, by gathering similar expression patterns or similar genes to a class (Yao et al., 2012). Functional classification by GO is an internationally standardized classification system for gene function that offers a dynamic, updated, and controlled vocabulary, employing strictly defined concepts to comprehensively describe the properties of genes and their products in any organism. GO classification encompasses three domains, as follows: molecular function, cellular component, and biological process (Yi et al., 2006). GO analysis revealed that DEmRNAs were mainly concentrated in genes related to oligosaccharide-binding, nucleotide-binding oligomerization domain containing one (NOD1) signaling pathway, peptide-transporting ATPase activity, fibrinogen binding, and marginal zone B cell differentiation ( and ).

Schematic diagram of the Gene Ontology (GO) category and gene enrichment of differentially expressed mRNAs (DEmRNAs) in the sciatic nerve of 12-month-old and 1-week-old Sprague-Dawley rats.

(A) GO classification: Red columns represent biological process; green columns represent cellular component; blue columns represent molecular function. (B) GO enrichment: Different shapes indicate different types of function. Circles represent biological process; triangles represent cellular component, and squares represent molecular function. The area of the circles represents the number of DEmRNAs involved in the terms. Different colors indicate different Q-values of DEmRNAs in the GO terms between adult rats and young rats. Green represents Q > 0.10; yellow represents Q > 0.05 and Q < 0.075; orange represents Q < 0.05 and Q > 0.025; and red represents Q < 0.025. The standard of the test was set as 0.05 (n = 4 rats per group).

The pathway enrichment of DEmRNAs was analyzed using KEGG. The KEGG pathway database comprises information about networks of molecular interactions for numerous organisms, thus permitting functional classification. In the adult rats, DEmRNAs were mainly enriched in steroid biosynthesis, Staphylococcus aureus infection, graft-versus-host disease, type 1 diabetes mellitus, and allograft rejection categories ( and ).

Schematic diagram of the pathway classification and gene enrichment of differentially expressed mRNAs (DEmRNAs) in the sciatic nerve of 12-month-old and 1-week-old Sprague Dawley rats.

(A) Kyoto Encyclopedia of Genes and Genomes (KEGG) classification: Orange columns represent cellular process; brown columns represent environmental information processing; green columns represent genetic information processing; blue columns represent metabolism; and purple columns represent organismal systems. (B) KEGG pathway enrichment: The area of the circles represents the number of DEmRNAs involved in the pathway. Different colors indicate different Q-values of DEmRNAs in the pathway between adult rats and young rats. Green represents Q > 0.15; yellow represents Q > 0.10 and Q < 0.15; orange represents Q < 0.10 and Q > 0.05; and red represents Q < 0.05. The standard of the test was set as 0.05 (n = 4 per group).

Validation of microarray sequencing accuracy using quantitative RT-PCR

Five groups of upregulated mRNAs, namely nerve growth factor receptor (NGFR), betacellulin (BTC), activin A receptor type 1C (ACVR1C), disks large homolog 2 (DLG2), and Jun proto-oncogene (JUN), and five groups of downregulated mRNAs, namely noncompact myelin associated protein (NCMAP), myelin basic protein (MBP), PDZ binding kinase (PBK), baculoviral IAP repeat containing 5 (BIRC5), and early growth response 2 (EGR2), were selected for quantitative RT-PCR analyses in rats with different genotypes. The housekeeping gene GAPDH was used as the internal control to compare the sequencing results. Among the 10 groups of DEmRNAs, only BTC expression differed from the sequencing results. The sequencing results showed that BTC expression was upregulated, whereas the quantitative RT-PCR results demonstrated that BTC expression was downregulated. In contrast, the fold changes of expression in the quantitative RT-PCR results of other DEmRNAs were similar to those of the RNA microarray sequencing results ().

Comparison of quantitative RT-PCR (n = 6/group) and RNA sequencing (RNA-seq; n = 4/group) results of differentially expressed mRNAs in the sciatic nerve of 12-month-old and 1-week-old Sprague-Dawley rats.

To achieve baseline consistency, the quantitative RT-PCR fold change data were expressed as 2-ΔΔCt, and RNA-seq data are expressed as log2FC. NGFR: Nerve growth factor receptor; ACVR1C: activin A receptor type 1C; DLG2: discs large MAGUK scaffold protein 2; JUN: Jun proto-oncogene; BTC: betacellulin; NCMAP: noncompact myelin associated protein; MBP: myelin basic protein; BIRC5: baculoviral IAP repeat containing 5; EGR2: early growth response 2; PBK: PDZ binding kinase; FC: fold change; RT-PCR: real time polymerase chain reaction.

A total of 18 groups of DEmRNAs were randomly selected, including the above 10 groups of genes, as well as KIF2C, MAP3K, PDGFB, MPZ, PMP22, MAL, WWC1 and CDC20. To evaluate the accuracy of NGS sequencing, the correlation between RNA sequencing results and quantitative RT-PCR results was evaluated using Spearman correlation coefficient analyses. The correlation coefficient rs was 0.74 and P < 0.05, indicating that the NGS results and quantitative RT-PCR results were positively correlated ().

Scatter map of Spearman correlation coefficient analyses between quantitative real time polymerase chain reaction and sequencing data (Seq) of differentially expressed mRNAs (DEmRNAs) in the sciatic nerve of 12-month-old (adult) and 1-week-old (young) Sprague-Dawley rats.

The correlation coefficient rs was 0.74, and P < 0.05 (n = 6/group).

Discussion

After peripheral nerve injury, neurons initiate axonal regeneration. Distal axons, cut from the neuronal cell body, develop Wallerian degeneration (Szepanowski and Kieseier, 2016; Xiao et al., 2016; Roberts et al., 2017; Cervellini et al., 2018). Proximal axons form growth cones and extend to distal ends, which grow into the target organ to achieve functional recovery (Scheib et al., 2013; Kızılay et al., 2016; Ma et al., 2016; Wang et al., 2016; Tallon and Farah, 2017). However, with increasing age, the speed of nerve regeneration and functional recovery greatly decreases (Kawabuchi et al., 2011). We previously reviewed literature on the repair efficacy of peripheral nerve injury in the upper limb from a 20-year period, to perform a univariate meta-analysis. The results indicated that a 1-year increase in age had an odds ratio of 0.97–0.98 for excellent recovery (He et al., 2014). For example, if the regeneration ability following nerve injury in a newborn is 100%, the nerve regeneration ability after injury during old age (e.g., in a 60-year-old male) is only 16–30% of that at birth. These results suggest that the “aging” factor has a relatively large influence on nerve injury repair efficacy (Esquisatto et al., 2014; Ugrenović et al., 2016; Scheib et al., 2016; Castelnovo et al., 2017).

Some researchers suggested that the major factor influencing peripheral nerve regeneration in mammals, such as mice, rats, and humans, is not neurons themselves but instead the decline of Schwann cell viability in peripheral nerves (Belin et al., 2017; Plaza-Zabala et al., 2017). Kang et al. (2013) used phase-contrast microscopy and immunohistochemistry to detect the rate of nerve regeneration in mice at different ages. The results showed that the axonal regeneration rate of old animals does not markedly differ from that of young animals. However, endoneurial tubes in old animals have more nerve and myelin sheath debris to block the growth of axons. Therefore, these researchers considered that the “dedifferentiation-proliferation-differentiation” timing of Schwann cells, and degree of Schwann cells impact the time course of Wallerian degeneration to further influence the rate of axonal regeneration. Painter et al. (2014) compared the differences in gene expression in Schwann cells between old and young mice. They demonstrated that the expression levels of growth-factor-related coding genes (BTC, NGFR, and BDNF) and mitosis-related genes (KIF2C, PBK, BIRC5, and CDC20) markedly decrease in old mice, while expression levels of myelination-associated genes (PMP2, MPZ, MAL, and EGR2) increase in old mice. These results suggest that the dedifferentiation ability of Schwann cells declines in old mice, thereby reducing the rate of nerve regeneration.

In the current study, we compared transcriptome gene expression in the sciatic nerve of rats at different ages and found that 3608 groups of mRNAs were differentially expressed between adult rats (12 months old) and young rats (1 week old). The upregulated genes were mainly concentrated in changes in development, anatomical morphology, cell differentiation, and the extracellular matrix. To validate the accuracy of microarray sequencing, we selected 10 groups of genes associated with proliferation, differentiation, myelination, injury, and synaptic transmission for quantitative RT-PCR. The results were compared with sequencing results to validate their accuracy. The proliferation- and differentiation-related genes included NGFR, ACVR1C, BTC, survivin (BIRC5) and EGR2. The injury-related genes included JUN, and the myelination-related genes included MBP and NCMAP. The synaptic transmission-related genes included DLG2. The PCR results demonstrated that the expression trends of all groups of genes, except for BTC, were consistent with the microarray results. Spearman correlation coefficient analysis showed that the results of these two groups showed a positive correlation, which further indicated the accuracy of the microarray sequencing. The growth, development, conduction, and repair of peripheral nerves is a complex process. Although more than 1000 groups of differential mRNA sequences were obtained from the NGS array of sciatic nerves of rats at different ages, further in-depth experimental studies, as well as the function and regulatory method of deep mining of differential sequences, are required to fill the gap between age and peripheral nerve structure and function.

Collectively, this study showed that gene expression in sciatic nerves differs between adult and young rats, and these differences occur in genes involved in cell activity, proliferation, differentiation, regeneration, myelination, and injury. These results suggest that, in peripheral nerves, both cells and the microenvironment change with age to influence the function and repair of peripheral nerves. These results thus provide important theoretical bases for further studies on the relationship between peripheral nerves, regeneration, and age.

Additional file: Open peer review reports

Click here for additional data file.

Click here for additional data file.

Table 1

Primer information for the genes that were detected

Table 2

Top 10 mRNA groups that were only expressed in the sciatic nerve of rats in young (1 week old, n = 4) and adult (12 months old, n = 4) Sprague-Dawley rats

Table 3

Top 10 differentially expressed mRNA groups in sciatic nerves that were upregulated or downregulated in young (1 week old, n = 4) and adult (12 months old, n = 4) Sprague-Dawley rats

Table 4

Categories in Gene Ontology (GO) analyses that corresponded to the top 20 groups of differentially expressed mRNAs (DEmRNAs) in the sciatic nerve of Sprague-Dawley rats at the ages of 1 week and 12 months

Table 5

Pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses that corresponded to the top 20 groups of enriched differentially expressed mRNAs (DEmRNAs, age of 12 months vs. age of 1 week)