Comprehensive transcriptional profiling of porcine brain aging

Comprehensive transcriptional profiling of porcine brain aging

Gene 693 (2019) 1–9 Contents lists available at ScienceDirect Gene journal homepage: www.elsevier.com/locate/gene Comprehensive transcriptional pro...

2MB Sizes 1 Downloads 66 Views

Gene 693 (2019) 1–9

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Comprehensive transcriptional profiling of porcine brain aging a,1

a,1

b

c

a

T b

Jianning Chen , Qin Zou , Daojun Lv , Muhammad Ali Raza , Xue Wang , Peilin Li , Yan Chena, Xiaoyu Xia, Anxiang Wena, Li Zhud, Guoqing Tangd, Mingzhou Lid, Xuewei Lid, ⁎ Yanzhi Jianga, a

Department of Zoology, College of Life Science, Sichuan Agricultural University, Ya'an, Sichuan 625014, China Sichuan Weimu Modern Agricultural Science and Technology Co., Ltd., Chengdu, Sichuan 611130, China c Department of Crop Cultivation and Farming System, College of Agronomy, Sichuan Agricultural University, Chengdu, Sichuan 611130, China d Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province, Sichuan Agricultural University, Chengdu, Sichuan 611130, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Aging Brain Pig Transcriptome

The brain as an important organ can be affected largely by aging, and the comprehensive transcriptional underpinnings of brain aging remain poorly understood. Here, we performed a high throughput RNA sequencing to evaluate the expression profiles of messenger RNA (mRNA), long non-coding RNAs (lncRNAs), micro RNAs (miRNAs), and circular RNAs (circRNAs) in porcine brain. We have identified 714 mRNAs, 38lncRNAs, 41miRNAs, and 148circRNAs were age-related genes in the porcine cerebral cortex. The lncRNAs, miRNAs and circRNAs have effect on the age of porcine brain due to the much changes of expression level as noncoding RNAs. The up-regulated genes were significantly enriched for stress response, reproductive regulatory process, immune response and metabolic process, and the down-regulated genes were related to neurologic function, stress response and signaling pathway. The synaptic transmission pathway may be the key role in aging of porcine brain that it was co-enriched for in both differentially expressed mRNAs and lncRNAs. Moreover, some lncRNAs and their target genes were also differentially expressed during brain aging. We further assessed the multi-group cooperative control relationships and constructed circRNA-miRNA co-expression networks during brain aging. We also selected 2 mRNAs, 2 lncRNAs, 2 miRNAs, and 1 circRNAs to perform the q-PCR, and the expression patterns were highly consistent between the two methods confirming the high reproducibility and reliability of the gene expression profiling in our study. In conclusion, our findings will contribute to understand the transcriptional underpinnings of brain aging and provide a foundation for future studies on the molecular mechanisms underlying brain aging.

1. Introduction Aging is a major risk factor for the development of many diseases especially for some prevalent neurodegenerative disorders (Yankner, 2000; Rodwell et al., 2002). The mammalian brain, organization and architecture are more complex compared with other organs and tissues, can be affected largely by aging (Persengiev et al., 2011). At the start of this decade, the brain aging is derived mostly from cross-sectional studies (Raz et al., 2005). At the same time, some studies also uncovered the changes of messenger RNA (mRNA) expression that appeared during normal brain aging using microarrays in many species, not only in humans (Lu et al., 2004) but also in monkeys (Fraser et al., 2005) and mice (Lee et al., 2000; Jiang et al., 2001). Moreover, some

studies also have shown that up-regulated genes were mainly related to stress and immune response, and the down-regulated genes were mainly involved in neuronal functions during the brain aging processes in human and macaque (Lu et al., 2004; Erraji-Benchekroun et al., 2005; Somel et al., 2010). These changes of mRNA expression were attributed to oxidative stress and associated with the accumulation of DNA damage during brain aging (Lu et al., 2004; Somel et al., 2010), and the expression reducing of selectively vulnerable genes involved in neuronal survival, memory, and learning may be also caused by DNA damage (Lu et al., 2004). However, brain development and differentiation are the multistep process (Persengiev et al., 2011), and the aging can not only affect the mRNA expression but also contribute to the changes of noncoding RNAs expression (Chen et al., 2018a), thus



Corresponding author at: Department of Zoology, College of Life Science, Sichuan Agricultural University, No. 46, Xinkang Road, Ya'an City 625014, Sichuan Province, China. E-mail address: [email protected] (Y. Jiang). 1 These authors contributed equally to this work. https://doi.org/10.1016/j.gene.2019.01.019 Received 6 November 2018; Received in revised form 27 December 2018; Accepted 22 January 2019 Available online 26 January 2019 0378-1119/ © 2019 Elsevier B.V. All rights reserved.

Gene 693 (2019) 1–9

J. Chen et al.

2.2. Samples preparation

