Digital Gene Expression Profiling Analysis of Aged Mice under Moxibustion Treatment

Aging is closely connected with death, progressive physiological decline, and increased risk of diseases, such as cancer, arteriosclerosis, heart disease, hypertension, and neurodegenerative diseases. It is reported that moxibustion can treat more than 300 kinds of diseases including aging related problems and can improve immune function and physiological functions. The digital gene expression profiling of aged mice with or without moxibustion treatment was investigated and the mechanisms of moxibustion in aged mice were speculated by gene ontology and pathway analysis in the study. Almost 145 million raw reads were obtained by digital gene expression analysis and about 140 million (96.55%) were clean reads. Five differentially expressed genes with an adjusted P value < 0.05 and |log⁡2(fold  change)| > 1 were identified between the control and moxibustion groups. They were Gm6563, Gm8116, Rps26-ps1, Nat8f4, and Igkv3-12. Gene ontology analysis was carried out by the GOseq R package and functional annotations of the differentially expressed genes related to translation, mRNA export from nucleus, mRNA transport, nuclear body, acetyltransferase activity, and so on. Kyoto Encyclopedia of Genes and Genomes database was used for pathway analysis and ribosome was the most significantly enriched pathway term.


Introduction
Aging is closely connected with death, progressive physiological decline, and increased risk of diseases, such as cancer, type 2 diabetes mellitus, all forms of arteriosclerosis, heart disease, hypertension, age-related macular degeneration, and neurodegenerative diseases [1]. Aging is a complex biological process and the process is modulated by multiple intrinsic and extrinsic factors such as hereditary, environmental, dietary, lifestyle, and stochastic factors throughout a person's lifetime [2]. Epigenetic alterations during aging are associated with modifications of DNA and histone proteins which straight affect chromatin structure and then affect gene expression and genomic stability [3]. The damage of various cell components such as organelles, especially protein and DNA damage, has been proposed as a factor causing aging on the cellular level [4]. Many changes of gene expression in the process of aging have been examined by some surveys based on the central dogma. And a number of pathways including the insulin/IGF-1 and TOR pathways and key factors including sirtuins and AMP kinase have been identified to regulate the aging process [5].
Moxibustion has been used for thousands of years in China and other Asian countries. Generally, moxibustion exerts its healing effects by burning the moxa herb (Artemisia vulgaris) on or above the skin at acupuncture points to warm the points [6]. Current studies have confirmed that moxibustion can efficiently control inflammation in bowel by 2 Evidence-Based Complementary and Alternative Medicine regulating physiological balance at multiple links and multiple targets in the body with the low adverse reactions, low recurrence, and good long-term efficacy. Therefore, moxibustion has been applied diffusely [7]. It is reported that moxibustion can treat more than 300 kinds of diseases including malposition, diarrhea, colitis, urinary incontinence, dysmenorrhea, knee osteoarthritis, temporomandibular joint disturbance syndrome, soft tissue injury, heel pain, asthma, urinary retention, herpes zoster, and aging related problems [8]. In addition, moxibustion can improve immune function and physiological functions including cerebral blood flow in animals, which were suggested by the preclinical studies [9]. These findings exposed the potential value of moxibustion for giving humans some padding against aging in experimental research as well as practical application.
Digital gene expression (DGE) is a tag-based transcriptome sequencing method. DGE can be used to analyze quantitative gene expression and compare expression profiles without potential deviations in order to achieve sensitive and accurate transcriptome analysis [13]. DGE technology has been used more commonly in examining the differences of transcriptional responses in different tissues and organs. This is because DGE can identify millions of differentially expressed genes before prior annotation and the organisms without much genetic information are also allowed to be analyzed by DGE [14].
In the present study, RNA sequencing technology was used to evaluate the mechanisms that moxibustion functioned on aged mice. Based on the theory of traditional Chinese medicine about "the spleen being the foundation of acquired constitution," only spleen samples were examined in DGE analysis. Differentially expressed genes in response to moxibustion were identified by assembling and annotating the transcriptome sequences which were identified in the spleen samples and analyzing the gene expression profiles.

Animals.
Sixteen Institute of Cancer Research (ICR) female mice (19 months old) weighting 45 ± 5 g were purchased from the Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China, Certificate number: SCXK (Beijing) 2012-0001). All the mice were housed in standard laboratory conditions (22 ∘ C ± 2 ∘ C, 50-60% indoor humidity) under a 12 h light-dark cycle (lights on from 8:00 am to 8:00 pm) and fed common chow. The experimental procedures were performed in strict accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. And all animal experiments in this study were approved by the Animal Ethics Committee of Beijing University of Chinese Medicine.

