Dynamic gene expression profiles during postnatal development of porcine subcutaneous adipose

Ethics statement

All research involving animals was conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and approved by the Institutional Animal Care and Use Committee in College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China under permit No. DKY-B20110801. Animals were fed according to the nutrition standards by National Research Council (NRC, 1998) and the feeding conditions (including temperature, humidity etc.) were strictly controlled. All animals were fasted 24 h prior to slaughter and were slaughtered simultaneously, with the exception of newborn animals which were not fasted before slaughter.

Animals and tissues collection

A total of 12 healthy female Jinhua pigs (a fatty Chinese native breed) were used in this study from four postnatal developmental stages: 0 d (birth), 30 d, 180 d, and 7 y. Each stage included three individuals, which were regarded as biological replicates. Venous blood (50 ml) was collected in pre-chilled tubes from each pig immediately before sacrifice. After sacrifice, the HLB, Abdominal Subcutaneous Adipose (ASA), Retroperitoneal Adipose (RAD), Mesenteric Adipose (MAD), Greater Omentum (GOM), Longissimus Dorsi Muscle (LDM), Psoas Major Muscle (PMM), Cardiac Muscle (CM), liver, spleen, lung and brain were rapidly separated from each carcass, immediately frozen in liquid nitrogen, and stored at −80 °C until RNA and DNA extraction.

Measurement of adipose-related phenotype

Measurements of concentrations of 8 serum-circulating indicators of metabolism and adipocyte volume are from our previous report (Li et al., 2012). Serum concentrations of Total Cholesterol (TC), Triglycerides (TG), High Density Lipoprotein (HDL), Low Density Lipoprotein (LDL), Very-Low Density Lipoprotein (VLDL), Lipoprotein a (Lip-a), Apolipoprotein A1 (Apo-A1) and Apolipoprotein B (Apo-B) were determined by using CL-8000 clinical chemical analyzer (Shimadzu, Kyoto, Japan) via standard enzymatic procedures. The adipocyte volume were measured using Hematoxylin-Eosin (H&E) staining method. The mean diameter of an adipocyte was calculated as the geometric average of the maximum and minimum diameter, and 100 cells were measured for each sample in randomly selected fields. The mean adipocyte Volume (V) was obtained according to the following formula: V = π/6 Σfi Di³/Σfi, where Di is the mean diameter; fi denotes number of cells with that mean diameter Di.

RNA and DNA extraction

Total RNA was extracted from frozen tissues using TRIzol (Invitrogen, Carlsbad, CA, USA) and further purified using an RNeasy column (Qiagen, Valencia, CA, USA) according to the manufacturer’s protocol. RNA integrity and concentration were analysis with the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). DNA was isolated from frozen tissues using the DNeasy Blood & Tissue Kit (Qiagen) follow the manufacturer’s instructions.

Digital Gene Expression (DGE) library preparation and sequencing

The DGE libraries preparation for HLB was carried out using Illumina Gene Expression Sample Prep Kit according to the manufacturer’s protocol. In brief, mRNA was isolated from 4 μg of total RNA by binding the mRNA to a magnetic Oligo dT beads. Double strand cDNA were synthesized using Oligo dT beads. Subsequently, the cDNA samples were digested using the restriction enzyme NlaIII, which recognizes and cuts the most 3′ “CATG.” Once digested into fragments, cDNA were ligated to Illumina adapter 1, which contains a recognition site for enzyme MmeI, and the enzyme MmeI was used to create the 17 bp tag. The Illumina adapter 2 was ligated after MmeI digestion and obtains 21 bp tag derived from a single transcript with Illumina adapters attached to both ends. A 15 cycle PCR was performed with two primers (primer GX1 and primer GX2) that anneal to the ends of the adapters to enrich the adapter-ligated cDNA construct. The amplified cDNA construct was purified from a 6% TBE PAGE gel. During the QC steps, Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System are used in quantification and qualification of the sample library. Finally, massively parallel sequencing by synthesis was performed on the Illumina HiSeq™ 2000. Image recognition and base calling were performed using the Illumina Pipeline.

