t i l e Differentially Expressed Genes Analysis of Brown Adipose Tissue During Cold Acclimation in Male Tree Shrews ( Tupaia belangeri ) based on RNA-Seq

pathway Thermogenic function of brown adipose tissue (BAT) was known to be markedly elevated when animals were exposed to the cold. In the present study, transcriptome sequencing of BAT in Tupaia belangeri between control and cold acclimation group were carried out by Illumina novaseq 6000 platform. Then the pathways and key genes related to lipid metabolism and energy metabolism were screened out by bioinformatics analysis. The results showed that there were significant differences in the expression of 2879 genes in cold acclimation group compared with that of control group, 1181 were up-regulated and 1698 were down-regulated in cold acclimation group. Pathway analyses showed that the differentially expressed genes were significantly enriched in oxidative phosphorylation, steroid biosynthesis, glycerophospholipid and fatty acid metabolism pathways related to energy metabolism and lipid metabolism. Differentially expressed genes ND3 , ND4 , ND5 , ND6 , ND4L , CYTB , ATP8 and ATP6 in oxidative phosphorylation pathway were significantly up-regulated. Key genes related to lipid metabolism were identified as ELOVL6 , ACACA , HSD17B12 , SQLE , LSS , DHCR7 , DHCR24 , CYP51A1 , HMGCS1 , LPL , and ACACA , etc, which were down-regulated in the lipid metabolism pathway. All of the above results indicated that T. belangeri could adapt to the cold environment by regulating genes expression in oxidative phosphorylation, steroid biosynthesis, glycerophospholipid and fatty acid metabolism pathways in BAT. Among them, expressions of genes in oxidative phosphorylation pathway were up-regulated, and most of the genes related to the lipid metabolism pathway were down-regulated, which promoted energy expenditure and thermogenesis.


INTRODUCTION
C old exposure imposes a metabolic challenge to mammals that must be met by the response of animals' body tissues to ensure homeothermy. Brown adipose tissue (BAT) was the main site of cold-induced non-shivering thermogenesis (NST) (Klingenspor, 2003), which has the remarkable ability to dissipate excess energy as heat in a process known as adaptive thermogenesis. Transcriptome sequencing technology, as the main method of studying transcriptomics, has high sensitivity in transcript identification and the quantification of gene expression (Conesa et al., 2016). Researchers used transcriptome sequencing technology to analysis the differences of BAT and white adipose tissue (WAT) in mice under cold exposure, the results indicated that the effects of cold exposure on the gene expression profiles of BAT were much greater than those of WAT (Watanabe et al., 2011). Other study have shown that expression of genes related to lipid metabolism and oxidative phosphorylation pathways in BAT of mice or rats significantly changed under cold exposure (Shore et al., 2013;Rosell et al., 2014).
Tree shrew, Tupaia belangeri (Mammalia: Scandentia: Tupaiidae), a squirrel-liked lower primate, is a unique Oriental species. This animal originates from the tropical O n l i n e

F i r s t A r t i c l e
island, in China; and is mainly distributed in Yunnan, Guizhou and Sichuan (Bremer et al., 2011). Judging from the distribution characteristics of existing tree shrews, the Yunnan-Kweichow Plateau and the southeastern part of the Qinghai-Tibet Plateau, that is, parts of the Hengduan Mountains may constitute the northern limit of the distribution of T. belangeri, and low temperature may limit its northward spread. Our group previously used classic physiological research methods to confirm that cold acclimation increased the BAT mass, uncoupling protein 1 (UCP1) content, and NST in T. belangeri to maintain body temperature, but the proportion of NST in the total heat production showed a gradual decline (Zhang et al., 2017). In the present study, we used RNA-seq to analyze the differential expression of genes in BAT of T. belangeri under cold acclimation, providing relevant information such as genes function and enrichment, and realizing the analysis of the function of thermogenesis genes and pathways, thus providing important reference for the study of the mechanism of cold adaptation heat-producing in small mammals.