the transcriptional underpinnings of brain aging remain poorly understood. Several non-coding RNAs, including long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs) have been believed to play an important role in regulating gene expressions (Eulalio et al., 2008; Mercer et al., 2009; Chen and Yang, 2015), and they were showed being closely related to many diseases (Chen et al., 2013; Lukiw, 2013; Feng et al., 2014; Zhao et al., 2016). In recent years, many studies about profiles of comprehensive transcriptome involved in noncoding RNAs during aging in mammals have been reported, these species include human, pig, and mice, and the research object covered stem cell, cardiac and skeletal muscles, and kidney (White et al., 2015; Zhou et al., 2015; Chen et al., 2018a). However, to the best of our knowledge, a very few researchers have studied the changes of comprehensive transcriptome during mammal's brain aging, therefore, the changes of comprehensive transcriptome including mRNA, lncRNA, miRNA, and circRNA need to be revealed during brain aging in mammals. It was reported that pigs have a genomic structure, physiologic and biochemical similarities with human except anatomic structure (Welsh et al., 2009), is more similar than mice (Wernersson et al., 2005), and has been found to be a valuable model for studying human diseases (Lunney, 2007; Zhu et al., 2015). Moreover, many human neurological and neuropsychiatric diseases are inadequately modeled in rodents (Bakken et al., 2016), and suggests that the pig may be an ideal model animal to understand the aging of the human brain at the transcriptional level. Yanan pig is an indigenous Chinese pig which is raised in the hilly area of western Sichuan Province in the past. And the Yanan pig production has been displaced in the last few decades due to the poor growth performance and carcass composition and is in danger of extinction (Jiang et al., 2012). Moreover, the Yanan pig has been included in the National Program for Farm Animal Resources since 2014 (Chen et al., 2018b). Though the number of Yanan pig has grown in recent years, the conservation research of Yanan pig still needs to be strengthened. To better understand the transcriptional regulatory mechanism of the brain aging and provide more reference information for the conservation research of Yanan pig, high-throughput sequencing was used to perform an integrated multi transcriptome-wide profiling (mRNA, miRNA, lncRNA, and circRNA) analysis. This procedure enabled us to identify the numerous age-related genes and unclear multi-group cooperative control networks in aging of brain. Taken together, these data may help to explain brain age-associated changes in the transcriptional pattern during aging.

All samples used in this study were collected according to the guidelines for care and use of experimental animals established by the Ministry of Science and Technology of China. The porcine brain was collected from each pig and cerebral cortex was rapidly separated from each brain, immediately frozen in liquid nitrogen, and stored at −80 °C until RNA extraction. Total RNA was extracted using TRIzol Reagent (Invitrogen, California, USA), and treated with DNase and purified using an RNeasy Mini Kit (Qiagen, California, USA). The concentration and quality of RNA were assessed on the Agilent Bioanalyzer 2100 system (Agilent, California, USA). 2.3. RNA sequencing and data analysis A total of 5 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs, Massachusetts, USA) following the manufacturer's guidelines and library quality was assessed using an Agilent Bioanalyzer 2100 system. After cluster generation, libraries were sequenced on an Illumina HiSeq 4000, and 150 bp paired-end reads were generated. Clean reads were obtained after removing reads containing the adapter, reads containing ploy-N and low quality reads from raw data. Clean reads were aligned to the Ensemble (Susscrofa 10.2) using TopHat2 (v2.0.14) (Kim et al., 2013) with default parameters. StringTie software (Pertea et al., 2015) was used to assemble the mapped reads of each sample, which at least existed in one of both replicates. Transcripts included BLASTed (evalue = 1e–10) to Ensembl, and mapped transcripts were directly described as known lncRNA or mRNA. Here, we used Salmon (v0.6.0) (Patro et al., 2015) to calculate transcripts per million (TPMs) of both lncRNAs and coding genes in each sample. Then, Coding Potential Calculator (CPC) (0.9) (Kong et al., 2007) and Pfam Scan (v1.5) (Punta et al., 2012) were performed to analyze the coding potential of transcripts. Transcripts predicted to have coding potential by either/all of the tow tools were filtered out, and those without coding potential served as our candidate set of novel lncRNAs. The circRNAs were predicted by CIRCexplorer2 (junction reads ≥2) (Zhang et al., 2016). 2.4. Small RNA sequencing and data analysis A total of 5 μg total RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (New England Biolabs, Massachusetts, USA) following the manufacturer's guidelines, and index codes were added to attribute sequences to each sample. Library quality was assessed using an Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina, California, USA) according to the manufacturer's guidelines. After cluster generation, the library preparations were sequenced on an Illumina MiSeq platform, and 50 bp single-end reads were generated. The miRBase21 was used as the reference, and software mirdeep2 (Friedländer et al., 2011) was used to obtain the potential miRNA, and to predict novel miRNA. In detail, we used miRanda (v3.3a) software (Betel et al., 2008) with default parameters and cutoffs (Score S ≥ 140 and Energy E ≤ –20.0) to predict miRNA target. Subsequently, miRNA expression levels were estimated by TPMs by the following criteria: normalization formula: normalized expression = mapped read count / total reads ∗ 1000000.