Analysis of DGE libraries

The raw data were filtered after data processing. The filtering pipeline is, as follows: 1) empty reads removal (reads with only 3′ adaptor sequences but no tags, 2) low quality tag removal (tags with unknown sequences, N), 3) removal of tags too long or too short, leaving only those 21 nt in length; 4) removal of tags with a copy number of 1 (probably resulting from a sequence error), and 5) generation of clean tags. The raw datasets have been submitted to NCBI Gene Expression Omnibus database with the accession number GSE46755. All possible CAGT+17 nucleotide tags were created by using sus scrofa UniGene from Ensembl. All clean tags were mapped to the reference sequences (Sscrofa 10.2) and only 1 bp mismatch was allowed. The numbers of mapped clean tags was calculated for each library and were then normalized to Transcripts Per Million tags (TPM). To identify DE genes (P < 0.01) for the clustering analysis, we used one-way repeated-measures ANOVA for comparisons. Resulting P-values of above tests were corrected with adjusted Bonferroni method (FDR < 0.01, 1,000 permutations).

Measurement of mitochondrion DNA (mtDNA) content

The relative mtDNA copy number was determined by quantitative PCR (q-PCR) with primers for the mitochondrial ATP synthase 6 (ATP6), Cytochrome c Oxidase 1 (COX1) and NADH Dehydrogenase 1 (ND1) genes. The nuclear encoding gene Glucagon (GCG) was simultaneously used as an endogenous control gene for normalization (see Table S1 for primer sequences). All reactions were performed in triplicates, and negative controls (without template) were always included. Relative mtDNA copy number per diploid cell = (No. of copies of the mtDNA gene)/(No. of copies of GCG), data are expressed as mean ± SD.

Analysis of development-related expression patterns

Gene expression patterns during the investigated four growth and development stages were identified using STEM software (Ernst, Nau & Bar-Joseph, 2005; Ernst & Bar-Joseph, 2006), which applied for the clustering, comparing, and visualizing gene expression data from short time series experiments (in general, ∼8 time points or fewer). STEM implements a novel method for clustering short time series expression data that can differentiate between real and random patterns. The expression patterns were judged to be statistically significant between the number of genes expected (n_(E)) and the number of genes assigned (n) and potentially contained genes that were coordinately regulated.

To visualize the correlations between genes within the profiles, we constructed colored heat maps by plotting pair-wise correlation values of expression of all the genes within the profiles was performed with MultiExperiment Viewer (MeV) (Saeed et al., 2003). The average gene expression value of each significant cluster profile was correlated with each of nine obesity phenotypic traits using a non-parametric Spearman rank correlation coefficient with Bonferroni correction.

Gene functional enrichment analysis

