Gene expression profile analysis of pig muscle in response to cold stress

ABSTRACT

Low-temperature injury is one of the most significant causes of livestock damage worldwide. To identify genes that respond to cold stress, the transcriptome profile of the skeletal muscle of the Min pig (Sus scrofa) after cold stress was analysed using the deep sequencing technique. One group was raised in the normal environment (18 ± 2°C, n = 3) as the thermos-neutral group; another group was raised in the low-temperature environment (4 ± 2°C, n = 3) as the cold stress group. After 24 h, the gene expression profile of skeletal muscle was analysed. The results showed that the expression levels of 662 genes were significantly up-regulated, and those of 1058 genes were significantly down-regulated in the thermos-neutral and the cold stress groups (P < 0.01). It was found through Gene Ontology analysis that the significantly differently expressed genes were mainly distributed in the cytoplasm, mitochondrion, organelle envelope and ribosomal subunit. They enriched ribosome, fatty acid metabolism, RNA transport and metabolic pathways. It means that cold stress could significantly change the gene expression profile of pig muscle.

KEYWORDS:

• Pig
• cold stress
• deep sequencing
• gene expression profile

1. Introduction

Temperature is one of the primary environmental factors that influences and limits the growth and development of livestock and poultry species. Low and freezing temperatures affect animals not only in feed intake, but also in mortality (Mendes et al. 1997). Pigs are rather cold susceptible and are usually kept in heated housing when raised in cold regions (Young 1981). This is mainly because pigs are unusual mammals in a number of ways: they are the only ungulates that build nests for birth, and have large litters, where each young is behaviourally well developed but with a poor thermoregulatory ability. It is also striking that piglets appear to lack brown adipose tissue and rely on shivering as the main mechanism for thermoregulation (Berg et al. 2006). So, the piglet is more susceptible to cold stress.

When cold stress occurs, the cold stressor can regulate the synthesis and release of glucocorticoids (GCs) through the hypothalamic–pituitary–adrenal (HPA) axis (Girotti et al. 2011). Neurones in the hypothalamic paraventricular nucleus release corticotrophin-releasing hormone and other secretagogues into the pituitary portal system to stimulate secretion of adrenocorticotrophic hormone from the anterior pituitary, which in turn elicits release of adrenal corticosteroids (Ma & Morilak 2005). The release of GCs influences immune cells by suppressing differentiation and proliferation, regulating gene transcription and reducing expression of cell adhesion molecules that are involved in immune cell trafficking (Guo & DiPietro 2010). In our previous research, we found that cold stress down-regulates the expression of some immune gene (Zhang et al. 2012).

Not only is HPA the main regulative system, but muscle, liver and brain are also main heat production organs. Because the volume of the liver and brain are much smaller than the muscle, the muscle is the body's most important heat production site, especially skeletal muscles. Skeletal muscles serve a variety of functions, including movement and posture maintenance, heat production, protection and pressure alteration to aid circulation. More than 75% of the energy used during a muscular contraction is released as heat, and muscles considered as the body's furnace (Peter 2005). Until now, there was almost no reports about molecular mechanism of pig muscle under cold stress, that limited the research on protect pig avoid low temperature injury, even, on cultivate a cold tolerance pig breed.

2. Materials and methods

2.1. Animals

Six 75-day-old Min piglets from one nest were randomly divided into two groups: the thermos-neutral group (C, n = 3) and the cold stress group (E, n = 3). The piglets of the thermos-neutral group were placed in a room which was heated to 18 ± 2°C and the piglets of the cold stress group were placed in a house which was cooled to 4 ± 2°C. The relative humidity was 40–45% in both rooms. Both groups were fed the same feed. After 24 h, piglets were killed and the skeletal muscle from different piglets were collected and stored at −80°C until RNA isolation. All experiments were approved by the Management Regulation for Laboratory Animals of China.

2.2. High throughout sequencing

Total RNA was extracted from muscle tissue and posted to BGI Company (Shenzhen). The Illumina Solexa platform was used to obtain the raw sequences. Raw sequences were transformed into clean Tags by removal of 3′ adaptor sequences, empty reads, low-quality Tags, Tags that were too long or too short and Tags with a copy number of 1.

2.3. Tag annotation and data normalization for gene expression level