Moxibustion Intervention.
Sixteen aged mice were randomly divided into the control and moxibustion groups ( = 8 mice per group). The normal mice were used as control group. The Guanyuan acupoint (RN4) was selected for treatment with moxibustion. The used moxa sticks (diameter 0.5 cm, length 20 cm) were provided by Nanyang Hanyi Moxibustion Technology Development Co., Ltd. (Henan, China). The mice were immobilized in fixators and the moxa stick was ignited and placed 1 cm above RN4 acupoint for 10 min to perform moxibustion. The mice in control group were immobilized for 10 min using the same method without moxibustion intervention. The treatment of moxibustion was carried out once in two days for a total of 60 days.

Sample Collection.
After 60 days of moxibustion treatment, 6 healthier mice were selected from two groups ( = 3 mice per group). They were anesthetized with chloral hydrate and spleen specimens were collected. The obtained spleen samples were rapidly frozen in liquid nitrogen and stored at −80 ∘ C until RNA extraction.

RNA Extraction and Qualification.
Total RNA was isolated from the spleen samples using Trizol reagent (Invitrogen, Canada) according to the protocol provided by manufacturer. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked with a NanoPhotometer5 spectrophotometer (IMPLEN, USA). RNA concentration was determined by a Qubit5 RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, USA) and RNA integrity was evaluated using RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, USA).

Library Preparation for DGE Sequencing.
Sequencing libraries were produced using NEBNext5 Ultra6 RNA Library Prep Kit (Illumina, USA) according to the recommendations of manufacturer. Three micrograms of RNA was used as the input material to sequence. Briefly, mRNA was purified from total RNA and fragmentation was carried out. Then first strand cDNA was synthesized with random hexamer primer and M-MuLV Reverse Transcriptase and second strand cDNA was subsequently synthesized. Following this, NEBNext Adaptor with hairpin loop structure was ligated. cDNA fragments of preferentially 150∼200 bp in length were selected by purifying the library fragments using AMPure XP system (Beckman Coulter, USA). Then PCR was carried out and library quality was evaluated on the Agilent Bioanalyzer 2100 system.

Quality Control and Mapping
Analyses. Raw reads of FASTQ format were processed with in-house Perl scripts. In this step, reads containing adapter sequences and ploy-N and low quality reads were removed from raw data to clean data, followed by calculating the Q20, Q30, and GC content of the clean data. The clean data with high quality were used for all the downstream analyses. Bowtie v2.0.6 was used to build index of the reference genome and single-end clean reads were aligned to the reference genome through TopHat v2.0.9 [15]. TopHat can create a database of splice junctions on the basis of the gene model annotation file and thus generate a better mapping result than other nonsplice mapping tools. Therefore, TopHat was selected as the mapping tool.

Quantification of Gene Expression Level and Differential
Expression Analysis. The number of reads mapped to each gene was counted using HTSeq v0.6.1 [16]. Then each gene's FPKM (fragments per kilobase of exon model per million mapped reads) value was calculated on the basis of the length of the gene and read count mapped to this gene. The effect of sequencing depth and gene length for the reads count was considered simultaneously by FPKM which is more accurate than RPKM (reads per kilobase of exon model per million mapped reads) in data processing of paired-end sequencing.
Before analysis of differential gene expression, edgeR program package was used to adjust the read counts for each sequenced library through one scaling normalized factor. The DESeq R package (1.12.0) was used to carry out the differential expression analysis of two groups [17]. There were three biological replicates per groups. The resulting values were adjusted according to the Benjamini-Hochberg procedure to control the false discovery rate [18]. Genes with an adjusted value < 0.05 and | log 2(fold change)| > 1 found by DESeq were considered to be of differential expression.