A Gene Ontology (GO) and KEGG pathway enrichment analysis was used DAVID (Database for Annotation, Visualization and Integrated Discovery) web server (http://david.abcc.ncifcrf.gov/) (Huang, Sherman & Lempicki, 2008). The P values (i.e. EASE score), which indicated the significance of the comparison, was calculated by Benjamini-corrected modified Fisher’s exact test. Only GO and pathway categories with a P value less than 0.05 were considered as significant and listed.

DE genes in QTLs region

QTL data were downloaded from the Pig Quantitative Trait Locus database

(PigQTLdb: http://www.animalgenome.org/QTLdb/pig.html) website (Hu et al., 2013). PigQTLdb release 23 (April 21, 2014) contains 10,497 QTLs from 416 publications representing 647 different pig traits. Here, we defined QTL genes as those that have an overlapping region with QTL regions, and the overlapping region is at least half the length of the gene or the QTL region, whichever is shorter. In this study, ∼282.57 Mb QTL regions of the 2,311 genes were used for analysis. These were assembled from 901 high confidence and narrowed (<2 Mb) QTL affecting fatness and fat composition.

q-PCR validation

Total RNA were treated with RNase-free DNase I (TaKaRa, Katsushika, Tokyo, Japan). cDNA synthesis and q-PCR was performed using the SYBR® Prime- Script® RT-PCR Kit (TaKaRa) on a CFX96 Real-Time PCR detection system (Bio-Rad, Hercules, CA, USA). The PCR conditions were 5 min at 42 °C, 10 s at 95 °C, and then 40 cycles of 5 s at 95 °C and 30 s at 65 °C. The primers of 12 genes (TUBA1B, SCD, COX1, COX2, COX3, ATP8, CYTB, ND1, ND2, ND3, ND4, and ND5) designed for q-PCR analyses are listed in Table S1. q-PCR for each RNA sample was performed in triplicate. β actin (ACTB), TATA box Binding Protein (TBP) and Topoisomerase II β (TOP2B) genes were simultaneously used as an internal gene for normalization. The mRNA expression was quantified using the 2^−ΔΔCt method, data are expressed as mean ± SD.

Phenotypic measurements

As shown in Fig. 1, the body weight (Mann–Whitney U test, P < 10⁻⁶) and adipocyte volumes (Student’s t-test, P < 10⁻⁴) were significantly different among the four stages. Additionally, measurement of eight representative serum adipose metabolism indicators gave the same ranking (One-way ANOVA, P < 0.05, Fig. S1). These phenotypic differences at various stages of HLB imply the existence of intrinsic molecular differences.

Analysis of DGE profiling libraries

To investigate gene expression changes during development, 12 porcine HLB DGE libraries were constructed using Illumina DGE methods. These DGE libraries generated 3.66 to 6.5 million raw tags for each of the 12 libraries. After filtering, the total number of clean tags per library produced ranged from 3.32 to 6.04 million and the number of distinct clean tags ranged from 141,865 to 270,124 (Table S2). To estimate the quality of the DGE data, the saturation and distribution of clean tag expression was analyzed (Figs. S2–S4).

For tag mapping, one reference tag database that included 22,293 sequences from Ensembl Sscrofa 10.2 was preprocessed. We obtained 177,693 total reference tag sequences and 164,561 unambiguous tag sequences. Tolerances were set to allow one mismatch in each alignment to take into account polymorphisms across samples. Among the distinct tags, the number that could be mapped to genes ranged from 54,177 to 69,786 (Table S2). A total of 12,618 genes were identified, with at least one tag in all analyzed samples. The level of gene expression was further analyzed by calculating the number of unambiguous tags for each gene and then normalizing this to the number of TPM. The results showed that the majority of genes transcribed were represented by fewer than 10 copies and only a small proportion of genes were highly expressed (Fig. S5). This result suggests that many genes have housekeeping cellular roles (and therefore low expression) and these may be the main regulatory mRNAs in adipogenesis and basal cellular metabolism (Guo & Liao, 2000).

Specifically, a set of genes (ATP8: ATP synthase protein 8; COX2: Cytochrome c Oxidase 2; COX3: Cytochrome c Oxidase 3; ND1: NADH Dehydrogenase 1; ND2: NADH Dehydrogenase 2; SCD: Stearoyl-Coenzyme A Desaturase; TUBA1B: Tubulin α1b) shared by the four libraries showed over 100-fold increased expression compared with the average level (total expression of all genes/gene number) in each library, which suggests essential roles for these genes in the growth and development of adipose tissue. Stearoyl-coenzyme A desaturase, which catalyzes the de novo synthesis of monounsaturated from saturated fatty acids, contributes to adipocyte differentiation by synthesizing monounsaturated fatty acids that are incorporated into the triglycerides of the mature adipocytes (Dobrzyn & Ntambi, 2005). Expression of the SCD gene increased gradually during porcine adipocyte differentiation, indicating a significant role of SCD in this process (Wang et al., 2004a). In addition, increased SCD gene expression is associated with obesity (Jones et al., 1996) and elevated SCD gene activity in adipose tissue is closely related to insulin resistance (Sjögren et al., 2008). Furthermore, the SCD gene has served as a postnatal marker of adipocyte differentiation in bovine subcutaneous adipose tissue (Martin et al., 1999). The TUBA1B gene has been shown to play an important role in the metabolic actions of insulin, and a mutation in this gene has been associated with insulin resistance (Kapeller et al., 1995). Cytochrome c oxidase and NADH dehydrogenase showed increased gene expression levels during adipocyte differentiation and intramuscular fat deposition, suggesting that these genes play important roles in fat development (Lee et al., 2012). Deficiencies in COX1 and COX2 gene activity have been correlated with a number of human aging diseases such as Alzheimer’s and muscle atrophy (Barrientos et al., 2002). Obese mice also have higher COX activity levels than lean mice (Mercer & Trayhurn, 1987; Hebert, Ratnasingham & deWaard, 2003). In addition, the mitochondrial ATP8 mutation induces mitochondrial Reactive Oxygen Species (ROS) generation, secretory dysfunction, and β-cell mass adaptation in mice and has been associated with type-2 diabetes and obesity (Weiss et al., 2012).

The expression of these genes was validated to ensure the changes were real, and not caused by technical error. q-PCR was applied to investigate the relative expression levels of these genes, and five other genes, in each of the four developmental stages. We correlated the DGE and q-PCR results, which showed that the expression patterns of these genes were consistent with each other (Pearson’s r > 0.7, P < 0.05; Fig. S6). The fold-change in expression among the stages was often lower when assessed by q-PCR compared with DGE profiling; however, these differences may derive from the intrinsic differences between these two technology platforms (Linsen et al., 2009).

Analysis of DE genes

To explore global transcriptional changes, we identified 3,274 DE genes across the four developmental stages. Gene ontology (GO) and pathway analysis revealed that DE genes were significantly enriched in the oxidative stress, immune process, apoptosis, energy metabolism process, hormone stimulus, cell cycle, angiogenic and translation categories (Fig. S7), such as ‘response to oxidative stress’ (P = 1.74 × 10⁻⁴, 40 genes), ‘regulation of immune system process’ (P = 3.19 × 10⁻⁴, 52 genes), ‘leukocyte activation’ (P = 4.43 × 10⁻³, 48 genes), ‘regulation of apoptosis’ (P = 1.94 × 10⁻², 128 genes), ‘phosphorylation’ (P = 1.65 × 10⁻², 128 genes), ‘response to insulin stimulus’ (P = 1.09 × 10⁻², 23 genes), ‘blood vessel development’ (P = 2.21 × 10⁻², 45 genes) and ‘transcription factor binding’ (P = 4.47 × 10⁻², 83 genes). Oxidative stress accumulating over a lifetime can inflict direct cellular damage and influence various signaling pathways and transcriptional programs regulating key development processes including immunity, energy metabolism, and apoptosis (Finkel & Holbrook, 2000). This suggests that interactions between biological processes affect porcine adipose development and growth. The increased levels of ROS observed in aging adipose tissue are likely a reflection of the modulatory action of age-related levels of hypoxia on adipose function. Hypoxia increases insulin-stimulated glucose uptake and decreases the lipid content of differentiated 3T3-L1 cells (Zhang et al., 2011). In addition, aging-associated oxidative stress reduces the fat mass in mice, while in vitro studies have found that increased oxidative stress through glutathione depletion alters pre-adipocyte differentiation (Zhang et al., 2011). Adipose tissue growth comprises of the enlargement of existing adipocytes and the formation of new cells from committed pre-adipocytes. This is accompanied by the coordinated development of a well-defined vascular system, which ensures that every adipocyte is surrounded by one or more capillaries. Angiogenesis plays an important role in adipose development not only by supplying nutrients and oxygen to nourish adipocytes, but also as a cellular reservoir of adipose precursor and stem cells that control adipose tissue mass and function (Cao, 2013). As part of the immune system, adipose tissue has also been identified to be involved in aging-associated changes in immune and inflammatory responses. Aged visceral fat in mice showed higher levels of the pro-inflammatory cytokines IL-1β, IL-6, TNF-α and COX-2, but lower expression of the anti-inflammatory PPAR-γ, chemokine and chemokine receptor (Wu et al., 2007; Lumeng et al., 2011). Starr, Evers & Saito (2009) revealed that IL-6 increased in white adipose tissue with age (Starr, Evers & Saito, 2009). The particular role of adipose tissue in calorie storage makes adipocytes well suited to the regulation of energy balance. Some adipocyte-secreted proteins, including leptin, adiponectin and resistin, play important roles in regulating whole-body metabolism (Galic, Oakhill & Steinberg, 2010; Kershaw & Flier, 2004). To meet the energy requirements of fat and other tissue development, adipocytes allocate energy by responding to hormones, cytokines, and other factors that are involved in energy metabolism (Frühbeck et al., 2001; Klaus, 2004).

To evaluate the reproducibility of DGE library sequencing, we analyzed the expression variations of the 3,274 DE genes by hierarchical clustering analysis to indicate the reliability of our experimental results as well as operational stability. As shown in Fig. 2A, the three biological replicates were highly correlated with each other (average Pearson’s r = 0.972, P < 10⁻¹⁶, Fig. S8) and all individuals could be clearly assigned to a group, suggesting experimental reliability and further highlighting the low variation in HLB mRNA profiles across different individuals. Interesting, the mRNA expression profiles show an obvious development-specific pattern. Two major branches are defined: one representing birth, 30 d and 180 d, and one representing 7 y. This different clustering pattern may correspond to the intrinsically different biochemical and physiological properties of porcine adipose tissue. Principal components analysis clearly recapitulated these findings. Furthermore, within the treatment groups, gene expression varied significantly by development and time point, as depicted by the well-defined clustering of treatment groups (Fig. 2B). Our previous comparative study of methylation and gene expression profiles of 180 d and 7 y porcine skeletal muscle indicated that in the aged (7 y) pig, symptoms of muscle atrophy have emerged, suggesting that porcine adipose tissue at this point may appear distinct or show specific gene expression patterns compared with earlier stages of growth and development (Jin et al., 2014); this also supports the findings of our current study.

We also examined the chromosomal distributions of the 3,274 DE genes, but did not find over- or under-representation of highly DE genes within a particular chromosome, except for the mitochondrion (P < 10⁻⁴, χ²-test, Table 1). Through BLAST analysis of the 3,274 DE gene sets against high confidence and narrowed (<2 Mb) QTLs affecting fatness and fat composition in the PigQTLdb (Hu et al., 2013), 327 (14.62%) genes were shown to overlap with the defined DE genes (Table 1). This highlights the potential of identifying candidate genes that may be involved in lipogenesis, development and growth. One such gene is the pig Complement Factor B (CFB) gene, which was localized to chromosome 7 within the interval of QTLs for backfat thickness, backfat weight, fat area and percentage of backfat in the carcass (Bidanel et al., 2001; Malek et al., 2001; Milan et al., 2002; Gilbert et al., 2007). Another gene is Activating Transcription Factor 6 (ATF6), which was located in the QTL region for body weight, backfat thickness, backfat weight, and linoleic acid percentage (Horák et al., 2010; Fontanesi et al., 2012). Meanwhile, obesity is associated with induction of the endoplasmic reticulum stress response signaling and that ATF6 is the important genes/proteins involved in the signaling regulation (Agouni et al., 2011), which supports the notion that pig might be a good model for further adipose studies.

DE genes involved in inflammatory and lipogenesis process

Numerous evidence indicates that inflammation is a possible underlying basis for the molecular alterations that link aging and age-related pathological processes (Chung et al., 2006). Based on the annotation of the Pathway Central database (SABiosciences, MD, USA), we identified 38 orthologous inflammation-related genes between pig and human, of which only five genes (5/38, 13.15%) were DE across porcine distinct developmental stages, indicating that there was no obvious change in porcine HLB tissue with the aspect of inflammation till middle age (7 y). Nonetheless, the expression level of two inflammatory genes, CRP and IL13RA1, which plays an important role in anti-inflammatory (Marnell, Mold & Du Clos, 2005; Alfarsi, Hamlet & Ivanovski, 2014), have been increasing along with development (Fig. S9A), which may serve as early biomarkers of aging-related inflammation. In addition, the expression of other three inflammatory genes, CX3CR1, CXCL9 and IL10RA, showed a steady increase during birth to 180 d followed the sharp decrease in 7 y (Fig. S9B), which also need to pay attention in further aging related studies.

Furthermore, previous studies in human have identified several key candidate genes involved in lipid metabolism, such as FASN, PPARG, CEBPA, SLC1A4, ADIPOR1 and ADIPOR2 genes (Wang et al., 2004b; Kersten, Desvergne & Wahli, 2000; Yamamoto et al., 2010). Surprisingly, these genes in our study were not DE across porcine distinct developmental stages in HLB tissue; however, they also showed a trend consistent of their biological function. Lipogenesis genes were high expression in 30 d and 180 d (Fig. S9C), in contrast, lipolysis genes were low expression in 30 d and 180 d (Fig. S9D). Of course, we also found DE genes related lipogenesis, such as CEBPG, CEBPZ and PRKAR1A genes that had higher expression level in Birth than other stages (Fig. S9E), and LIPA, HEXB and ACE genes in 180 d had highest expression level across four stages (Fig. S9F), suggesting that the governance of lipogenesis in distinct developmental stage are different genes.

Analysis of STEM cluster

To explore temporal gene expression patterns, we used the STEM algorithm to search the most probable set of clusters generating the observed time series. The 3,274 DE genes were clustered into 24 distinct expression patterns, of which 1,097 DE genes were significantly over-enriched in six expression patterns (Fig. 3A), which have significantly more genes assigned under the true ordering of four stages than the average number assigned to the model profile in the permutation runs. To investigate whether the change of significantly clustered genes affected phenotype variation, we averaged the expression of each significant cluster profile and performed association analysis with the phenotypic traits (Fig. 3B). We identified that these significant cluster profiles significantly correlated with the amount of Lip-a (Spearman r = 0.81, P = 1.61 × 10⁻³), VLDL (Spearman r = 0.78, P = 2.59 × 10⁻³; Spearman r = −0.82, P = 1.14 × 10⁻³), TG (Spearman r = 0.83, P = 7.85 × 10⁻⁴; Spearman r = −0.85, P = 5.21 × 10⁻⁴), LDL (Spearman r = 0.79, P = 1.9 × 10⁻³) and TC (Spearman r = −0.79, P = 2.33 × 10⁻³) in serum.

As shown in Fig. 3C, we found that profiles 1, 2, 3 and 4 are highly positive correlated to each other, as are profiles 5 and 6. We therefore combined these genes for further functional enrichment analysis. Significant GO categories of up-regulated genes (profile 5 and 6) correlated with phosphorylation and angiogenesis (Fig. 3D). Various well-known genes involved in phosphorylation and angiogenesis of adipose were identified. For example, Pyruvate Dehydrogenase Kinase 2 (PDK2) and Pyruvate Carboxylase (PC) are involved in phosphorylation, which continually increased with adipose growth and development (Fig. S10), and may play a role in diabetes, lipogenesis and insulin secretion (Lee, 2014; Xu et al., 2008). Angiopoietin 2 (ANGPT2) and Epidermal Growth Factor-Like domain 7 (EGFL7) are increased with adipose growth and development in general (Fig. S10), and could be regulating vasculogenesis (Wang, Huebert & Shah, 2014; Bambino et al., 2014). These results are consistent with these genes and their associated processes having important roles in the process of adipose development. Significant GO categories of down-regulated genes (profiles 1, 2, 3 and 4) correlated with cell cycle and cytoskeleton organization (Fig. 3D), consistent with the knowledge that inhibition of the cell cycle occurs during the later stages of growth and development. Genes involved in cell cycle were identified. For example, minichromosome maintenance complex component 2 and 3 (MCM2 and MCM3), which are involved in the initiation of eukaryotic genome replication (Kang, Warner & Bell, 2014); Cyclin D3 (CCND3), which functions as a regulator of CDK kinases, whose activity is required for cell cycle G1/S transition (Li et al., 2013); F-Box Protein 5 (FBXO5), which is a mitotic regulator that interacts with CDC20 and inhibits the anaphase promoting complex (Reimann et al., 2001); both of those genes are decreased with adipose growth and development (Fig. S10).

DE genes in middle-aged pig

The results of STEM analysis revealed that gene expression patterns are up- or down-regulated in general; however, the gene expression pattern was changed between 180 d and 7 y in profile 3, suggesting that the growth and development of 7 y porcine adipose tissue represents a distinct change. Previous studies mainly carried out in pigs are in neonatal or before they reach the age of 1 year (Zhang et al., 2013; Zhou et al., 2013). Limited studies have been carried out using young adult pigs (2 years old). Therefore, we used middle-aged pigs (7 years old) as representative of slightly older individuals. A comparative study of methylation and gene expression profiles of 180 d and 7 y porcine skeletal muscles indicated that 7-year-old pigs have emerging symptoms of aging (Jin et al., 2014). We therefore investigated whether known age-related genes from the Human Ageing Genomic Resources (HAGR) database (Tacutu et al., 2013) were included in the genes found to be DE (7 y vs. birth, 7 y vs. 30 d, 7 y vs. 180 d, P < 0.01, Student’s t-test) in our data set. Within the HAGR, GenAge is a database of 288 genes potentially associated with human age, 61 of which were found to be DE with age in our study (P = 0.017, χ²-test, Table S3). GO analysis revealed that these 61 genes are associated with DNA repair, apoptosis, transcription and immune processes. Various well-known genes were identified, such as Insulin-like Growth Factor 1 (IGF1), and Peroxisome Proliferator Activated Receptor α (PPARα), a Mechanistic Target of Rapamycin (MTOR) and Uncoupling Protein 2 (UCP2).

Highly expressed DE genes were related to mitochondrial energy metabolism

Surprisingly, five of the highly expressed DE genes were mitochondrial genes. It is well known that energy metabolism is closely related to mitochondrial gene expression (Patti & Corvera, 2010). To verify that the mitochondrial genes were not in different stages of DE by chance, we performed a χ²-test and found that out of 13 mitochondrial genes, nine (69.23% P < 0.05) exhibited dynamic expression changes during the postnatal development stages in HLB (Table 1). Moreover, a previous study reported that the mtDNA copy number was positively correlated with energy metabolism (Tang et al., 2000). We therefore used the mtDNA copy number of HLB as a means of partially reflecting the energy metabolism. The mtDNA copy number of HLB was significantly increased during the process of development (One-way ANOVA, P = 9.24 × 10⁻⁴), which was consistent with the expression of mitochondrial genes in our study (Pearson’s r = 0.89, P = 1.15 × 10⁻⁴, Fig. 4A).

We also quantified the expression level of six mitochondrial genes in various adipose tissues between 180 d and 7 y. Expression of all six genes was higher in SATs than in VATs (P < 0.05, Student’s t-test, Fig. 4B), which was comparable with the relative levels of mtDNA copy number between the two types of adipose tissue (P < 0.01, Student’s t-test, Fig. 4C). These results revealed distinct functional properties of SATs and VATs, in which SATs mainly store energy but VATs are involved with inflammatory and immune processes (Ibrahim, 2010). A comprehensive survey of mtDNA copy number across different tissues including brain and muscle was performed. Of all the tissues examined, muscle showed the highest mtDNA copy number (Fig. 4D), which is consistent with the higher energy demand of muscle due to its motor function. The mtDNA copy number of brain was relatively lower than that of muscle, but higher than all the visceral organs and adipose tissues (P < 0.05, Student’s t-test, Fig. 4D), suggesting that muscle and brain have the highest energy requirements (Belanger, Allaman & Magistretti, 2011; Magistretti & Allaman, 2013).