All clean Tags were mapped to the reference sequences (Sscrofa10.2) and were allowed no more than 1 nt mismatch. Clean Tags mapped to reference sequences from multiple genes were filtered. The remaining clean Tags were designated as unambiguous clean Tags. The number of unambiguous clean Tags for each gene was calculated and then normalized to TPM (number of transcripts per million clean Tags).

2.4. Statistical evaluation

Denote the number of unambiguous clean Tags from gene A as x; as every gene's expression occupies only a small part of the library, the p(x) is in the Poisson distribution.

The total clean Tag number of the C group is N1, and total clean Tag number of E group is N2; gene A holds x tags in the C group and y tags in the E group. The probability of gene A expressed equally between the two groups can be calculated with:

False discovery rate (FDR) ≤ 0.001 and the absolute value of log2 Ratio ≤ 1 as the threshold for significance of gene expression difference.

2.5. Functional annotation

All differentially expressed genes (DEGs) were mapped to terms in the Gene Ontology (GO) database (http://www.geneontology.org/), looking for significantly enriched GO terms in DEGs compared to the genome background. The formula used to calculate this was:

where N was the number of all genes with GO annotation; n was the number of DEGs in N; M was the number of all genes annotated to the GO terms and m was the number of DEGs in M. After the calculated P-value was corrected by Bonferroni correction, with the corrected P-value ≤ 0.01 as the threshold, comply this condition GO term was defined as the DEGs enriched GO term.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in DEGs was performed by comparing to the whole genome background. The calculating formula was the same as for GO analysis. Here, N was the number of all genes with KEGG annotation, n was the number of DEGs in N, M was the number of all genes annotated to specific pathways and m was the number of DEGs in M. The pathway was defined as a significantly enriched pathway when its Q value ≤ 0.05. Through pathway enrichment analysis, one can determine the most important biochemical metabolic pathways and signal transduction pathways which the significantly DEGs are involved in.

3. Results

3.1. Sequencing statistics and saturation analysis

In total, 6,150,674 and 5,917,130 sequence reads were obtained from two libraries (C and E), which corresponded to clean Tag numbers of 5,860,627 and 5,664,509, respectively. When the sequencing amount reaches 2 M or higher, the number of detected genes almost ceases to increase. In the two libraries, the saturation trend is the same (Figure 1).

Figure 1. The distribution of novel transcripts.

3.2. Tag mapping to transcript databases

The Tag sequences from C and E groups were mapped to porcine transcript data sets from the UniGene. All Tag mapping to gene number and per cent, unambiguous Tag mapping to gene number and per cent, mapping to genome number and per cent, and unknown Tag number and per cent were statistically analysed (Table 1). As shown in Table 1, there were no significant differences between the C and E groups.

Table 1. The efficiency of digital gene expression tag mapping to the transcript databases and the genome sequence.

Item	Control group	Cold stress group
All tag mapping to gene number (% of clean tags)	4763013 (81.27)	4731048 (83.52)
All tag mapping to gene number (% of distinct tags)	95833 (56.17)	90133 (58.69)
All tag-mapped gene numbers (% of ref genes)	25402 (54.50)	24971 (53.57)
Unambiguous tag mapping to gene number (% of clean tags)	3438161 (58.67)	3398164 (59.99)
Unambiguous tag mapping to gene number (% of distinct tags)	78614 (46.08)	74438 (48.47)
Unambiguous tag-mapped gene number (% of ref genes)	18952 (40.66)	18401 (39.48)
Mapping to genome number (% of clean tags)	569898 (9.72)	491863 (8.68)
Mapping to genome number (% of distinct tags)	40658 (23.83)	33938 (22.10)
Unknown tag number (% of clean tags)	387473 (6.61)	321789 (5.68)
Unknown tag number (% of distinct tags)	23663 (13.87)	20723 (13.49)

3.3. Identification of DEGs after cold stress

To increase the robustness of the data, only the TPM of more than 10 at least in one of the libraries was considered further. The transcripts with at least a 2-fold difference in expression after cold stress between two groups (FDR < 0.001) were analysed. The results indicate that after cold stress, 662 genes were significantly up-regulated and 1058 genes were significantly down-regulated. In total, 549 genes changed at least 2-fold, including 257 up-regulated genes and 292 down-regulated genes. In total, 136 genes changed at least 5-fold, including 64 up-regulated genes and 72 down-regulated genes.

3.4. Expression annotation of antisense transcripts

Sense–antisense plays an important role in gene expression regulation (Chen et al. 2005). Antisense transcripts may regulate the expression of their respective sense transcripts through diverse mechanisms, such as transcriptional interference, splicing control and mRNA stability. Sense and antisense transcripts might also form a dsRNA duplex to elicit transcriptional and/or post-transcriptional gene silencing (Ni et al. 2010). In this experiment, 1100 and 1099 Tags were mapped to the antisense transcripts in the C and E groups (TPM ≥ 10), representing 22.7% (4851) and 24.0% (4586) of expression gene, respectively. It means that there was a relative low proportion of sense–antisense regulation in the pig genome.

3.5. Novel transcripts' detection

The clean Tags which were not mapped to mRNA and mitochondria were mapped to the nuclear genome, providing a start position that can be uniquely mapped by those Tags. 1611 and 1413 Tags were mapped to the nuclear genome (TPM ≥ 10). The distribution of all novel transcripts in chromosomes is shown in Figure 1. It showed that the novel transcripts were mainly found on chromosomes 1, 2, 3, 4, 6 and 13. The distribution trends are basically the same in the two groups.

3.6. GO functional enrichment analysis

The significantly different transcripts between the C and E groups were enriched in the structural molecule activity, cytokine-mediated signalling pathway, cellular response to cytokine stimulus and small molecule metabolic process (terms from the cytoplasmic part and cytoplasm, macromolecular complex, mitochondrion, organelle inner membrane, ribosomal subunit) (corrected P-value ≤ 0.01) (Figure 2).

Figure 2. Terms from the component ontology with corrected P-value as good or better than 0.01.

3.7. Pathway enrichment analysis

In this experiment, the different transcripts enriched four pathways: ribosome, fatty acid metabolism, RNA transport and metabolic pathways.

4. Discussion

Although many organs and systems participate in resisting cold stress (Fu et al. 2013; Wang et al. 2013; Zhao et al. 2013), the most effective regulation is to generate heat quickly. Heat can be produced by all body cells through the conversion of metabolic energy into mechanical and thermal energy. But this process is very inefficient. As much as 30–70% of body heat is produced by muscle tissue through contractions (Jośe 2012). A single muscle fibre may contain 15 billion thick filaments. When muscle fibre is actively contracting, each thick filament breaks down roughly 2500 ATP molecules per second. Because even a small skeletal muscle contains thousands of muscle fibres, the ATP demands of a contracting skeletal muscle are enormous. In this experiment, we selected the skeletal muscle as the material. Because it is one of the most important heat producers, in cold stress, muscle shivering is also one of the main thermogenic patterns.

In this experiment, the result showed that after cold stress, the significantly changed genes were mainly distributed in the cytoplasm, mitochondrion, organelle envelope and ribosomal subunit. The mitochondrion is sometimes described as ‘cellular power plant', because it generates most of the cell's supply of ATP, used as a source of chemical energy (Campbell et al. 2006). The ribosomal subunit is the site of protein synthesis. The synthesis of protein promoted or suppressed after cold stress was still unknown, but at least it was found that protein synthesis had been significantly affected. There may be have some changed in cytoplasm and organelle envelope, but is unknown now, it need to further analysis.

Among the significant up-regulated genes, Thymosin β4 was the most significant changed gene, up-regulated 12.9-fold. Early research showed that Thymosin β4 could significantly induce hair growth (Philp et al. 2007), it consistent to phenomenon of Min pig rose a lot of fluff in winter. Some subunits of Nicotinamide adenine dinucleotide (NADH) were detected and changed more than 10-fold, too. NADH is a substance that helps the functionality of enzymes in the body. It plays an irreplaceable role in glycolysis, gluconeogenesis, the citric acid cycle and respiratory chain. Cold stress significantly stimulated the expression of NADH, which means the metabolism velocity increased. Among the significant down-regulated genes, voltage-dependent anion-selective channel protein 3 and glycerol-3-phosphate dehydrogenase were the most changed.

5. Conclusion

The RNA-seq of pig muscles revealed that gene expression profile changed by cold stress, especially those genes related to energy metabolism, protein synthesis and hair growth. More research of a longer term is needed to identify the genes involved in cold stress resistance.

Additional information

Funding

This work was supported by the National Nature Science Foundation of China [grant number 31301955]; Modern agricultural technology system [grant number CARS-36].