Animals
T. belangeri were captured (25°25′-26°22′N, 102°13′-102°57′E, at 1679 m altitude) at the boscage of Luquan County, Yunnan Province, and maintained at the School of Life Sciences, Yunnan Normal University, Kunming (1910 m altitude). Average body mass of all animals was 115.36±3.65g, after one month of adaptation at room temperature; the animals were divided into 2 groups according to their body mass, with 6 male individuals in each group: One group was kept at 25℃±1℃ for 28 days; these animals comprised the control group. A second group was kept at 5℃±1℃ for 28 days comprised the cold acclimatized group. Two groups were healthy adults, and housed individually in a wire cage (40 cm×40 cm×40 cm), The cage environment was maintained at 12L:12D (lights on at 08:00), and 65-92% relative humidity. Water and foods were provided ad libitum. T. belangeri were fed a food mixture containing 30% corn flour, 20% wheat flour, 5% fish meal, 6% wheat bran, 3.6% milk powder, 10% sugar, 2% yeast, 3% electrolytic multivitamin, 0.4% edible salt, and 20% eggs, as well as apples, pears and other fruits twice weekly. After 28 days of acclimatization and testing, all animals were anesthetized with ether and then killed. Samples of BAT was taken from each animal and immediately frozen in liquid nitrogen followed by storage at -80℃ prior to analysis.

Total RNA extraction and RNA-Seq preparation
Total RNA was extracted from BAT according to miRNeasy Mini Kit instructions. The integrity, purity and concentration of total RNA were confirmed using a Bioanalyzer 2100 (Agilent Technologies). The total RNA content > 3μg, and the integrity (RIN) > 7.0 was qualified, which can be used in the next experiment. cDNA synthesis was done using TruSeq PE Cluster Kit v4-cBot-HS (Illumia). RNA-seq libraries were constructed according to the manufacturer's instructions and sequenced using the Illumina novaseq 6000 platform.

RNA-Seq data analysis
Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC content and clean data were calculated. All the downstream analyses were based on clean data with high quality. These clean reads were then mapped to the tree shrew reference genome sequence (version tupBel1). Hisat2 tools soft were used to map with reference genome. Quantification of gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped. The formula is shown as follow: FPKM= cDNA Fragments/ mapped fragments (millions) × transcript length (kb) Through PCA principal component analysis of two groups of data under different conditions, the difference and distance between groups were reflected on the coordinate axis. Differential expression analyses of two groups were performed using the DEseq. The screening criteria were fold change (FC) ≥ 2 and false discovery rate (FDR) < 0.01 (Shore et al., 2013), and count of differentially expressed genes (DEGs) of log 2 FC >4. Functional enrichment analysis was performed with KOBAS using pathways related to metabolism from the KEGG database annotation (Mao et al., 2005). The significant enrichment Q-value (Q-value is the P-value corrected for multiple hypothesis testing) was calculated to identify the pathways with significant enrichment, and key genes were screened based on FPKM values and related literature reports.

Quantitative real time PCR
To verify the reliability of transcriptome sequencing results, 8 DEGs were selected for RT-PCR quantification (Table I), and the reference gene primer sequence was F: GAGAGGGAAATCGTGCGTGAC R: CATCTGCTGGAAGGTGGACA The design of DEGs primer and the steps of real-time quantitative operation are detailed Zhang et al. (2017).

Transcriptome data evaluation
All samples were qualified (total RNA content> 3 μg, RIN≥7). The proportion of clean data of all samples with Q30 was more than 93.62%, and the GC content were about 50%, the quality of the data obtained by sequencing was reliable (Table II). The results of principal component analysis showed that the data composition of the two groups was significantly different, the samples in the groups were clustered, and the obtained results were consistent with the expected results ( Fig. 1).