2. Material and methods 2.1. Animals Four healthy female Yanan pigs, including two 180-day-old young sows and two 8-year-old sows, were used in our study. No direct and collateral blood relationships were been within the last 3 generations among the four pigs. The piglets were weaned at 28 ± 1 day of age. Then a starter diet, including 3.40 Mcal kg−1 metabolizable energy with 20.0% crude protein (11.5 g/kg lysine), was administered from day 30 to day 60 after weaning. From day 61 to day 120, the pigs were fed a diet, containing 14.0 MJ/kg of metabolizable energy consisting of 18.0% crude protein (9.0 g/kg lysine). From day 121, pigs received a diet containing 13.5 MJ/kg of metabolizable energy and 16.0% crude protein (8.0 g/kg lysine). Animals were allowed access to food and water ad libitum and were maintained under the same conditions. The night before slaughtering, pigs were fastened and given 2 h for rest after transportation. Then, pigs were stunned electrically at 90 V and 50 Hz for 10 s and exsanguinated as essential to ameliorate pain.

2.5. Differential expression analysis Differentially expressed mRNAs, lncRNAs, miRNAs, and circRNAs 2

Gene 693 (2019) 1–9

J. Chen et al.

were grouped separately, indicating an age-dependent expression pattern. Furthermore, a closer distance value was found between old individuals compared to that of younger individuals. The Principal Component Analysis (PCA) showed that the mRNAs and lncRNAs transcripts from YB had a high degree of variance between each biological replicate, whereas the OB were clustered (Fig. 2B). We also investigated the four samples used unsupervised Euclidean matrix plots based on miRNAs (Fig. 2C), we observed that a closer distance value between old individuals compared with younger individuals. Moreover, this result was also reflected by the result of PCA analyses for miRNA that there was only a little bit high degree of variance between each biological replicate (Fig. 2D). Thus, the results indicated that the transcriptional levels of mRNA and lncRNA were more variable than that of their younger counterparts in the brain.

were discovered using the edgeR package (Robinson et al., 2010) running in the R programming environment. The edgeR user guide was followed as detailed in the “classic analysis” section. A differential expression q-value < 0.05 and fold change > 2 assigned as differentially expressed in the old group vs. young group. For DECs, we used a standard p-value < 0.05 and fold change > 2, fold change was log2 transformed ratio of average expression between two groups, calculated by edgeR. 2.6. Gene ontology enrichment analysis GO (Gene Ontology) enrichment analysis was implemented by the GOseq R package, GO terms with corrected p-value < 0.05 were considered significantly enriched by differentially expressed genes.

3.2. Differentially expressed mRNAs during brain aging

2.6.1. CircRNA-miRNA co-expression network The circRNA-miRNA co-expression network was constructed based on the correlation analysis between the differentially expressed circRNA and miRNAs. The expressions of differentially expressed circRNAs and miRNAs were analyzed by Pearson's correlation coefficient. The absolute coefficient value of 0.8 between a circRNA and a miRNA was considered relevant for network construction. A p-value of < 0.05 was considered statistically significant.

From our results, though most genes had been shared for both ages during brain development, we also detected a proportion of age-dependent differentially expressed genes (DEGs) (Methods). We first screened DEGs in old brains vs. young brains (OYB) and identified 714 DEGs, of which 646 genes were up-regulated and 68 genes were downregulated (Fig. 3A, Dataset S3). To explore the biological functions of these age-related genes, we performed the Gene ontology (GO) analysis and found that the up-regulated genes were significantly enriched for stress response, reproductive regulatory process, immune response, locomotion and metabolic process, as response to cold (GO: 0009409, p = 0.000842875, q = 1), embryonic process involved in female pregnancy (GO: 0060136, p = 0.000139674, q = 0.634354322), negative regulation of viral genome replication (GO: 0045071, p = 0.00058449, q = 1), locomotion (GO: 0040011, p = 9.73222E–05, q = 0.634354322) and glucose catabolic process (GO: 0006007, p = 0.000986792, q = 1) (Fig. 3C, Dataset S4). For example, the gene SPNS2 (ENSSSCG00000017879) was identified as up-regulated (log FC = 1.605893177) (Dataset S3) during brain aging, and it was GO enriched for locomotion (Dataset S4). Previously, the SPNS2 had been reported that it can regulate heart development and lymphocyte trafficking in zebrafish and mice respectively acting as S1P transporter (Kawahara et al., 2009; Fukuhara et al., 2012). More importantly, it also had been confirmed that the deficiency of SPNS2 could cause earlyonset progressive hearing loss, and the SPNS2 was required for the maintenance of normal endocochlear potential and normal auditory function by playing the supporting role for S1P signaling in hearing (Chen et al., 2014). While, the down-regulated genes were related to neurologic function, stress response, signaling pathway and immune response, as axon ensheathment (GO:0008366, p = 0.000125859, q = 0.449367113), synaptic transmission, dopaminergic (GO:0001963, p = 0.000695973, q = 0.862057371), cellular response to nerve growth factor stimulus (GO:1990090, p = 0.000164905, q = 0.449367113), response to pain (GO:0048265, p = 0.000449156, q = 0.679972649), tachykinin receptor signaling pathway (GO:0007217, p = 0.00023851, q = 0.464242968), G-protein coupled adenosine receptor activity (GO:0001609, p = 0.0000798, q = 0.449367113), and immunoglobulin production (GO:0002377, p = 0.00000158, q = 0.021460161) (Fig. 3C, Dataset S4). Here, we found that the adenosine A(2A) receptor gene (ADORA2A, ENSSSCG00000010054) was significantly down-regulated (log FC = −2.277155734) (Dataset S3), and it was GO enriched for synaptic transmission, dopaminergic (Dataset S4). The ADORA2A gene is closely related to neuropathic diseases, for example, ADORA2A variants may increase the autistic symptoms and anxiety in autism spectrum disorder (Freitag et al., 2010). And some studies also revealed that the genetic variation (rs5751876) in ADORA2A alters age at onset in Huntington's disease (Dhaenens et al., 2009), it was also found closely related to Parkinson's disease (Albasanz et al., 2008) and Alzheimer's disease (Horgusluoglu et al., 2016).