GO and Pathway Analysis.
In order to investigate the biological process, cellular component, and molecular function of differentially expressed genes, GO enrichment analysis was carried out by the GOseq R package [19]. Functional categories were enriched within genes and GO terms with value < 0.05 were considered significantly enriched by differential expressed genes.
Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used for pathway analysis of differential expression genes by the KOBAS software [20]. KEGG is the most commonly used database resource to understand high-level functions and utilities of the biological system, including the cell, the organism, and the ecosystem (http://www.genome.jp/kegg/).

Analysis of DGE Libraries.
Six DGE libraries were sequenced from 6 independent samples of female mice with or without moxibustion treatment to constitute a complete quantitative and qualitative gene expression response profile. Three individual samples marked as C1, C2, and C3 were included for the control group and another 3 individual samples marked as M1, M2, and M3 were included for the moxibustion group. The main characteristics of libraries were summarized and presented in Table 1. Approximately 26.83, 22.68, 21.59, 20.01, 25.58, and 27.05 million raw reads were obtained from the C1, C2, C3, M1, M2, and M3 libraries, respectively. And more than 96% of raw reads were clean reads in each library.

Mapping
Reads to the Transcriptome. The clean reads from the six libraries were aligned to the reference genome to establish gene expression profiles. More than 84% of the clean reads mapped to the reference genome (Table 2). Less than 6% of the clean reads were multiple mapped reads which primarily came from the rRNA and intergenic regions [21].

Analysis of Differential Gene
Expression. The number of unambiguous clean tags for each gene was calculated in different FPKM value intervals for analysis of gene expression. Results showed that most FPKM of reads were between 0 and 1 ( Table 3). Both of control group and moxibustion group showed significant and high correlation between any two replicates. The correlation coefficients were 0.912, 0.906,   Figure 1. The logarithmic FPKM +1 values were assigned as coordinate values of two axes and all data points exhibited a distribution in the region of the diagonal. An adjusted value < 0.05 and | log 2(fold change)| > 1 were chosen as the cutoff criteria to identify differentially expressed genes. Consequently, only 5 differentially expressed genes were eligible, including 4 upregulated and 1 downregulated genes (Table 4).

GO and KEGG Enrichment Analysis of Differentially
Expressed Genes. Gene ontology of differentially expressed genes was analyzed based on biological process, cellular component, and molecular function. Twenty-two GO terms were significantly regulated by the expressed genes. The highly enriched GO terms in biological process, cellular component, and molecular function included translation, mRNA export from nucleus, mRNA transport, transcription export complex, nuclear speck, nuclear body, acetyltransferase activity, N-acyltransferase activity, and structural constituent of ribosome. The GO analysis results of differentially expressed genes are shown in Figure 2. The pathway significantly (corrected value < 0.05) affected by differentially expressed genes was enriched through KEGG pathway analysis. Ribosome (corrected value = 0.0059) was the most significantly enriched KEGG pathway term (Figure 3).

Discussion
In the present study, the female mice were treated with moxibustion and transcriptome analysis was carried out by using digital gene expression profiling. We investigated the influence of moxibustion on aged mice. It is believed that moxibustion can increase the generation of protective proteins including heat shock protein 70 (HSP70) and activate heat-sensitive neural release of nitric oxide to exert its efficacy [22]. According to biomedical theories, various pathways including elimination of free radical damage and regulation of immunity, neuroendocrine system, lipid metabolism, rheological characteristics of blood flow, trace elements, telomerase, and carbonyl poisoning are supposed to be the mechanisms by which moxibustion impact aging processes [23]. However, the mechanism by which moxibustion affects the genes expression of aged mice was still unknown. Therefore, transcriptome sequencing was carried out and compared for each sample.
The transcriptomes of spleen samples obtained from aged ICR mice given moxibustion for 60 days were sequenced to reveal the molecular mechanisms of moxibustion. In total, almost 145 million raw reads were obtained. Among them, approximately 140 million (96.55%) were clean reads. Five differentially expressed genes with an adjusted value < 0.05 and | log 2(fold change)| > 1 were identified between the control and moxibustion groups. They were Gm6563, Gm8116, Rps26-ps1, Nat8f4, and Igkv3-12.
In the light of GO enrichment analysis, the differentially expressed genes were enriched in molecular function, cellular components, and biological processes. According to the GO classifications, the differentially expressed genes were involved in structural constituent of ribosome, nucleus, mRNA export, mRNA transport, and translation. Results suggested that moxibustion plays a role in aged mice by influencing molecular function including acetyltransferase activity, n-acyltransferase activity, and structural constituent of ribosome. Moxibustion also affects cellular component including transcription export complex, nuclear speck, and nuclear body and biological process, such as translation, mRNA transport, and nuclear export. MRNA is released from nuclei and the number obviously decreases with age, which results from an impairment of polyadenylation of mRNA, hnRNA processing, release of mRNA from nuclear matrix, and translocations of mRNA from nuclear to cytoplasmic compartment [24]. Therefore, the level of mRNA released from nuclei is low in aged mice. But, this situation might be improved by moxibustion. It has been found that normal aging is connected with the exceptional enhancement in the generation and maturation of many mRNAs in human peripheral blood leukocytes [25]. A good deal of evidence proves that mRNA translation plays an important role in modulating aging. Studies have shown that lifespan can be increased by interventions resulting in decreased mRNA translation. Besides that, mRNA translation also affects various aging related processes and impacts on multiple longevity pathways [26]. Although the potential mechanisms for such regulation are still a mystery, it is definitely important to give the strong connection between mRNA translation and aging in different species [27]. In the present study, the functions of moxibustion target involved mRNA translation, which suggested that moxibustion might affect mRNA translation of aged mice. In addition, nuclear body, nuclear speck, and transcription export complex were the results of GO enrichment analysis in cellular component. Intriguingly, some evidence shows that nucleoli have additional functions connected with the cell cycle, cellular aging, signal  recognition particle biosynthesis, small RNA processing, and mRNA transport [28]. Cellular aging is generally related to the instability of the nuclear and mitochondrial genomes and oxidative protein damage. That may directly lead to changes in transcription factors and chromatin structure and function. These changes gradually occur with aging and make for alterations of nuclear gene expression, mRNA stability, or both [29]. The finding of lamin A-dependent nuclear defects in human physiological aging further evidenced the correlation between nuclear envelope alterations and the process of aging [30]. Therefore, moxibustion may have a potential activity in antiaging because of intervention in molecular function, cellular components, and biological processes.
KEGG pathway analysis of the 5 differentially expressed genes was completed and result exhibited that a KEGG pathway named ribosome was the most significantly enriched.   The ribosome, the integral but primarily passive participant in the synthesis of proteins across all kingdoms of life, is one of life's most ancient molecular machines. The ribosome has been considered as the decoder of the genome with high precision in the flow of biological information from mRNA to protein [31,32]. It has been confirmed that the ribosome is huge machinery referred to translation of the genetic code into proteins and is one of the major targets of inhibitors for protein synthesis [33]. It revealed that moxibustion might function on aged mice by affecting protein synthesis through ribosome. There is increasing evidence revealing that the human pathological conditions characterized by an upregulated ribosome biogenesis can increase the risk of tumorigenesis [34]. Moxibustion may decrease the risk of tumorigenesis in aged mice by inhibiting the upregulation of ribosome biogenesis. Disturbances of ribosome homeostasis and hematopoietic dysfunction can be disproportionately caused by haploinsufficiencies of specific ribosomal proteins and defects in ribosome biogenesis [35].
The changes of morphology and function in the nucleolus were widely observed in cancer tissues, which have been indicated by the available evidence. These changes were caused by both the raised demand for ribosome biogenesis (the characteristic of proliferating cells) and the changes in the mechanisms which control cell proliferation [36]. The mRNA that carries the message of gene from the DNA leaves the nucleus to move to a ribosome in production of the particular protein. The ribosome moves along the mRNA for translating a ribonucleotide sequence into a corresponding amino acid sequence that constitutes the protein. Therefore, moxibustion may intervene in mRNA export and transport by ribosome pathway. According to the results of GO and KEGG enrichment analysis of differentially expressed genes, it was speculated that moxibustion functioned on aged mice by intervening in the biological process, such as translation, mRNA export from nucleus, and mRNA transport and pathways involved in the ribosome.

Conclusions
In the present study, we investigated the digital gene expression profiling of aged mice with or without moxibustion treatment and speculated the mechanisms of moxibustion in aged mice by analyzing the biological properties, molecular functions, and enriched pathways of differentially expressed genes. There were 5 significantly differentially expressed genes, namely, 4 upregulated genes and 1 downregulated gene. GO enrichment analysis exhibited that these genes were connected with structural constituent of ribosome, nucleus, mRNA transport, and translation. KEGG pathway analysis showed that they were enriched in the ribosome. According to the results of GO and KEGG enrichment analysis and considering the function of the ribosome, it was speculated that interventions to some biological process, such as translation, mRNA export from nucleus, and mRNA transport and pathways involved in the ribosome might be the mechanisms by which moxibustion functioned on aged mice.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.

Authors' Contributions
Nan Liu and Yunyao Jiang contributed equally to this paper.