Bioinformatics analysis
Using a significance level of FDR < 0.05, FC ≥ 2, we found a total of 2,879 DEGs, of which 1,181 and 1,698 genes showed increased and reduced expression under 28-day cold exposure, respectively, Notably, the UCP1 (log 2 FC=1.316063) gene encoding UCP1 reported previously as having increased expression in BAT after cold exposure. The present study (Tables III) showed that the expression of 28 genes were significantly up-or down-regulated more than four folds by cold exposure, these include ATP8, ATP6, ND3, ND4L, ND4, ND5 and CYTB genes in oxidative phosphorylation pathway, as well as ACACA, ELOVL6 and HSD17B12 genes that played an inhibitory role in fatty acid metabolism, and so on. Functional enrichment analyses of the genes with the most pronounced cold-induced expression using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways revealed a significant enrichment of metabolic pathways, including oxidative phosphorylation, steroid biosynthesis, fatty acid metabolism, glycerophospholipid metabolism (Fig. 2). The most highly enriched pathway for up-regulated genes were oxidative phosphorylation pathway, which includes several genes encoding NADH dehydrogenase (ND3, ND4, ND5, ND6 and ND4L), cytochrome c reductase (CYTB) and ATP synthase (ATP6, ATP8). Most genes related to lipid metabolism pathway were repressed in BAT, including SQLE, LSS, DHCR24, DHCR7, CYP51A1, SCD, ACSL4, ELOVL6, ACACA, HSD17B12, GPD2, LPGAT1, and SELENOI (Table IV).

O n l i n e F i r s t A r t i c l e
Genes Analysis of BAT During Cold Acclimation in Tupaia belangeri  Figure  3. The up and down regulation trends of gene expression between the two results were consistent, indicating that the quality of transcriptome sequencing was reliable.

DISCUSSION
Small mammals in the wild maintain body temperature by increasing thermogenesis when they are exposed to low temperatures (Wang and Wang, 1990). Studies showed that there were significant differences in gene expression in BAT of yak between cold and warm seasons (Zhu et al., 2021). In the present study, it was found that 2879 genes were up-regulated or down-regulated in the cold acclimation group, suggesting that T. belangeri may adapt to the cold environment by changing the expression pattern of genes. Studies suggested that UCP1 was the main substance that determines the heat production of BAT (Mei et al., 2019). In our study, expression level of heat-producing gene UCP1 increased after cold acclimation, which was consistent with the results of previous studies, but the up-regulated multiple was not significant (log 2 FC =1.3), indicating the reason why the proportion of NST in total heat-producing decreased gradually (Zhang et al., 2017). In the KEGG enrichment analysis of this study, DEGs were significantly enriched in oxidative phosphorylation related to energy metabolism, steroid biosynthesis, fatty acid metabolism, and glycerolipid metabolism pathways related to lipid metabolism, suggesting that lipid metabolism pathways in T. belangeri may play an important role in the mechanism of cold adaptation. Previous analysis of transcriptome pathways showed that the genes most significantly induced in mouse BAT under cold exposure were those involved in glycerophospholipid synthesis and fatty acid elongation (Rosell et al., 2014). Mitochondrial protein-coding genes played a major role in the oxidative phosphorylation system (Leonard and Schapira, 2000), and the four enzymes encoded by them were the driving force for ATP synthesis. Studies showed that after long-term cold exposure in mice, gene expression in BAT showed significant changes, and the oxidative phosphorylation pathway was significantly up-regulated in BAT (Rosell et al., 2014). In present study, compared with the control group, mitochondrial protein-coding genes NADH Dehydrogenase subunit 3(ND3), NADH dehydrogenase subunit 4 (ND4), NADH dehydrogenase subunit 5 (ND5), NADH dehydrogenase subunit 6 (ND6),