2.7. Quantitative PCR validation cDNA was synthesized using the oligo (dT) and random 6-mer primers provided in the PrimeScript RT Master Mix kit (TaKaRa, Shiga, Japan). q-PCR was performed using the SYBR Premix Ex Taq kit (Takara, Shiga, Japan) on a CFX96 Real-Time PCR detection system (Bio-Rad Laboratories, California, USA.). The q-PCR validation was carried out using three biological replicates. The primer pairs used are presented in Table S3. Three endogenous control genes (porcine GAPDH, ACTB, and U6 snRNA) were used in this assay. The ΔΔCt method was used to determine the expression level of objective mRNAs, miRNAs, lncRNAs, and circ RNAs. 3. Results 3.1. Transcriptome profile of aging brain To observe the extent that aging induces changes of the transcriptome during brain aging, we collected cerebral cortex of porcine brain samples from four pigs including two pigs in young brain (YB) group and two pigs in old brain (OB) group. Then we performed RNAseq (RNA sequencing) to investigate the transcriptome-wide profiling, including mRNA, miRNA, lncRNA, and circRNA. For RNA-seq libraries, we obtained an average of ~69.17 million clean reads from each of the four samples and 69.80%~74.60% of these reads were uniquely aligned to the reference genome Ensemble Susscrofa 10.2 (Table S1). Furthermore, for the small RNA-seq libraries, we obtained approximately 6.38–9.66 million clean reads from each of the four samples and 63.61%–77.32% of these reads were uniquely aligned to the reference genome Ensemble Susscrofa 10.2 (Table S2). In the four samples, we detected a total of 19,591 mRNAs, 4813 lncRNAs, 1217 miRNAs and 9541 circRNAs expressed (Fig. 1, Dataset S1). In addition, 589 miRNAs of the 1217 miRNAs were identified as novel miRNAs (Fig. 1E). It's worth noting that, there are > 80% of these mRNAs, lncRNAs and miRNAs were expressed in both young and old brains, but only 21.9% of these identified circRNAs were expressed in both young and old brains. Moreover, the number of these identified expressed mRNAs, lncRNAs, miRNAs, and circRNAs in the old group were more than those which expressed in the young group (Fig. 1, Dataset S2). The four samples, we used unsupervised Euclidean matrix plots for expressed mRNA and lncRNA transcripts to determine the variance between old and young individuals (Fig. 2A), the young and the old 3

Gene 693 (2019) 1–9

J. Chen et al.

Fig. 1. Global mRNA, lncRNA, miRNA, and circRNA expression pattern across brain samples. (A) mRNA. (B) lncRNA. (C) miRNA. (D) circRNA. (E) Known and novel miRNA. YB, young brain: OB, old brain.

q = 4.66268E–09) (Fig. 3D, Dataset S5). Moreover, the target genes of down-regulated lncRNAs were significantly related to not only neurological functions but also stress response and ion transport, as synaptic transmission (GO:0007268, p = 2.87356E–09, q = 7.83044E–06), neuropeptide signaling pathway (GO:0007218, p = 4.86556E–06, q = 0.005066913), neurological system process (GO:0050877, p = 5.42653E–07, q = 0.000821517), response to pain (GO:0048265, p = 0.000208398, q = 0.118309297) and ion transport (GO:0006811, p = 9.43816E–11, q = 6.42975E–07) (Fig. 3D, Dataset S5). To investigate the co-regulated relationships of lncRNA-mRNA during brain aging, we selected 3 significantly up-regulated lncRNAs (MSTRG.122718, MSTRG.177553, MSTRG.210181; logFC > 8.0) and 3 down-regulated lncRNAs (MSTRG.189460, MSTRG.189455, MSTRG.102658; logFC < −2.0) to screen their predicted target genes, and found some of their target genes were differentially expressed during brain aging (Table 1). Here, we found the target genes of MSTRG.122718, which had the biggest change (logFC = 9.132076302) in expression, were significantly enriched for neurological function,