D-M. Hou et al.
NADH dehydrogenase subunit 4L (ND4L), cytochrome b (CYTB), ATP synthase F0 subunit 6 (ATP6) and ATP synthase F0 subunit 8 (ATP8) were highly expressed in the oxidative phosphorylation pathway of BAT in the cold acclimatization group, all of which were unregulated by more than 4 fold. These results indicated that the synthesis of ATP was promoted in the response to cold environment, which provided more energy for animals. Studies have shown that mitochondrial protein-encoding genes defects could lead to increased mitochondrial damage and decreased function, which ultimately leaded to a decrease in the ability of oxidative phosphorylation (Luoma et al., 2015). NADH dehydrogenase encoded by ND3, ND4, ND5, ND6 and ND4L was one of the largest membrane-bound enzymes in cells. Studies have shown that the variation of ND3 could cause a loss of the role of mitochondrial oxidative phosphorylation. If the mitochondrial oxidative energy supply was increased, Might cause the expression level of ND3 to also increase (Hinokio et al., 1995;Wang et al., 2013). In this study, ND3 was up-regulated 6 folds, suggesting that T. belangeri under cold acclimationthe may promote mitochondrial oxidative energy supply through the regulation of this gene. Inhibition of protein synthesis, transcription, and ion transport can lead to inhibition of energy expenditure. Among the 13 protein-coding genes, the structure and function of CYTB gene has been most thoroughly studied, it was an important mediator of electron transfer in oxidative phosphorylation, and could conduct transcription and protein synthesis in mitochondria itself (Schoepp et al., 2000). Therefore, our data demonstrate that T. belangeri might adapt to the cold environment by the regulation of this gene to maximize energy expenditure. ATP6/8 of the genes encoding ATP synthase protein 6/8. Studies have shown that the expression of ATP6/8 in BAT during hibernation in black-line hamsters into also up-regulated by 4 folds. If the expression of ATP6 gene was down-regulated, ATP production can reduce, thereby affecting energy metabolism (Hittel et al., 2013). Lipid was the main component of biofilm and important structure of organism. It not only participates in energy supply and storage of organism, but also can be transformed into various derivatives to participate in metabolic activities. Studies have performed transcriptome analysis of BAT in mice under cold exposure; they proposed the response to cold stress involves decreased gene expression in a range of lipid metabolism in order to maximize pathways involved in heat production (Shore et al., 2013). Our data demonstrate that the general trend was for a reduction in the expression of a large number of genes involved in steroid biosynthesis, glycerophospholipid and fatty acid metabolism pathways related to lipid metabolism, and far fewer genes were up-regulated by the cold challenge. This might represent a mechanism for reducing metabolic reactions so that metabolism could be focused on the cell thermogenic functions, particularly in BAT, that are essential to prevent hypothermia. Steroids were mainly composed of cholesterol. Studies have shown that insufficient cholesterol content could cause more free fatty acid (FFA) to remain in the circulatory system of mice to provide energy sources for other tissues and organs of the animal body (Dugail et al., 2003). In our study, compared with the control group, the squalene epoxidase (SQLE), lanosterol synthase (LSS), 24-dehydrocholesterol reductase (DHCR24), 7-dehydrocholesterol reductase (DHCR7), and lanosterol synthase (LSS) genes encoding proteins related to cholesterol synthesis were significantly down-regulated in the cold acclimation group. These results suggested that T. belangeri might reduce the synthesis of cholesterol in BAT by down-regulating SQLE, LSS, DHCR24, DHCR7, and CYP51A1 in cold environments, which leaved more FFA in the circulatory system to provide energy source for other tissues and organs of animal body. SQLE and LSS were located in the endoplasmic reticulum, and squalene epoxidase encoded by SQLE was one of the key ratelimiting enzymes that catalyze the conversion of squalene to squalene 2, 3-epoxides in cholesterol biosynthesis. Studies have shown that inhibition of SQLE could effectively reduce cholesterol synthesis (Ferdinandusse et al., 2017). LSS was a member of the gene super family OSC (3S-2,3 epoxy-squalene cyclase), and the production of lanosterol catalyzed by LSS was a key step in the synthesis of sterols in animals. CYP51 was one of the most O n l i n e