3.3. Differentially expressed noncoding RNA during brain aging To observe the expression profiles of age-related long noncoding RNA, we also screened the differentially expressed lncRNAs (DELs) during the brain aging. We identified 38 DELs during brain aging, 30 lncRNAs were up-regulated and 8 lncRNAs were down-regulated (Fig. 3B, Dataset S3). In addition, GO analysis was performed to observe the biological functions that the target genes of DELs involved in. The results showed that these target genes of up-regulated lncRNAs were significantly enriched for neurological function and ion transport, as synaptic vesicle exocytosis (GO:0016079, p = 4.03E–10, q = 7.85133E–07), synaptic transmission (GO:0007268, p = 7.83679E–13, q = 3.55921E-09), synapse organization (GO:0050808, p = 7.65355E–09, q = 9.47997E–06), neuron projection (GO:0043005, p = 1.16519E–14, q = 1.58757E-10), axonogenesis (GO:0007409, p = 1.07128E–09, q = 1.62181E–06), regulation of ion transmembrane transport (GO:0034765, p = 1.28188E–09, q = 1.74657E–06) and ion transport (GO:0006811, p = 1.36886E–12, 4

Gene 693 (2019) 1–9

J. Chen et al.

Fig. 2. Variance between old and young individuals. (A) Pearson correlation coefficient heat map of mRNA and lncRNA between YB and OB, values closer to 1 is less variable. (B) Principal Component Analysis (PCA) plot based on normalized expression level (log2 (TPM)) of identified mRNAs and lncRNAs. (C) Pearson correlation coefficient heat map of miRNAs between YB and OB, values closer to 1 are less variable. (D) Principal Component Analysis (PCA) plot based on normalized expression level (log2 (TPM)) of identified miRNAs. YB, young brain: OB, old brain.

(Fig. S1B, Dataset S3). It's worth noting, we identified the miRNA-144 (logFC = 4.28792888) and miRNA-101 (logFC = 1.800018685) were significantly up-regulated in the aging brain (Dataset S3). Studies have shown that miRNA-144 were up-regulated in aging cortex and cerebellum of humans, chimpanzees and rhesus macaques (Persengiev et al., 2011). Moreover, past report has revealed that the miR-144 and -101 were increased in the cerebellum and cortex of SCA1 and Alzheimer's disease in patients relative to healthy aged brains, and miR144 and -101 inhibition increased ATXN1, the disease-causing gene for the SCA1, levels in human cells (Persengiev et al., 2011).

such as the SLC30A3 (related to neuron projection), CPLX1 (related to neurotransmitter transport) and CRTC1 (related to memory). Moreover, the target genes of down-regulated MSTRG.189455 were also significantly enriched for neurological function, such as CALB2 (related to synapse), TAC1 (related to cellular response to nerve growth factor stimulus) and PDYN (related to neuropeptide signaling pathway). However, the most significantly down-regulated lncRNA was MSTRG.189460 (logFC = −3.813550884) and its target gene was enriched for positive regulation of GTPase activity (Table 1). A total of 41 differentially expressed miRNAs (DEMs) were identified in OYB of which 26 miRNAs were up-regulated and 15 miRNAs were down-regulated (Fig. S1A, Dataset S3). Moreover, 148 differentially expressed circRNAs (DECs) were identified in OYB, of which 133 circRNAs were up-regulated and 15 circRNAs were down-regulated

3.4. Construction of the circRNA-miRNA co-expression network The circRNAs could regulate the function of miRNAs acting as 5

Gene 693 (2019) 1–9

J. Chen et al.

(caption on next page) 6

Gene 693 (2019) 1–9

J. Chen et al.

Fig. 3. Differentially expressed mRNAs and lncRNAs during aging. (A) Differentially expressed mRNAs in OYB. (B) Differentially expressed lncRNAs in OYB. (C) Function enrichment analysis for differentially expressed mRNAs in OYB. Only the most enriched (p < 0.05) and meaningful Gene Ontology terms are presented here. (D) Function enrichment analysis for differentially expressed lncRNAs in OYB. OYB, old brains vs. young brains; red dots and blue dots represent, respectively, up-regulated and down-regulated mRNAs during aging. FDR, false discovery rate; FC, fold change. BP, biological process; CC, cellular component; MF, molecular function. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Zhang et al., 2018), which also indicates that the regulatory mechanism is not conservative at circRNA level during aging in mammals. Moreover, we identified the 714 differentially expressed genes during brain aging in our study, the number of age-related genes are more than that of cardiac and skeletal muscle tissues(Chen et al., 2018a), and relationships of the circRNA-miRNA co-expression network were more complex compared with cardiac and skeletal muscle tissues, we reported in our previous study (Chen et al., 2018a). This may be related to complex neural regulation in the brain and that many of differentially expressed genes were related to neurological function. More changes of gene expressions in the brain compared with cardiac and skeletal muscle tissues during aging also reflect that the brain has more complex biological functions than cardiac and skeletal muscle tissues at the molecular level. Our results showed that the up-regulated genes were related to stress and immune response, and the down-regulated genes were mainly involved in neuronal functions. Similar findings were also reported in human and macaque brain aging (Lu et al., 2004; ErrajiBenchekroun et al., 2005; Somel et al., 2010). In addition, we also identified some age-related genes altered during porcine brain aging, such as the ADORA2A gene, miRNA-144 and -101 were found expression changes during aging, and these genes have been shown to be closely related to the neurodegenerative diseases (Albasanz et al., 2008; Dhaenens et al., 2009; Freitag et al., 2010; Persengiev et al., 2011; Horgusluoglu et al., 2016). Previous study reported that pigs have a genomic structure, physiologic and biochemical similarities with human (Welsh et al., 2009), is more similar than mice (Wernersson et al., 2005). Here, we demonstrated the similarities between pig and humans at the comprehensive transcriptome level during brain aging process. A recent study showed a knockin pig model of Huntington's disease was successfully established using CRISPR/Case9 (Yan et al., 2018), so the pig may be an ideal and valuable model for studying human neurodegenerative diseases. In our study, we observed a closer distance value between old individuals compared to younger individuals for expressed mRNA and lncRNA transcripts, which is contrary to what we have observed in the aging of cardiac and skeletal muscle tissues (Chen et al., 2018a), it may be related to gene expressions, both in progenitor cells and maturing neurons, changes more rapidly before birth (Bakken et al., 2016), and in brain many expression changes associated with aging initiate in early childhood, such as down-regulation of genes involved in neural

competing for endogenous RNAs (ceRNAs) suggests that circRNAs and their target miRNAs may be co-expressed in the aging brain (Sanger et al., 1976; Li et al., 2015). Here, we used the Pearson correlation coefficient based on DECs and DEMs data to predict a co-expression network between circRNAs and their target miRNAs (Dou et al., 2016). We constructed a network map, and it contained 11 DEMs, 54 DECs, and 62 relationships (Fig. 4, Dataset S6). Within the network, 4 miRNAs were co-related with simultaneous up-regulation of 50 circRNAs, and 7 miRNAs were co-related with simultaneous down-regulation of 4 circRNAs; miRNA-411, with 43 circRNAs, showing most relations, and circRNA 005253, showing relations with 6 miRNA (Fig. 4, Dataset S6). 3.5. The expression of genes validated by q-PCR Here, we performed q-PCR to validate the high throughput sequencing results. We selected some gene which were found related to some diseases or important biology functions, or randomly selected some highly expressed genes to performed q-PCR to detect these genes expression level in the four samples (two 180-days pigs and two 8-years pigs), and these genes including 2 mRNAs (DVL1 and RGS16), 2 lncRNAs (MSTRG.102415 and MSTRG.102658), 2 miRNAs (miR-101 and miR-144), and 1 circRNAs (circRNA005253) (Fig. S2). The data indicated that the expression patterns were highly consistent between the two methods (Fig. S2). The q-PCR results confirmed the high reproducibility and reliability of the gene expression profiling in our study. 4. Discussion Here, we analyzed the profiles of comprehensive transcriptome expression, including mRNA, lncRNA, miRNA, and circRNA, during porcine brain aging. These results showed that 90.5% of identified mRNAs, 81.5% of identified lncRNAs, and 95.2% identified miRNAs, were co-expressed in young and old brains, and such a high proportion of genes were co-expressed during porcine brain aging indicates the regulatory mechanism may be very conserved at the mRNA, lncRNA and miRNA levels during aging in mammals. But for the circRNA, only 21.9% of identified circRNAs were co-expressed, which indicates that the circRNAs can not only regulate the function of miRNAs acting as ceRNAs, but also other roles working during aging, as they can translate to proteins and play important roles in organisms (Yun et al., 2017;

Table 1 Gene ontology annotations of differentially expressed lncRNA and their differentially expressed targets in OYB. lncRNA OYB MSTRG.122718↑

MSTRG.177553↑

MSTRG.210181↑ MSTRG.189460↓ MSTRG.189455↓

MSTRG.102658↓

Targets (corrected p-value)

Pearson correlation

Adjusted p-value

Gene ontology annotations

SLC30A3 (ENSSSCG00000021534)↑ CRTC1 (ENSSSCG00000013915)↑ CPLX1 (ENSSSCG00000029481)↑ OLIG1 (ENSSSCG00000012040) ↑ BFSP2 (ENSSSCG00000028736)↑ CDK5R2 (ENSSSCG00000021584)↑ ENSSSCG00000030009↑ RGS16 (ENSSSCG00000015550)↓ CALB2 (ENSSSCG00000002730)↓ TAC1 (ENSSSCG00000023592)↓ PDYN (ENSSSCG00000007180)↓ DMRT-1 (ENSSSCG00000005236)↓ DDX3X (ENSSSCG00000026430)↓

0.947507269 0.94232565 0.930057161 0.937530875 0.943278611 0.92945121 0.944636825 0.981200354 0.949631344 0.916870278 0.933030538 0.946186716 0.943408411

7.96E–06 3.54E–05 0.000676104 0.000122258 2.73E–05 0.000768796 1.86E–05 6.67E–14 4.09E–06 0.008377081 0.000352436 1.18E–05 2.63E–05

Neuron projection Memory Neurotransmitter transport Neuron fate commitment Lens fiber cell development Hippocampus development Protein binding Positive regulation of GTPase activity Synapse Cellular response to nerve growth factor stimulus Neuropeptide signaling pathway Developmental process involved in reproduction ATP-dependent RNA helicase activity

Note: abbreviations: OYB: old brain vs. young brain; ↑: up-regulated; ↓: down-regulated. 7

Gene 693 (2019) 1–9

J. Chen et al.

Fig. 4. Construction of the circRNA-miRNA co-expression network. Construction of the circRNA-miRNA co-expression network in the brain. Red color and blue color represent up and down-regulation, respectively. The black solid line represents the positive correlation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

functions (Somel et al., 2010). Study also shown approximately 50% of the total expression changes observed between newborns and 40-yearolds occurred within the first year of life in human brain aging (Somel et al., 2009). These results suggest us that alters of the age-related transcriptome of the brain may predate other tissues of the body in big mammals. In addition, we also observed that there was only a little bit high degree of variance between each biological replicate from expressed miRNA. And a past study has shown the miRNA expression in the frontal cortex appeared to be less affected by aging and showed more uniform patterns when young and old subjects were compared (Persengiev et al., 2011). These suggest the regulation functions of miRNA are more conserved than mRNA during brain aging.

locomotion and metabolic process, and the down-regulated genes were related to neurologic function, stress response and signaling pathway. Moreover, the circRNA-miRNA co-expression networks during brain aging have been constructed. In brief, our findings will contribute to understand the transcriptional underpinnings of brain aging and provide a foundation for future studies on the molecular mechanisms underlying brain aging.

5. Conclusions

Acknowledgments

We performed a high throughput RNA sequencing to evaluate the expression profiles of mRNA, lncRNAs, miRNAs, and circRNAs. There were numerous age-related genes have been identified in the porcine cerebral cortex. The up-regulated genes were significantly enriched for stress response, reproductive regulatory process, immune response,

This study was supported by the National Key R&D Program of China (2018YFD0501204), Key Program of National Science Foundation of China (Grant No. 31530073), and Chengdu Local Good Varieties of Livestock and Poultry Resources Protection and Exploitation and Utilization of Construction Projects.

Conflict of interest statement The authors declare that they have no competing interests.

8

Gene 693 (2019) 1–9

J. Chen et al.

Availability of data and materials

in migration of zebrafish myocardial precursors. Science 323, 524–527. Kim, D., Pertea, G., Trapnell, C., et al., 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14, R36. Kong, L., Zhang, Y., Ye, Z.-Q., et al., 2007. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. Lee, C.K., Weindruch, R., Prolla, T.A., 2000. Gene-expression profile of the ageing brain in mice. Nat. Genet. 25, 294–297. Li, J., Yang, J., Zhou, P., et al., 2015. Circular RNAs in cancer: novel insights into origins, properties, functions and implications. Am. J. Cancer Res. 5, 472. Lu, T., Pan, Y., Kao, S.-Y., et al., 2004. Gene regulation and DNA damage in the ageing human brain. Nature 429, 883–891. Lukiw, W., 2013. Circular RNA (circRNA) in Alzheimer's disease (AD). Front. Genet. 4, 307. Lunney, J.K., 2007. Advances in swine biomedical model genomics. Int. J. Biol. Sci. 3, 179–184. Mercer, T.R., Dinger, M.E., Mattick, J.S., 2009. Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159. Patro, R., Duggal, G., Kingsford, C., 2015. Salmon: Accurate, Versatile and Ultrafast Quantification from RNA-seq Data Using Lightweight-alignment. bioRxiv 021592. Persengiev, S., Kondova, I., Otting, N., et al., 2011. Genome-wide analysis of miRNA expression reveals a potential role for miR-144 in brain aging and spinocerebellar ataxia pathogenesis. Neurobiol. Aging 32, 2316.e17–2316.e27. Pertea, M., Pertea, G.M., Antonescu, C.M., et al., 2015. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. Punta, M., Coggill, P., Eberhardt, R., et al., 2012. The Pfam protein families database. Nucleic Acids Res. 40, D290–D301. Raz, N., Lindenberger, U., Rodrigue, K.M., et al., 2005. Regional brain changes in aging healthy adults: general trends, individual differences and modifiers. Cereb. Cortex 15, 1676–1689. Robinson, M.D., McCarthy, D.J., Smyth, G.K., 2010. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. Rodwell, G.E.J., Sonu, R., Zahn, J.M., et al., 2002. A transcriptional profile of aging in the human kidney. PLoS Biol. 12, 1566–1573. Sanger, H.L., Klotz, G., Riesner, D., et al., 1976. Viroids are single-stranded covalently closed circular RNA molecules existing as highly base-paired rod-like structures. Proc. Natl. Acad. Sci. U. S. A. 73, 3852. Somel, M., Franz, H., Yan, Z., et al., 2009. Transcriptional neoteny in the human brain. Proc. Natl. Acad. Sci. U. S. A. 106, 5743–5748. Somel, M., Guo, S., Fu, N., et al., 2010. MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain. Genome Res. 20, 1207. Welsh, M.J., Rogers, C.S., Stoltz, D.A., et al., 2009. Development of a porcine model of cystic fibrosis. Trans. Am. Clin. Climatol. Assoc. 120, 149. Wernersson, R., Schierup, M.H., Jørgensen, F.G., et al., 2005. Pigs in sequence space: a 0.66X coverage pig genome survey based on shotgun sequencing. BMC Genomics 6, 70. White, R.R., Milholland, B., Macrae, S.L., et al., 2015. Comprehensive transcriptional landscape of aging mouse liver. BMC Genomics 16, 1–15. Yan, S., Tu, Z., Liu, Z., et al., 2018. A Huntingtin Knockin pig model recapitulates features of selective neurodegeneration in Huntington's disease. Cell 173, 989–1002 e13. Yankner, B.A., 2000. A century of cognitive decline. Nature 404, 125. Yun, Y., Fan, X., Mao, M., et al., 2017. Extensive translation of circular RNAs driven by N6-methyladenosine. Cell Res. 27, 626. Zhang, X.-O., Dong, R., Zhang, Y., et al., 2016. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. In: Genome Res. Zhang, M., Huang, N., Yang, X., et al., 2018. A novel protein encoded by the circular form of the SHPRH gene suppresses glioma tumorigenesis. Oncogene 37. Zhao, Y., Alexandrov, P.N., Jaber, V., et al., 2016. Deficiency in the ubiquitin conjugating enzyme UBE2A in Alzheimer's disease (AD) is linked to deficits in a natural circular miRNA-7 sponge (circRNA; ciRS-7). Gene 7, 116. Zhou, Z., Gao, M., Liu, Q., et al., 2015. Comprehensive transcriptome analysis of mesenchymal stem cells in elderly patients with osteoporosis. Aging Clin. Exp. Res. 27, 595–601. Zhu, J., Chen, C., Yang, B., et al., 2015. A systems genetics study of swine illustrates mechanisms underlying human phenotypic traits. BMC Genomics 16, 88.

All data generated in this study are included in the main article and its supplementary files. Requests for raw data should be made to the corresponding authors. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.gene.2019.01.019. References Albasanz, J.L., Perez, S., Barrachina, M., et al., 2008. Adenosine receptors are increased and sensitized in the frontal cortex from Alzheimer's disease cases. Purinergic Signal 4, S32. Bakken, T.E., Miller, J.A., Ding, S.L., et al., 2016. Comprehensive transcriptional map of primate brain development. Nature 535, 367. Betel, D., Wilson, M., Gabow, A., et al., 2008. The microRNA. org resource: targets and expression. Nucleic Acids Res. 36, D149–D153. Chen, L.L., Yang, L., 2015. Regulation of circRNA biogenesis. RNA Biol. 12, 381–388. Chen, G., Wang, Z., Wang, D., et al., 2013. LncRNADisease: a database for long-noncoding RNA-associated diseases. Nucleic Acids Res. 41, D983–D986. Chen, J., Ingham, N., Kelly, J., et al., 2014. Spinster homolog 2 (Spns2) deficiency causes early onset progressive hearing loss. PLoS Genet. 10, e1004688. Chen, J., Zou, Q., Lv, D., et al., 2018a. Comprehensive transcriptional landscape of porcine cardiac and skeletal muscles reveals differences of aging. Oncotarget 9, 1524–1541. Chen, Y., Wei, Y., Chen, J., et al., 2018b. Growth, carcass characteristics and meat quality of Chinese indigenous Yanan pig crossbred with Duroc and Berkshire genotypes. Anim. Prod. Sci. https://doi.org/10.1071/AN17450. Dhaenens, C.M., Burnouf, S., Simonin, C., et al., 2009. A genetic variation in the ADORA2A gene modifies age at onset in Huntington's disease. Neurobiol. Dis. 35, 474–476. Dou, C., Cao, Z., Yang, B., et al., 2016. Changing expression profiles of lncRNAs, mRNAs, circRNAs and miRNAs during osteoclastogenesis. Sci. Rep. 6, 21499. Erraji-Benchekroun, L., Underwood, M.D., Arango, V., et al., 2005. Molecular aging in human prefrontal cortex is selective and continuous throughout adult life. Biol. Psychiatry 57, 549–558. Eulalio, A., Huntzinger, E., Izaurralde, E., 2008. Getting to the root of miRNA-mediated gene silencing. Cell 132, 9–14. Feng, X., Wang, Z., Fillmore, R., et al., 2014. MiR-200, a new star miRNA in human cancer. Cancer Lett. 344, 166–173. Fraser, H.B., Khaitovich, P., Plotkin, J.B., et al., 2005. Aging and gene expression in the primate brain. PLoS Biol. 3, e274. Freitag, C.M., Agelopoulos, K.E., Rothermundt, M., et al., 2010. Adenosine A(2A) receptor gene (ADORA2A) variants may increase autistic symptoms and anxiety in autism spectrum disorder. Eur. Child Adolesc. Psychiatry 19, 67–74. Friedländer, M.R., Mackowiak, S.D., Li, N., et al., 2011. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. Fukuhara, S., Simmons, S., Kawamura, S., et al., 2012. The sphingosine‑1‑phosphate transporter Spns2 expressed on endothelial cells regulates lymphocyte trafficking in mice. J. Clin. Invest. 122, 1416–1426. Horgusluoglu, E., Nho, K., Risacher, S.L., et al., 2016. A meta-analysis identifies ADORA2A associated with hippocampal volume in Alzheimer's disease. Alzheimers Dement. 12, P636. Jiang, C.H., Tsien, J.Z., Schultz, P.G., et al., 2001. The effects of aging on gene expression in the hypothalamus and cortex of mice. Proc. Natl. Acad. Sci. U. S. A. 98, 1930–1934. Jiang, Y.Z., Zhu, L., Li, F.Q., et al., 2012. Carcass composition and meat quality of indigenous Yanan pigs of China. Genet. Mol. Res. 11, 166–173. Kawahara, A., Nishi, T., Yu, H., et al., 2009. The sphingolipid transporter Spns2 functions

9