F i r s t A r t i c l e
Genes Analysis of BAT During Cold Acclimation in Tupaia belangeri widely distributed members of the CYP superfamily of cytochromases. The CYP51A1 gene encoding this enzyme originated from prokaryotes and was widely distributed in eukaryotes (Lepesheva and Waterman, 2011). Glycerophospholipid were the most abundant phospholipids in animals. The metabolism of each phospholipid formed a complex network and was regulated by a variety of metabolic enzymes (Liu and Huang, 2013). Fatty acid metabolism mainly includes de novo synthesis, oxidation, desaturation and lengthening of fatty acids (FA) to produce FA with different saturation degrees and different carbon chain lengths. In our study, compared with a control group, stearoyl-CoA desaturase (SCD), acyl-CoA synthetase long-chain family member 4 (ACSL4), ELOVL fatty acid elongase 6 (ELOVL6), acetyl-CoA carboxylase alpha (ACACA) and hydroxysteroid 17-beta dehydrogenase 1 (HSD17B12) genes related to fatty acid metabolism and glycerol-3-phosphate dehydrogenase 2 (GPD2), lysophosphatidylglycerol acyltransferase 1 (LPGAT1) and selenoprotein I (SELENOI) involved in the glycerolphospholipid metabolism pathway were significantly down-regulated in cold acclimation group. This result was consistent with previous studies that showed significant changes in genes involved in glycerolphospholipid synthesis and fatty acid extension in BAT of cold-stimulated mice (Marcher et al., 2015). Among them, ACACA was a target gene of steriod response element binding protein-1 (SREBP-1), and the enzyme encoded by it was a key rate-limiting enzyme regulating de novo synthesis of fatty acids. Its expression was also regulated by a variety of hormones. Studies have shown that the down-regulation of the ACACA gene expression leaded to a decrease in fatty acid de novo synthesis (Knobloch et al., 2013;Bakhtiarizadeh and Alamouti, 2020). ACSL4 belonged to the long-chain family of acyl-CoA synthase, which mainly catalyzed the activation of long-chain fatty acids to synthesize cell lipids. Down-regulation of ACSL4 reduced the intracellular triglyceride concentration (Fan et al., 2020). ELOVL6 belonged to the sixth member of the ultralong-chain elongase gene family. It was an important rate-limiting enzyme for the extension and synthesis of long-chain fatty acids. Its main function in fatty acid metabolism was the extension of fatty acids above C16. Studies have shown that ELOVL6 gene down regulation of expression inhibited the extension of C16:1 to C18:1 (Corominas et al., 2013;Green et al., 2010). Fatty acid desaturation index in fat cells of ELOVL6 knockout mice decreased significantly, and the content of linoleic acid also decreased significantly (Sunaga et al., 2013). 17β-SHD12 belonged to the short-chain dehydrogenase/ reductase (SDR) superfamily and has been discovered relatively late in humans (Peltoketo et al., 1999), recently, it has been found that HSD17B12 gene regulates fatty acid synthesis and the elongation of long-chain fatty acyl-CoA in fatty acid metabolism (Moon and Horton, 2003). SCD was a key enzyme that catalyzes the desaturation of saturated fatty acids. Studies have shown that over expression of SCD can easily lead to fat accumulation (Nagao et al., 2019). The above results in present study suggested that T. belangeri may adapt to cold environment by down-regulating SCD, ACSL4, ELOVL6, ACACA, and HSD17B12 genes to regulate fatty acid synthesis in BAT.
The GPD2 gene encodes glycerol 3-phosphate dehydrogenase 2, which was the key enzyme linking glycolysis and oxidative phosphorylation, and its play an important role was as part of the glycerophosphate shuttle (Orr et al., 2014). Studies have shown that the decreased activity of GPD2 limited the gluconeogenesis of lactic acid and glycerol, and reduced endogenous glucose production, and could also inhibit the proliferation of tumor cells (Madiraju et al., 2014). In our study, GPD2 gene was significantly down-regulated, which could not only regulate energy balance, but also provided a target for disease treatment. As a remodeling enzyme, LPGAT1 could reconstitute lysophosphatidylglycerol into PG. Studies have found that the variation of LPGAT1 gene played an important role in the process of weight regulation (Traurig et al., 2013). Selenoprotein I (SELENOI) was also known as ethanolaminephosphotransferase 1 (EPT1). Studies have shown that mutations in the SELENOI gene significantly reduced the activity of encoding EPT1, leading to impeded PE synthesis and increased ATP synthesis, which interfered with the metabolic pathway (Ahmed et al., 2017;Ma et al., 2021). In this study, SELENOI gene was down-regulated, indicating that PE synthesis in BAT of T. belangeri was reduced under cold acclimation, thus affecting energy metabolism of the animals.
We measured the adaptive changes of genes in BAT in T. belangeri under cold acclimation at the transcriptome level. In summary, cold acclimation induced a significant up-regulation of genes expression related to oxidative phosphorylation, and significant down-regulation of genes expression related to steroid biosynthesis, glycerophospholipid metabolism, and fatty acid metabolism, which ultimately promoted energy expenditure and heat production. In addition, we also found that the heat-producing gene UCP1 was up-regulated after cold acclimation, but it was not significantly up-regulated, which might be the reason why the proportion of NST to total heat-producing decreased gradually in the process of cold acclimation and limited its northward spread the reason.