**Identification of Molecular Mechanisms Related to Pig Fatness at the Transcriptome and miRNAome Levels**

**Katarzyna Ropka-Molik 1, \* , Klaudia Pawlina-Tyszko 1 , Kacper Zukowski ˙ 2 , Mirosław Tyra 3 , Natalia Derebecka 4 , Joanna Wesoły 4 , Tomasz Szmatoła 1,5 and Katarzyna Piórkowska 1**


Received: 13 April 2020; Accepted: 27 May 2020; Published: 29 May 2020

**Abstract:** Fat deposition and growth rate are closely related to pork quality and fattening efficiency. The next-generation sequencing (NGS) approach for transcriptome and miRNAome massive parallel sequencing of adipocyte tissue was applied to search for a molecular network related to fat deposition in pigs. Pigs were represented by three breeds (Large White, Pietrain, and Hampshire) that varied in fat content within each breed. The obtained results allowed for the detection of significant enrichment of Gene Ontology (GO) terms and pathways associated directly and indirectly with fat deposition via regulation of fatty acid metabolism, fat cell differentiation, inflammatory response, and extracellular matrix (ECM) organization and disassembly. Moreover, the results showed that adipocyte tissue content strongly affected the expression of leptin and other genes related to a response to excessive feed intake. The findings indicated that modification of genes and miRNAs involved in ECM rearrangements can be essential during fat tissue growth and development in pigs. The identified molecular network within genes and miRNAs that were deregulated depending on the subcutaneous fat level are proposed as candidate factors determining adipogenesis, fatness, and selected fattening characteristics in pigs.

**Keywords:** pig; fatness; obesity; extracellular matrix; fat deposition; lipid metabolism; NGS

#### **1. Introduction**

The quantity of fat is a critical factor influencing meat quality and fattening efficiency in pig production. The fat level and mass of adipocyte tissue control food intake and body weight via leptin secretion [1,2]. Leptin is considered an essential hormone regulating adipose tissue metabolism and energy expenditure. Leptin acts as an appetite regulation factor (hunger inhibition and satiety stimulation) through interaction with hypothalamic receptors [3] and controls glucose/lipid metabolism

and body weight homeostasis through gluconeogenesis regulation [4]. Moreover, pigs can be a suitable animal model for studies of fat accumulation and obesity [5]. The identification of genetic factors determining fatness levels in pigs can be a valuable resource for medical research in humans.

In humans, genetic predisposition is a key contributing factor in obesity [6]. It has been estimated that the genetic basis of phenotypic variations in obesity can range from 40% to 70% [6,7]. In pigs, the high heritability coefficient (h<sup>2</sup> = 0.5) indicates the genetic background of fatness traits and the ability to improve them by selection [5]. The identification of genome regions, genes, or single nucleotide polymorphisms (SNPs) related to fatness in pigs has been conducted using different methodological approaches. In 2011, the usage of expression and single nucleotide polymorphism microarrays allowed for the detection of the panel of porcine genes relevant to fatness-associated traits [8]. In 2015, the genome-wide association study (GWAS) method was applied to identify QTLs (quantitative trait loci) located on Sus scrofa chromosome 7 (SSC7) and SSC4 strongly associated with growth and fatness traits based on Chinese pig breeds [9]. Guo et al. [10] using the GWAS method also indicated the involvement of the SSC7 region in the regulation of fatness and growth features. Another approach that can allow for the identification of molecular networks related to fat metabolism and growth traits is transcriptomic research concentrated on examining the hypothalamus–pituitary axis. The hypothalamus is a key brain structure controlling food intake and fat accumulation [11]. The pituitary is a crucial gland involved in fetal adipose metabolism via controlling the fatty acid synthesis and receptor-mediated lipolytic response [12]. Pérez-Montarelo et al. [13] confirmed that the usage of network interaction based on hypothalamic transcriptome analysis could pinpoint candidate genes for fatness in pigs. The transcriptomic profiling of the swine pituitary showed the strong association of this gland with postnatal growth and development [14], and also with fatness-associated traits [15].

A significant number of studies have confirmed an important involvement of adipocyte tissue in food intake [16] and body-weight regulation [17]. Such findings indicated the necessity of the analysis of adipocyte tissue to obtain the full view of metabolic processes related to fatness in pigs. To date, in pigs, transcriptomic research has been applied to identify the regulation processes determining feed efficiency [18] and to examine the body's response to different nutritional treatments [19]. Processes based on adipocyte tissue gene expression profiling in two Portuguese local pig breeds to determine the fat deposition via controlling of de novo fatty acid synthesis have been proposed, [20]. Furthermore, gene expression profiling in combination with genome resequencing has allowed for the detection of genes potentially related to backfat thickness [21].

The present study aimed to perform comprehensive whole transcriptome and miRNAome analyses [20] of adipocyte tissue to identify the molecular networks strongly related to the most crucial production traits, namely, fatness and feed intake in pigs.

#### **2. Materials and Methods**

#### *2.1. Animals, Phenotype Data Collection, and Tissue Sampling*

Animals used for gene and miRNA expression analyses were selected from a large pig population, based on fatness phenotypic traits measured after dissection. All pigs, representing three breeds (Pietrain, Pi; Hampshire, Hp; and Large White, LW) were unrelated (at least three generations back), females, and were maintained under the same housing and feeding conditions according to a procedure previously described [22]. The pigs were kept in individual pens and fed to a weight of 105 kg (±2.5 kg) (based on the Pig Test Station (SKURTCh) of the National Research Institute of Animal Production in Chorzelów). Half-carcasses were dissected 24 h after slaughter and several carcass characteristics, including fatness traits, were evaluated. The traits included average backfat thickness (ABT), calculated from five measurement points (cm) (at the thickest point over the shoulder; on the back above the joint between the last thoracic and first lumbar vertebrae; and at three points over the loin: above the rostral

edge (loin I), above the middle (loin II), and the caudal edge (loin III) of gluteal muscle cross-section), and the weight of peritoneal fat (kg), using a procedure described by Tyra and Zak [23]. ˙

Pigs for next-generation sequencing (NGS) analysis were selected based on the most important fatness parameter, i.e., ABT. In each analyzed breed (Pietrain, 8, Hampshire, 8, and Large White, 8), pigs with different fatness phenotypes were selected (4 pigs with high and 4 with low fatness). The separation of the two groups with extreme fatness within each breed allowed for avoidance of any potential breed effect on obtained genes or miRNA expression profiles (Table 1). Immediately after slaughter, the samples of fat tissue (subcutaneous fat) were collected into the RNAlater solution (Ambion, Thermo Fisher Scientific, Waltham, MA, USA) and stored at −20 ◦C.

**Table 1.** The differences in fatness traits detected in all three breeds and obtained pig groups used in genetic analyses.


L, low fatness; H, high fatness. Data are presented as means ± standard error; means with different letters (a,b differ significantly, *p*-value = 0.05).

The research did not require the approval of the Animal Experimentation Committee due to the fact that all biological material was collected from animals which were slaughtered and dissected, and, after carcass evaluation, meat was intended for consumption. The pigs were maintained and slaughtered according to Directive 2010/63/EU of the European Parliament and the Council of 22 September 2010.

#### *2.2. Whole Transcriptome Sequencing (RNA-seq)*

Extracted total RNA (Direct-zol RNA Mini Prep kit (Zymo Research, Orange, CA, USA)) was subjected to quantity and quality controls using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and a TapeStation 2200 instrument (Agilent, Santa Clara, CA, USA), respectively. Samples with RIN value (RNA integrity number) above 7 were chosen for further analyses. A quantity of 300 ng of total RNA was used for cDNA library preparation according to the TruSeq RNA Kit v2 kit protocol (Illumina, San Diego, CA, USA). The quality and quantity of cDNA libraries were checked using Qubit 2.0 and TapeStation 2200. The individual cDNA libraries were ligated with adaptors with different indexes and pooled. After pooling, each cDNA library was loaded into at least four flow cell lanes as four technical replicates. Sequencing of RNA was performed on the HiScanSQ system (Illumina, San Diego, CA, USA) using the TruSeq SR Cluster Kit v3- CBOT-HS kit and TruSeq SBS Kit v 3 - HS chemistry (Illumina, San Diego, CA, USA) according to the standard protocols and with 81 single-end reads.

The quality of raw data was analyzed using FastQC software (overrepresented sequences, sequence duplication levels, adapter contents, sequence length distribution, sequence quality scores, and GC content). The adapter's sequence reads shorter than 36 bp and/or reads with a quality score lower than 20 were removed using Flexbar software v2.5 [24]. The estimation of transcript abundance was conducted using RSEM software (v1.2.22) [25] supported by STAR aligner software (STAR\_2.4.2a) [26]. The reads were aligned to the Sus scrofa reference genome (assembly Sscrofa11.1, Ensemble release 94). The whole procedure was followed by the alignment parameter evaluation using ENCODE3's STAR-RSEM pipeline. Differentially-expressed genes (DEGs) were identified according to the following criteria: fold change ≥ |1.0| and adjusted *p*-value < 0.05 separately for each breed using DESeq2 software (version 1.12.4) [26].

#### *2.3. MicroRNA Sequencing*

MicroRNA libraries were prepared with the use of the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, MA, USA) according to the protocol. The libraries were barcoded with different indexed primers to allow multiplexing of the samples during NGS sequencing. The amplified libraries were purified and size-selected on a Novex 6% TBE PAGE gel (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). The quality and quantity of obtained libraries were measured with a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, USA) and size-assessed with a 2200 TapeStation instrument (Agilent, Santa Clara, CA, USA). The whole miRNA profile sequencing was performed in 36 cycles on a HiScanSQ (Illumina, San Diego, CA, USA) instrument according to the manufacturer's protocol.

The raw miRNA reads were quality controlled using FastQC software (parameters as for the transcriptome) [27]. The adaptor trimming and filtering were performed with the UEA sRNA Workbench v3.2 helper tool [28], while the miRCat tool v1.0 was used to detect known and novel miRNA sequences. The identification was carried out on the basis of the *Sus scrofa* genome (assembly Sscrofa 10.2) and miRBase v21.0 (Griffiths-Jones lab at the Faculty of Biology, Medicine, and Health, University of Manchester, USA) [29,30], due to the lack of 11.1 assembly in miRbase dataset. The obtained miRNAs genome localization based on 10.2 reference were converted to 11.1 reference using NCBI Genome Remapping tool. The default animal parameters were set except for minimum abundance (6 reads), minimum length (17 nt), and maximum length (25 nt). In the next step, potential other non-coding RNAs were identified and excluded employing the RNAcentral database (the RNAcentral Consortium, 2017). Moreover, isomiR-SEA software (version 1.60) [31] with the default settings was used to identify microRNA length and sequence variants (isomiRs). miRNAs differentially expressed between pigs with low and high fatness in each breed were detected using DESeq2 software (v.1.12.4) [26]. The identification of experimentally-validated target genes and enriched biological processes (KEGG, GO) was performed using the mirPath v.3.0 DIANA Tools web application (DIANA-Lab, Athens, Greece) [32], including TarBase v7.0 (DIANA-Lab, Athens, Greece) as a reference [33]. The analysis was performed based on human homologues deposited in miRBase (21.0) due to the lack of data for pig microRNAs.

#### *2.4. Validation of NGS Results*

Validation of the RNA-seq results was carried out using a real-time PCR method, for 10 selected DEGs. The primers for selected DEGs were designed based on reference sequences showed in Supplementary Table S1 and using Primer3 Input (version 0.4.0) The detailed information about the chosen genes is presented in Supplementary Table S1. The cDNA was prepared from 300 ng of total RNA using a Maxima First Strand cDNA Synthesis Kit for RT-qPCR (Thermo Fisher Scientific, Waltham, MA, USA) according to the protocol. The exact gene expression levels were measured on a QuantStudio 7 flex instrument (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA) using Sensitive RT HS-PCR EvaGreen Mix (A&A Biotechnology, Gdynia, Poland). The reaction was carried out in three replications for each gene. The expression was calculated using the delta-delta CT method, according to Pfaffl [34] and based on two reference controls *OAZ1* [35] and *RPS29* [36].

Validation of the detected expression levels of 10 selected microRNAs was carried out using the qPCR method. A quantity of 10 ng of total RNA was reverse transcribed to cDNA using a TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific, USA) following the protocol. The quantitative estimation of miRNAs was performed on a QuantStudio 7 flex instrument using TaqMan® Advanced miRNA Assays labelled with VIC or FAM fluorescent dye and with TaqMan® Fast Advanced Master Mix (Thermo Fisher Scientific, USA) according to the standard protocol (Supplementary Table S1). The relative expression levels of selected miRNAs were calculated using the delta-delta CT method according to Pfaffl [34] and based on mir-451a as a reference control. The comparison between NGS data (RNA-seq, miRNA-seq) and relative quantity obtained by the real-time PCR method was performed using the Pearson correlation (SAS software, version 8.02).

#### **3. Results**

#### *3.1. RNA Sequencing Results—Di*ff*erentially-Expressed Genes*

The whole transcriptome sequencing of adipocyte tissue allowed for the identification of differentially-expressed genes (DEGs) between pigs varying in backfat thickness. The investigation of three pig breeds enabled the detection of genes related to fatness traits regardless of breed factor. According to the NGS approach, 167 DEGs were identified between Large White pigs with different fatness phenotypes (85 up-regulated and 82 down-regulated genes in pigs with thicker backfat). A total of 32 of these were recognized as novel genes. In turn, in the Pietrain breed, 247 DEGs were detected (53 up-regulated and 194 down-regulated), of which 50 were novel genes or uncharacterized proteins. The highest number of DEGs was identified for the Hampshire breed: 128 up-regulated and 173 down-regulated genes, for a total of 301 genes (46 novel).

The comparison of the detected sets of DEGs across breeds enabled the identification of eight known genes (*SCD, SFRP2, MMRN1, PCK1, TNC, C2, CXCl10,* and *CXCL9*) and one novel gene common to all three breeds (Figure 1). The largest group of common genes identified as differentially-expressed was found for Large White and Pietrain breeds (44 DEGs and 12 novel genes). The lowest similarity in the identified gene panel was observed for Hampshire pigs and the other two breeds.

**Figure 1.** Comparison of differentially-expressed genes (DEGs) between subcutaneous fat tissue varying in thickness in each of three breeds (LW, Large White; Pi, Pietrain; and Hp, Hampshire).

α The Gene Ontology (GO) analysis performed for the DEGs identified in at least two breeds showed significantly-overrepresented GO terms (Table 2). The highest number of genes was assigned to innate immune response (GO:0045087), 14 genes; inflammatory response (GO:0006954), 11 genes; and positive regulation of apoptotic process (GO:0043065), 10 genes. The DEGs were also included in GO terms related to lipid metabolism: fatty acid biosynthetic process, long-chain fatty acid biosynthetic process, positive regulation of fat cell differentiation, positive regulation of mast cell degranulation, and response to excessive food intake (Table 2; David software). Furthermore, GO and pathways analyses using over-representation analysis (WebGestalt software) confirmed that detected DEGs were involved in the biosynthesis of unsaturated fatty acids, fatty acid metabolism, and regulation of fat cell differentiation and brown fat differentiation (not significant GO) (Figure 2). The DEGs associated with several identified GO terms were as follows: *ACACA* (acetyl-CoA Carboxylase α), *SCD* (stearoyl-CoA desaturase), *SCD5* (stearoyl-CoA desaturase 5), *FASN* (fatty acid synthase), *LEP* (leptin), and *CTGF* (connective tissue growth factor) (Table 3; Figure 3).



 FDR (false discovery rate); *p*-values are shown after Benjamini correction and according to David software; *N*, number of identified DEGs.−

\*

− − −

− − −

− − −

− − −

− −

− −

−

− −

**Figure 2.** The significant enrichment GO terms and pathways identified based on backfat thickness differences and related to lipid metabolism (WebGestalt software).

**Table 3.** The fold-change values of DEGs related to lipid metabolism and detected as significant for at least two breeds.


FDR (false discovery rate); FC, fold change; ns, not significant.

*—* **Figure 3.** The detected differentially-expressed genes involved in GO terms related to lipid metabolism processes (String software based on Sus scrofa reference; detected genes showed no more than five interactions. Line shape indicates the predicted mode of action: red, interactions that were experimentally determined; blue, interactions from curated databases; black, co-expression; green, text mining associations and interactions based on relevant publications mentioning a transfer from other organisms; yellow, transcriptional regulation).

#### *3.2. MiRNA Sequencing Results—Di*ff*erentially-Expressed miRNAs*

The comparison of whole miRNAome profiles between groups of pigs with significant-different fatness traits allowed for the identification of differentially-expressed (DE) miRNAs: 46 for LW (21 up-regulated and 25 down-regulated genes in pigs with thicker backfat), 61 for Pietrain (36 up-regulated and 25 down-regulated), and 41 for Hampshire (8 up-regulated and 31 down-regulated). As for DEGs, analogous comparison across breeds was also carried out for DE miRNAs and showed the presence of 14 miRNAs common to all three breeds (Figure 4). The lowest number of common miRNAs was detected for Hampshire and Large White pigs. Detailed data on the chromosomal localization of the identified miRNAs according to Sscrofa10.2 and Sscrofa11.1 assemblies are in the Supplementary Table S2. *—*

**Figure 4.** Venn diagram of differentially-expressed miRNAs in subcutaneous fat tissue dependent on fat content in three analyzed breeds (LW, Large White; Pi, Pietrain; and Hp, Hampshire).

The GO analysis performed using the mirPath v.3 DIANA tool showed that differentially-expressed miRNAs were involved in extracellular matrix organization (GO:0030198), extracellular matrix disassembly (GO:0022617), and innate immune response (GO:0045087), where the highest number of predicted targeted genes for miRNAs was detected. Eight miRNAs were associated with cellular lipid metabolic process (GO:0044255) (Table 4). Moreover, miR-100-5p, miR-143-3p, miR-10b-5p, and miR-24-3p were involved in fatty acid metabolism and fatty acid biosynthesis (Figure 5). The predicted targeted genes for differentially-expressed miRNAs are summarized in Supplementary Table S3.

**Figure 5.** Venn diagram of differentially-expressed miRNAs in subcutaneous fat tissue dependent on fat content in three analyzed breeds (LW, Large White; Pi, Pietrain; and Hp, Hampshire).



\**p*-values are presented after FDR correction and according to the mirPath v.3 tool; *N* target genes, the number of predicted target genes of differentially-expressed miRNAs;*N*miRNAs, number of DE miRNAs.

#### *3.3. Analysis of Pathways Common for DEGs and DE miRNAs*

The pathways analysis was performed for both DEGs and DE miRNAs identified in at least two breeds. Furthermore, the predicted genes regulated by the identified miRNAs, and involved in selected pathways, detected using mirPath v.3 DIANA tool are shown in Table 5.

An enriched pathway significant for both DEGs and DE miRNAs was fatty acid metabolism (KEGG ID: hsa01212 for miRNA/ ssc01212 for DEGs) for which four miRNAs and five genes were detected. Three of all DEGs involved in fatty acid metabolism overlapped with predicted genes regulated by detection miRNAs, namely: *SCD, FASN,* and *ACACA*. A similar situation was observed for the fatty acid biosynthesis pathway (*p*-value for miRNAs < 1.0 × 10−325; not significant for DEGs); *FASN* and *ACACA* differentially-expressed genes overlapped with miRNA-predicted targets. The pathways identified as significant and common for DE miRNAs and DEGs were also the Hippo-signaling pathway, ECM−receptor interaction, cell cycle, and P53-signaling pathways.

The pathway analysis of DE miRNAs also showed significant enrichment for the set of pathways related to toll-like receptors (toll-like receptor 2, 3, 4, 5, 9, and 10 signaling pathways; toll-like receptor TLR1:TLR2 signaling pathways; and toll-like receptor TLR6:TLR2 signaling pathways), as well as TRIF-dependent toll-like receptor signaling pathways (Figure S1).


 were used to highlighted

 the genes identified

 in both–DEGs

 and targeted genes groups.

> as

differentially-expressed

 and as predicted

 genes regulated

 by miRNAs.

 The bold and underline

SignificantpathwaysdetectedforbothDEGs(significantforleastbreeds)andDEmiRNAs(14miRNAsfor

*Genes* **2020**, *11*, 600

#### *3.4. Validation of Obtained Data*

The validation of RNA-seq results was performed for eight DEGs and six miRNAs (Table 6). For the majority of analyzed DEGs, the RNA-seq data and RQ results showed high and significant correlation coefficients (from 0.66 to 0.95), which confirmed the reliability of the NGS results. The highest similarity of results was obtained for *PCK1, ACACA*, and *LEP* genes. For four miRNAs, high correlation coefficients were detected, but these were significant only for miR-26a-5p, mir-100-5p, and mir-103a-3p.


**Table 6.** Correlation coefficients for NGS results and qPCR data for both miRNAs and DEGs.

\* *p*-value < 0.05; \*\* *p* < 0.01.

#### **4. Discussion**

In pigs, fatness traits are one of the important production features due to their strong relationship with meat quality and fattening performance. The inverse correlation between lean meat content and fatness traits results in the decrease of the fatness level in carcasses as meatiness increases [37]. Moreover, the proper fat content, including of intramuscular fat (IMF), is critical to achieving the desired meat quality parameters. The adipocyte tissue, as a secretory organ, also determines the regulation of food intake; thus, the body fat content can influence food intake and body weight [2]. In turn, the food intake and feed conversion ratio (FCR) are primary factors which determine the economic efficiency of pig production.

The high phenotypic variability of fatness characteristics, as well as the heritability coefficient (about 0.5), strongly indicate the genetic background of this group of traits and the possibility of their improvement or modification using genetic markers [5]. From a breeding perspective, the most important step would be to establish such a genetic marker that enables the controlling of pig fatness without losing the already-achieved level of meatiness. The latest successes of high throughput genetic methods, such as NGS sequencing technology and the GWAS approach, have become more applicable and create new possibilities in the search for and identification of the molecular background of phenotypic variations.

In the presented report, comprehensive whole mRNA and miRNA profiling of adipocyte tissue was applied to detect the molecular network related to fat deposition in pigs. The use of three pig breeds representing different usage types (maternal line, Large White; sire lines, Pietrain and Hampshire) allowed indication of a universal genetic basis of fatness characteristics that was not associated with any breed. The transcriptome sequencing of fat tissue delivered from pigs diverse in terms of backfat thickness enabled the pinpointing of significant, enriched Gene Ontology terms and pathways possibly related to fat deposition. For both DEGs and DE miRNAs, the fatty acid metabolism pathway was detected as significant. This pathway involved three DEGs, namely, *ACACA, FASN*, and *SCD,* and differentially-expressed miRNAs that recognized these genes as their targets. Moreover, the stearoyl-CoA desaturase gene (*SCD*) was identified as significant for all three analyzed breeds, and the *SCD* expression level was higher in pig groups with thicker fat cover. The stearoyl-CoA desaturase plays a vital role during the biosynthesis of monounsaturated fatty acids, within palmitoyland oleoyl-CoA, making them readily available to the body [38]. The higher *SCD* expression observed

in this study, in pigs with higher fat content was in accord with previously-found human data. Hulver et al. [39] showed that an increased *SCD* level is correlated to obesity. Similarly, decrease of *SCD* transcription leads to a reduction of adiposity [40]. The deficiency of stearoyl-CoA desaturase results in leanness and leads to an increase of metabolic rate as well as insulin sensitivity [38]. Moreover, a study performed on ob/ob mice indicated that *SCD* is involved in metabolic response to leptin and down-regulation of *SCD* can be a key element of the metabolic actions of leptin [41].

Reports of many authors have indicated the association of the *SCD* gene and selected fatness traits in pigs. Using the GWAS approach, Yang et al. [42] showed a possible relationship between the *SCD* locus and C18:0 content. Similarly, results from the expression Genome-Wide Association Study (eGWAS) indicated that the *SCD* gene is related to intramuscular fat composition in pigs [43]. Xing et al. [21] compared transcriptomes of subcutaneous adipose tissue of Landrace pigs depending on variable backfat levels, which showed significantly-higher *SCD* expression in animals with increased fatness. The authors proposed stearoyl-CoA desaturase as a candidate gene for fat deposition. Our results, based on three different pig breeds, confirmed previous observations, thus strongly supporting an essential role of the *SCD* gene in the determination of fat-associated traits.

Moreover, the present miRNAome sequencing identified that expression of miR-24-3p is also dependent on subcutaneous fat level. In turn, the mirPath v.3 and TarBase v7.0 DIANA bioinformatic tools found that this microRNA modulates expression of the *SCD* gene, as well as fatty acid synthase (*FASN*), acetyl-CoA carboxylase α (*ACACA)*, acetyl-coA acyltransferase 1 and 2 (*ACAA1* and *ACAA2*), and malonyl coA-acyl carrier protein transacylase (*MCAT*). miR-24-3p has been established to play a role in adipogenesis via activating adipocyte differentiation by targeting genes engaged in MAPK7-signaling pathways [44,45]. To date, only a few reports have been published concerning the possible involvement of miR-24-3p in lipid metabolism processes and obesity development. The present study indicates that this miRNA can be related to the regulation of fat deposition processes by targeting the key genes of fatty acid metabolism.

The analogues gene expression profile—increased transcript level in all fatty pigs independent of breed—was also observed for *ACACA* and *FASN* genes. A recent study applied the RNA-seq method to analyze the whole transcriptome profile of subcutaneous fat in native Puławska pigs that differed in backfat tissue [46]. The authors identified the higher transcript level of all three genes—*ACACA, FASN*, and SCD—in pigs with increased backfat thickness. Previous research showed the significant association of the *FASN* gene with meat quality and the *ACACA* gene with IMF content. In Iberian pigs, the *ACACA* gene was proposed as a candidate gene responsible for intramuscular content of saturated fatty acids and monounsaturated fatty acids fatty acids [47]. Furthermore, the up-regulation of *FASN* and *SCD* genes were reported in Alentejano pigs and related to the higher fat deposition observed in this breed [20]. The research performed on Iberian × Landrace pig cross lines indicated a significant association of the *ACACA* gene and IMF palmitic fatty acid percentage [48]. Similarly, the *FASN* gene was previously correlated to meat quality, fatty acid content, and composition [49,50]. Our research, supported by previous studies, strongly confirms that the identified DEGs determine fat deposition in pigs.

#### *Enriched Metabolic Process and Pathways Associated with Fat Deposition*

Genes whose expression varied between thick and thin backfat groups belonged to several processes and pathways directly related to fatty lipid metabolism, within the fatty acid biosynthetic process, the long-chain fatty acid biosynthetic process, positive regulation of fat cell differentiation, and response to dietary excess. The molecular network of food intake controlled by adipocyte tissues has been widely investigated and described. One of the main factors determining the energy balance by hunger inhibition is leptin [51]. The DEGs set comparison showed that *LEP* gene expression was increased in fat tissue of more obese pigs independent of breed. These findings are in accordance with observations made in humans, where the overexpression of the leptin gene in subcutaneous and omental adipose tissues in obese patients was detected [52]. Adipocyte cells accumulate triglyceride

and increase the synthesis of leptin during their growth. Leptin acts via the hypothalamus to control energy balance and feeding behavior [53] Thus, leptin, also called an anti-obesity hormone, reduces food intake and increases energy expenditure [54]. Likewise, the other two DEGs identified in the present study—*PPARGC1A* (peroxisome proliferator-activated receptor-γ co-activator-1alpha) and *PCSK1N* (proprotein convertase subtilisin/kexin type 1 inhibitor)—have been associated with obesity in humans. A polymorphism in the *PPARGC1A* gene is related to obesity and type 2 diabetes [55], while the down-regulation of *PCSK1N* gene expression is also associated with obesity [56]. In the present research, *PPARGC1A* expression was up-regulated in thicker backfat tissue in all three investigated pig breeds, whereas expression of the *PCSK1N* gene was lowered in obese pigs represented by Hampshire and Large White breeds. This suggests that all identified DEGs related to mast cell (MC) degranulation can affect the increase of fat tissue mass and lipid metabolism.

Research performed in humans and animals strongly indicates that MCs are involved in activation of adipocytes and recruitment of inflammatory cells [57]. The positive regulation of the MC degranulation-GO term was identified as significantly deregulated in fat tissue dependent on fat cover thickness. The research allowed for the detection of *FGR, FCER1A, FCER1G,* and *ZAP70* as DEGs associated with MC degranulation. The function of MC degranulation is associated with the release of multiple enzymes which can influence adipocyte tissue modification. In vitro studies on human white adipocyte tissue showed that MCs can promote adipose initiation in response to cold [58]. Furthermore, MC cells induced remodeling of adipose tissue extracellular matrix [57].

Interestingly, the obtained data showed the significant enrichment of ECM organization (GO:0030198) and ECM−receptor interaction pathways. The panel of DEGs was identified including genes coding for collagens, thrombospondins, and laminin, which belong to extracellular matrix remodeling. A recent report highlighted that ECM plays a critical role in adipose tissue expansion through controlling cell migration during body growth and development [59]. It has been proven that adipose tissue expansion is strongly related to ECM remodeling at the level of matrix individual components—collagens, thrombospondins, metalloproteinases, and their inhibitors (Tissue Inhibitor of Metalloproteinase - TIMPs) [60]. The present research allowed the identification of DE miR-145-5p and its target gene, tenascin C (*TNC*), both of which are involved in the ECM−receptor interaction pathway. Tenascin C, as a glycoprotein member of ECM, is related to tissue developmental processes, while miR-145 is involved in abdominal obesity, insulin resistance, and lipid metabolism [61]. The exact role of miR-145 has not been well established but the increased expression of miR-145 was observed in adipose and liver tissues in obese patients [62]. It has been proposed that this miRNA plays a role in lipolysis and our results supports this thesis, and also confirm that miR-145 overexpression in fat tissue is associated with obese individuals, which was independently observed in the investigated pig breeds.

Another interesting observation is the down-regulation of the *TNC* gene identified in pigs with thicker fat cover. Tenascin C determines the formation of the collagen network in adipocyte tissue controlling cell migration and proliferation [63]. On the other hand, the excessive cross-linking of adipocytes by fibrotic components can lead to the reduction of their metabolic activity [64]. We hypothesized that the increased expression of the *TNC* gene in thinner backfat can be related to significant modification of extracellular matrix components, which may result in reduced adipocyte activity and a decreased rate of fat tissue development. Tenascin C is also related to activation of the toll-like receptor 4 (TLR4) signaling pathway, which triggers the obesity-induced inflammatory response [65]. Our results also showed the significant enrichment of both toll-like receptor TLR1:TLR2 signaling and toll-like receptor TLR6:TLR2 signaling pathways for DE miRNAs. Gene expression profiling of human peripheral blood mononuclear cells (PBMCs) showed the increased transcript level of *TLR2* and *TLR4* genes in obese patients [66]. The increase in expression of the *TLR4* gene, reported in differentiating adipocytes in db/db mice suggests that this gene is critical during obesity development processes [67]. In the present study, in addition to the detection of GO terms related to *TLR2* and *TLR6*, the set of DE miRNAs targeting the *TLR4* gene was identified, namely, the members of the let-7 family

(let-7a-5p, let-7i-5p, and let-7d-5p). To date, many reports have indicated a significant role of let-7 miRNAs in regulation of adipogenesis, lipolysis, and insulin resistance [68]. The differential expression of several let-7 family members may confirm that, beside *TLR* genes, this miRNA family plays a key role in the fat deposition rate and lipid metabolism.

The other proven role of ECM receptor pathways during fat tissue development is the regulation of adipocyte apoptosis and/or necrosis processes [60]. The results obtained in our study also indicated expression changes of genes involved in the positive regulation of apoptotic process GO term and cell cycle pathways. Other members of the EC matrix—thrombospondins—such as those identified in our study (*THBS2, THBS3,* and *THBS4*), are mainly detected in visceral adipose tissues and their increased expression level is associated with obesity. Moreover, thrombospondin 1 has been proposed as a novel marker related to obesity and metabolic syndrome [69]. The identification of the panel of genes involved in ECM receptor pathways and extracellular matrix assembly confirmed that such signaling is essential for fat deposition processes in pigs irrespective of a breed factor.

#### **5. Conclusions**

The backfat thickness and growth rate are closely related to pork quality and fattening efficiency, as well as reproductive and immune characteristics. Taking into account whole transcriptome profiling of pigs varying in fat content across three pig breeds, the significant enrichment of Gene Ontology terms and pathways associated directly and indirectly with fat deposition was detected (*ACACA, ACOX3, FASN, SCD5, SCD,* and *PLPP1).* The differentially-expressed miRNAs and genes were involved in fatty acid metabolism, positive regulation of fat cell differentiation, and the inflammatory response. Moreover, the results showed that adipocyte tissue content regulated the expression of leptin and other genes related to a response to dietary excess. Our results confirmed previous findings in humans that showed ECM organization and disassembly are fundamental for fat tissue growth and development. The modification of genes and miRNAs involved in ECM rearrangements can also be essential during fat tissue growth and development in pigs. We also suggest that mast cell degranulation can be one of the significant processes associated with adipocyte tissue development. The pinpointed molecular networks within genes and miRNAs deregulated by subcutaneous fat level are proposed as candidate factors determining adipogenesis, fatness, and selected fattening characteristics in pigs. Identified DEGs and DE miRNAs should be investigated in the future in terms of their use as genetic markers associated with pig production traits.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/6/600/s1, Table S1. Detailed information about DEGs and miRNAs whose expression was estimated using real-time PCR methods. Table S2. Detailed data on the chromosomal localization of the identified miRNAs according to Sscrofa10.2 and Sscrofa11.1 ("LW" stands for Large White samples, "Pi" denotes Pietrain samples, and "Hp" stands for Hampshire samples; "NA" denotes potentially novel miRNA identified in this study). Table S3. TarBase experimentally supported interactions—the predicted targeted genes (N) for each of the 14 detected miRNAs. Figure S1. The significant enrichment molecular pathways in which 14 microRNAs commonly detected in all three pig breeds were detected (DIANA-miRPath v3.0, DIANA-Lab, Athens, Greece). Figure S2. The significant enrichment Hippo-signaling pathway (KEGG database, hsa04390) involved genes targeted for identified differentially-expressed miRNAs. Targeted genes are marked orange (genes included in more than one pathway) and yellow (genes included only in one pathway). Figure S3. The significant enrichment cell cycle (KEGG database; hsa04110) involved genes targeted for identified differentially-expressed miRNAs. Targeted genes are marked orange (genes included in more than one pathway) and yellow (genes included only in one pathway).

**Author Contributions:** For research articles with several authors, a short paragraph specifying their individual contributions must be provided. The following statements should be used "Conceptualization, K.R.-M., K.P.-T., and K.P.; methodology, K.R.-M. and K.P.-T.; software, K.Z. and T.S.; validation, K.P. and K.R.-M.; formal analysis, ˙ K.R.-M., K.P.-T., N.D., and J.W.; investigation, K.R.-M., K.P.-T., N.D., and J.W.; resources, M.T.; data curation, K.Z. ˙ and T.S.; writing—original draft preparation, K.R.-M.; writing—review and editing, K.P.-T., K.P., K.Z., T.S., M.T., ˙ N.D., and J.W.; visualization, K.R.-M.; supervision, K.R.-M.; project administration, K.R.-M. and K.P.; funding acquisition, K.R.-M. and K.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by statutory activity of the National Research Institute of Animal Production, grant number 04.-0.14.1 (cost of biological material collection and storage) and grant number 01-18-05-21 (cost of genetic analyses).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


muscle contributes to abnormal fatty acid partitioning in obese humans. *Cell Metab.* **2005**, *2*, 251–261. [CrossRef]


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Genome-Wide Transcriptome and Metabolome Analyses Provide Novel Insights and Suggest a Sex-Specific Response to Heat Stress in Pigs**

**Krishnamoorthy Srikanth 1,**† **, Jong-Eun Park 1,**† **, Sang Yun Ji 2 , Ki Hyun Kim <sup>2</sup> , Yoo Kyung Lee 2 , Himansu Kumar 1 , Minji Kim <sup>2</sup> , Youl Chang Baek 2 , Hana Kim <sup>1</sup> , Gul-Won Jang 1 , Bong-Hwan Choi <sup>1</sup> and Sung Dae Lee 2, \***


Received: 5 March 2020; Accepted: 6 May 2020; Published: 11 May 2020

**Abstract:** Heat stress (HS) negatively impacts pig production and swine health. Therefore, to understand the genetic and metabolic responses of pigs to HS, we used RNA-Seq and high resolution magic angle spinning (HR-MAS) NMR analyses to compare the transcriptomes and metabolomes of Duroc pigs (*n* = 6, 3 barrows and 3 gilts) exposed to heat stress (33 ◦C and 60% RH) with a control group (25 ◦C and 60% RH). HS resulted in the differential expression of 552 (236 up, 316 down) and 879 (540 up, 339 down) genes and significant enrichment of 30 and 31 plasma metabolites in female and male pigs, respectively. Apoptosis, response to heat, Toll-like receptor signaling and oxidative stress were enriched among the up-regulated genes, while negative regulation of the immune response, ATP synthesis and the ribosomal pathway were enriched among down-regulated genes. Twelve and ten metabolic pathways were found to be enriched (among them, four metabolic pathways, including arginine and proline metabolism, and three metabolic pathways, including pantothenate and CoA biosynthesis), overlapping between the transcriptome and metabolome analyses in the female and male group respectively. The limited overlap between pathways enriched with differentially expressed genes and enriched plasma metabolites between the sexes suggests a sex-specific response to HS in pigs.

**Keywords:** heat stress; pig; Duroc; RNA-Seq; NMR; metabolome

#### **1. Introduction**

Heat stress affects animal husbandry worldwide and is a major environmental factor that affects animal health and production [1,2]. The increasing environmental temperature due to global warming will particularly affect pigs due to their lack of functional sweat glands, which would have helped in endogenous heat dissipation [3,4]. Added to this, the thick subcutaneous adipose tissue in the pig, impedes effective radiant heat loss [5]. Moreover the increase in metabolic rate due to rapid lean tissue accretion increases endogenous heat production, exacerbating the innate inability of porcine animals to tolerate heat [3,6–9]. The estimated annual economic loss due to heat stress for the swine industry in the US alone is nearly US \$300 to \$450 million [5,10,11], mainly due to increased mortality and morbidity, altered carcass composition, decreased feed efficiency, inconsistent growth, reduced fecundity and poor sow performance [6], so heat stress in a pressing issue for worldwide pig production.

Quiniou and Noblet [12] have shown that the growth rate is normally between 18 and 25 ◦C and that the thermoregulatory response in pigs is activated at about 25 ◦C. Animals respond to heat stress by regulating physiological and metabolic changes, such as redistribution of blood flow from the body core to the periphery, and by reducing feed intake [13] to reduce metabolic heat production [14]. This has a significant effect on animal health and productivity [3]. The heat stress mitigation options available at present include expensive heat abatement processes, such as spray or floor cooling and nutritional supplementation strategies. Since the response to heat stress varies within a population, due to genetic variations, selection for thermal tolerance could be a better alternative [1,15]. However traditional breeding programs have limitations due to the difficulties in collecting precise records of thermal states of individual animals. To overcome this, production traits that are correlated with heat stress are used for selecting animals with higher thermal tolerance [16]. They, however, mask the actual effects of heat stress; therefore, a precise and efficient method for identifying genetic tolerance to heat stress is necessary. Identifying genetic and metabolic markers that are sensitive to heat stress will be the most important step for managing heat stress in pigs. These markers could also be putative candidates for identifying genomic variants for improving heat tolerance in pigs.

Previous studies have shown that heat stress triggers a complex array of gene expression and metabolic changes in livestock, resulting in the modulation of a wide variety of pathways involved in the heat stress response, the inflammatory response, DNA damage repair, chaperone use, etc. [5,17–19]. Similarly, a variety of metabolites involved in carbohydrate, amino acid, fatty acid and amine metabolism are modulated, including glucose, lactate, alanine, cysteine, isoleucine, etc. [20–23]. However, these studies did not involve integrated transcriptome and metabolome profiling, which could lead to better understanding of the effect of heat stress on swine physiology, gene expression and metabolism, which will help in designing and implementing innovative strategies to limit the economic losses due to heat stress on pig farm profitability. Therefore, in this study, we exposed 3-month-old pigs to heat stress and performed transcriptome and metabolome analysis using RNA-Seq and high resolution magic angle spinning (HR-MAS) NMR (nuclear magnetic resonance) to understand the gene expression and metabolome perturbations in response to heat stress.

#### **2. Materials and Methods**

#### *2.1. Animals*

The experiments were performed following the ethical guidelines laid down by the Institutional Animal Care and Use Committee (IACUC) of the National Institute of Animal Science (NIAS). The experimental procedures were reviewed and approved by IACUC of NIAS (NO 2017-1070). A total of 6 (3 barrows and 3 gilts) unrelated, 3-month-old Duroc pigs, weighing 56 ± 2.62 Kg, were used in this study. Animals were fed a corn and soybean meal based diet (Table 1), containing 18% crude protein, 4.9% crude fat, 4.6% crude fiber and 4.4% crude ash; the whole diet was formulated to provide 3450 Kcal/Kg. The animals had ad libitum access to feed and water. Water and feed intake were measured (Table S1). One-way analysis of variance (ANOVA) was performed in R to test for significance. Differences were considered to be significant at *p* < 0.05. Other than the typical heat stress responses, no significant differences in behavior were noted.


**Table 1.** Composition of the total mixed ration (TMR) used during the experimental period.

#### *2.2. Heat Stress Experimental Setup*

The heat stress (HS) experiment was carried out in an environmental chamber which had controls for temperature and humidity in NIAS. After acclimation to the chamber for 3 days, the experiments were performed. The animals were first kept at 25 ◦C and 60% relative humidity for 24 h, and at the end, 10 mL of blood was drawn via venipuncture and 5 mL of blood was stored in a Tempus Blood RNA tube (Life Technlogies, Carlsbad, CA, USA) for RNAseq analysis and the rest of the blood was stored in a vacutainer tube containing anticoagulants (EDTAK2), for the metabolomic analysis. The same animals were then exposed to 33 ◦C temperature and 60% relative humidity for 24 h and blood was drawn and stored in the same way.

#### *2.3. RNA Isolation and Sequencing*

The RNA was isolated using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's guidelines. The integrity and quality of the RNA was assessed using 2100 Bioanalyzer and RNA 6000 Nano LabChip Kit (Agilent Technologies, Palo Alto, CA, USA). Only samples with RIN values >8 were used for library construction. The concentration of RNA was determined using NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The sequencing library was constructed using Illumina TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA) following processes previously described [16]. The sequencing was performed on an Illumina HiSeq 2000 sequencer. The raw reads are available for download from sequence read archive (SRA), NCBI under the accession number SUB7043048.

#### *2.4. Sample Preperation and Metabolite Analysis with <sup>1</sup>H-NMR*

36 µL of Plasma was mixed with 4 µL of D2O (20 mM TSP-d4), and the sample was then transferred to a 4 mm NMR nanotub, and <sup>1</sup>H-NMR spectra were obtained by high-resolution magic-angle-spinning (HR-MAS) NMR spectra using an Agilent 600 MHz NMR spectrometer (Agilent Technologies, Palo Alto, CA, USA) with a 4 mm gHX Nanoprobe. A Carr–Purcell–Meiboom–Gill (CPMG) pulse sequence was used to reduce signals generated by macromolecules and water. The <sup>1</sup>H-NMR spectra were measured as described previously [24]. The TSP-d<sup>4</sup> peak at 0.0 ppm was used as a reference for calibrating chemical shifts.

#### *2.5. Data Analysis*

#### 2.5.1. RNA-Seq Data

The processing of the data was previously described [25]. Briefly, after checking the quality of the raw reads using FastQC (version 0.11.5) [26] and trimming the adaptors and low-quality bases using TRIMMOMATIC (version 0.36) [27], the reads were aligned to the Pig reference genome (*Sus scrofa* 11.1) using HiSAT2 (version 2.05) [28]. The aligned reads were then counted using FeatureCounts (version 1.5.0) [29]. After correcting for batch and unknown effects using Svaseq [30], differential expression analysis was performed using DESeq2 [31]. Significant genes (FDR < 0.1) were identified, and functional enrichment analysis based on Gene Ontology (GO) under Biological Process and Molecular Function was performed with DAVID [32]. KEGG (Kyoto Encylopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using ClueGO plugin [33] in cytoscape version 3.7.2 [34]. Goseq [35] was used for metabolic pathway enrichment analysis.

#### 2.5.2. Metabolome Data

The metabolite quantification was performed using Chenomx NMR suite 7.1 (Chenomx, Edmonton, Canada). The generated spectra were binned with a binning size of 0.001 ppm, and normalized to the total area, and binned data were aligned with the icoshift algorithm of MATLAB R2013b (MathWorks, Natick, MA, USA). Principle component analysis (PCA) and statistical analyses were performed using MetaboAnalyst 4.0 [36]. Only features that were detected in at least 50% of samples were used.

#### 2.5.3. Pathway Enrichment Analysis of Differentially-Enriched Metabolites

Significantly differentially-enriched metabolites (*p* < 0.05) were subjected to KEGG pathway analyses using the Pathway analysis module in MetaboAnalyst 4.0 [36]. A combination of quantitative enrichment and topology analysis using only curated metabolic pathways from the KEGG database was used in the analyses.

#### *2.6. Real-Time PCR Validation*

Real time reverse transcriptase PCR (qRT-PCR) was performed using gene-specific primes (Supplementary File S1). The PCR was performed on an ABI 7500 Real Time PCR system using Fast SYBR green master mix (Applied Biosystems, Foster City, CA, USA). A total of 18 genes (9 each from male and female differentially expressed gene (DEG) set) were analyzed. β-actin and GADPH (Glyceraldehyde-3-phosphate dehydrogenase) were used as endogenous controls. The stability of expression for each of those two genes was checked using GeNORM (https://genorm.cmgg.be/) against the same concentration of RNA from different samples. β-actin was the most stable and was used for normalizing the expression data of the target genes.

#### **3. Results**

#### *3.1. Transcriptome Alignment, Mapping and Principle Component Analysis*

We investigated the impact of heat stress on Duroc pigs using high throughput RNA-Seq analysis. A total of 423.3 million 100 bp Paired-End (PE) reads corresponding to an average of 35.2 million reads per individual was generated. After trimming for adapters and low-quality reads, 410.5 million reads remained. The reads were mapped to the *Sus scrofa* reference genome at an average alignment rate of 96.8% (File S2). The reads were mapped to a total of 19,283 genes. Principal component analysis (PCA)

(Figure S1) suggested that sex had a large effect on transcriptome difference between the groups, so we decided to compare the heat stress effect on male and female pigs separately. PCA showed that 30% and 36% of the expression variation was due to heat stress in female and male pigs respectively (Figure 1a,c); however, considerable within group variation was also observed, confirming previous reports that the heat stress response varies within populations due to underlying genomic variation [1,15].

**Figure 1.** Summary of the transcriptome analysis: (**a**) PCA of female group samples; (**b**) volcano plot showing the number and distribution of significantly differentially expressed genes in the female group; (**c**) PCA of male samples; (**d**) volcano plot showing the number and distribution of significantly differentially expressed genes in the male group; Venn diagrams showing differentially expressed genes DEGs common between male and female groups: (**e**) up-regulated and (**f**) down-regulated.

#### *3.2. Porcine Transcriptome Response to Heat Stress*

Heat stress resulted in 552 and 879 genes being significantly (FDR < 0.1) differentially expressed in male and female groups respectively. Out of those, 236 and 540 genes were up-regulated and 316 and 339 genes were down-regulated in female and male pigs respectively (Figure 1b,d, File S3). The genes modulated in response to heat stress were considerably different between the male and female groups, with only 19 up-regulated and 11 down-regulated genes being common between the groups (Figure 1e,f). Among the top up-regulated DEGs in the female group were *HIST2H2AA4 (*histone cluster 2H2A

family member a4), *MYH4* (myosin heavy chain 4), *ABCA6* (ATP binding cassette subfamily A member 6), *HSF4* (heat shock factor 4) and *IL18* (interleukin 18), while top down-regulated DEGs included *CXCL10* (C-X-C motif chemokine 10), *IDO1* (indoleamine 2,3-dioxygenase 1), *HES4* (Hes family BHLH transcription factor 4) and *IFI6* (γ-interferon-inducible protein Ifi-16); in the male group the up-regulated genes included *DES* (Desmin), *RAG1* (recombination activating 1)*, LSAMP* (limbic system associated membrane protein), *CD1E* (CD1e molecule), *HSPG2* (heparin sulfate proteoglycan 2) and down-regulated genes included *CEBPE*(CCAAT/enhancer binding protein), *ZNF316* (Zinc Finger Protein 316)*, PXDC1*(PX domain containing 1) and *COX4I1* (cytochrome c oxidase subunit 4l 1).

Gene Ontology (GO) enrichment analysis of the significantly DEGs (FDR < 0.1) in the female group (Figure 2a) showed that among the up-regulated DEGs, "angiogenesis," "apoptosis," "response to heat," " heme biosynthetic process" and "ATPase activity" were enriched, whereas "negative regulation of innate immune response," "positive regulation of signal transduction," "response to cytokine," "MAPK cascade," "response to lipopolysaccharide," "threonine-type endopeptidase activity," etc., were enriched amongst the down-regulated genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of the female group DEGs (Figure 3a) showed that "apoptosis," "non-small cell lung cancer," "Toll-like receptor signaling pathway," "necroptosis," "NOD-like receptor signaling pathway" and "Influenza A" pathways were enriched.

**Figure 2.** Gene Ontology analysis: (**a**) enriched GOs in female group; (**b**) enriched GOs in male group. The GO enrichment was performed under Biological Process and Molecular Functions categories.

In the male group, heat stress resulted in the up-regulation of genes functioning in "inflammatory response," "innate immune response," "cell adhesion," "integrin-mediated signaling pathway," "cell cycle arrest," "positive regulation of angiogenesis," "apoptosis," "response to oxidative stress," "cellular response to heat" and "glycoprotein binding." Genes functioning in the "cellular protein catabolic process," "regulation of endocytosis," "negative regulation of astrocyte differentiation," "cytochrome-c oxidase activity" and "N-acetylgalactosamine-6-sulfatase activity" were down-regulated (Figure 2b). Among KEGG pathways (Figure 3b), genes part of "complement and coagulation cascades," "platelet activation," "fc γ r-mediated phagocytosis," "thyroid hormone signaling," "pathways in cancer," "proteoglycans in cancer," and "Notch signaling" pathways were enriched. The top 10 DEGs from the female and male group were validated using q-RT PR analysis (quantitative Real-time PCR) (Figure 4); the correlation between RNA-Seq and q-RT PCR analyses was 0.861.

**Figure 3.** KEGG pathway enrichment analysis: ( **a**) enriched KEGG pathways in female group; ( **b**) enriched KEGG pathways in male group. Nodes are genes and edges represent pathways. Up-regulated genes are represented as rectangles and down-regulated genes are represented as hexagons.

‐ ‐ ‐ ‐ ‐ ‐ **Figure 4.** qRT-PCR validation of RNA-Seq results, using the top 10 DEGs from male and female groups respectively. The fold change from q-RT PCR was plotted on the *x*-axis, and the log<sup>2</sup> fold change from the RNA-Seq analysis was plotted on the *y*-axis; and the correlation between the two methods is given in the figure.

#### *3.3. Porcine Metabolome Response to Heat Stress*

 ‐ ‐ We detected and identified a total of 46 metabolites (Figure 5b,d, File S4). After normalizing the sample (quantile normalization) and scaling the data (auto scaling) in MetaboAnalyst 4.0 [36], PCA and a heatmap were generated. The overall metabolic profiles of the groups resolved well in PC1 (Figure 5a,c), and hierarchical clustering based on metabolite concentration showed that the control and heat stress group were well separated. However, PC2 showed that there was considerable within-group variation in both male and female pigs. The fold change was calculated as a ratio between the two groups' means. A Wilcoxon rank-sum test was performed, and 30 and 31 metabolites were found to significantly differ (*p* < 0.05) between the control and heat stressed groups in female and male groups respectively (File S4).

‐ ‐*δ*‐ ‐ ‐ ‐ In the female group, out of the 30 metabolites that significantly differed between the heat stress and control groups, the concentrations of 21 metabolites increased while nine metabolites decreased in concentration. The concentrations of tryptophan, arginine, pyruvate, alanine, acetate and lactate were higher, while DL-plus *allo*-δ-hydroxylysine, creatine, urea, 3-aminoisobutric acid, 4-aminobutyric acid and 3-methylhistidine were lower in the heat stress group (File S4). In the male group samples, 18 of the 31 significantly differing metabolites increased, among which where aspartic acid, pyruvate, betaine, creatine and lactate, while 13 metabolites, including tryptophan, carnosine, valine, alanine and creatine, were significantly reduced in the heat stress group compared to the control group (File S4).

**Figure 5.** Metabolome analysis: (**a**) metabolic profiles of the female group visualized with principal component analysis; (**b**) heat map of all metabolites identified in the female group; (**c**) metabolic profiles of the male group visualized with principal component analysis; (**d**) heat map of all metabolites identified in the male group.

#### *3.4. Metabolic Pathways Enriched in Response to Heat Stress*

‐

Pathway analysis was performed with metabolites that significantly differed between the heat stress and control group using the pathway enrichment module in MetaboAnalyst 4.0 [36]. Twelve pathways were significantly (*p* < 0.05) enriched in the female group (Figure 6a), including "arginine and proline metabolism," "glycine, serine and threonine metabolism," "cysteine and methionine metabolism," "alanine, aspartate and glutamate metabolism" and "histidine metabolism". We used goseq [35] to identify metabolic pathways in KEGG database that were significantly altered in the transcriptome, in response to heat stress. Among the twelve pathways significantly enriched (Figure 6b), four pathways were found to be commonly enriched with DEGs and DEMs (differentially enriched metabolites) (Figure 6c); these were "arginine and proline metabolism," "glutathione metabolism," "selenoamino acid metabolism" and "histidine metabolism." Out of these, the first three pathways have been found to be enriched in response to oxidative stress induced by chronic HS in the pig [37–39].


**Figure 6.** Metabolic pathways altered in the female group: (**a**) Metabolic pathway enrichment analyses of metabolites enriched in response to heat stress. (**b**) Metabolic pathway enrichment analyses of genes differentially expressed in response to heat stress. (**c**) Number of overlapping pathways enriched with DEGs and significantly differing metabolites in the heat stress group.

Among the male animals, 10 metabolic pathways were enriched with differently enriched metabolites, including "arginine and proline metabolism," "glycine, serine and threonine metabolism," "cysteine and methionine metabolism," "taurine and hypotaurine metabolism" and "histidine metabolism" (Figure 7a). The metabolic pathways enriched with DEGs (Figure 7b) included "inositol phosphate metabolism," "amino sugar and nucleotide sugar metabolism," "pantothenate and COA biosynthesis" and "riboflavin metabolism." Comparative analysis of the metabolic pathways enriched with DEMs and DEGs showed that three metabolic pathways were common (Figure 7c); those were "pantothenate and CoA biosynthesis," "sulfur metabolism" and "arginine proline metabolism." Pantothenate and CoA biosynthesis was previously found to be enriched in response to HS in the pig [5,20].


**Figure 7.** Metabolic pathways altered in the male group: (**a**) Metabolic pathway enrichment analyses of metabolites enriched in response to heat stress. (**b**) Metabolic pathway enrichment analyses of genes differentially expressed in response to heat stress. (**c**) Number of overlapping pathways enriched with DEGs and significantly differing metabolites in the heat stress group.

#### **4. Discussion**

Heat stress occurs when the rate of heat accumulation exceeds the rate of heat loss; i.e., when the animal is removed from the thermal comfort zone. HS causes the core body temperature of a pig to rise [40]. Temperature increase due to global warming will be detrimental to swine health and production. The global temperature has already been raising at an average of 0.13 ◦C over the past 50 years [41]. Pigs are particularly vulnerable due to their lack of functional sweat glands, and the thick layer of subcutaneous adipose tissue that they possess which impedes effective radiant heat loss [42]. Heat stress triggers an adaptive stress response by regulating a complex array of genes and metabolites [16,18,21]. In this study we found significant differences in water and feed intake post heat stress, and a significant difference in water intake post heat stress treatment between male and female pigs (Table S1).

‐ ‐ ‐ ‐ Several studies have suggested that gender might have an effect on the heat stress response in pigs [40,43–45]. The sex-specific response to heat stress could be due to the difference in surface area to mass ratio or due to body compositional difference (lean meat to fat mass) which affects radiant heat loss [43]. Barrows are reported to generally have greater backfat thickness than gilts [46]. Though barrows consume more feed and gain body weight more rapidly, gilts deposit proportionally more muscle and less fat in their carcass [43,47]. This sex-specific difference in body composition could significantly alter the heat stress response, as it affects the ability of the animal to maintain the core body temperature through radiant heat loss. Moreover, several hormones, such as prostaglandins [48], endogenous opioids [48] and glucocorticoids have roles in the thermoregulatory mechanism [49]. Cortisol, a glucocorticoid, is the mostly widely used biomarker for detecting stress [50,51]. The concentration of cortisol is 15% higher in barrows than in gilts [52]. HS-induced increases in the concentrations of catecholamines and glucocorticoids have been observed in many species [53,54]; the effect of HS on immune cells is suggested to be dependent upon these hormones [55–57]. Similarly, Rudolph et al. [40], found that oxidative stress did not result in muscle injury in barrows, while previous studies in gilts

‐

reported that oxidative stress causes muscle injury [58–61]. This led them to suggest that muscle response to heat stress could be partially dependent upon sex of the pig, and that muscle tissues in males could be more resistant to heat stress than in females. The result of our study concurs with all these reports; our study shows that HS in pigs, as evidenced from transcriptome expression and metabolome concentration, is sexually dimorphic.

#### *4.1. Transcriptome Regulation in Response to Heat Stress*

Transcriptome analysis showed that heat stress caused a profound shift in gene expression in the male group: 879 genes were found to be significantly differentially regulated in the male group, while 552 genes were DEs (differentially expressed) in the female group. Among the DEGs were several genes involved in Heat shock response. Heat shock proteins (HSPs) are expressed upon encountering stress [62,63]. Heat is proteotoxic and causes proteins to denature [64], HSPs act as chaperones in assisting in protein folding, thereby helping in avoiding protein aggregation [65] and maintaining cellular homeostasis [66]. *HSPA5* (heat shock 70 kDa protein 5)*, HSP90AA1* (heat shock protein 90 α family class a member 1) and *HSF4* (heat shock transcription factor 4) were up-regulated, in both male and female groups, while *FKBP5* (FKBP prolyl isomerase 5) was up-regulated only in the male group. *HSPA5* is a member of the HSP70 (heat shock protein 70) family of genes; it serves as an endoplasmic reticulum chaperone and as a sensor of protein misfolding [67]. *HSP90AA1* a member of the HSP90 family takes part in several cellular functions such as regulating protein activity, transport and also in activating different signaling pathways by forming complexes with steroid receptors [68]. *HSP90AA1* is essential for normal spermatogenesis in pigs [69], and regulation of *HSP90AA1* protects cells from heat shock [70]. *FKBP5* a member of immunophilin family, is a glucocorticoid receptor (GR)-regulating co-chaperone of heat shock protein 90, their expression are regulated by progestins and glucocorticoids [71], and are significantly correlated with plasma cortisol concentration in pig [72,73], suggesting that the male pigs were under considerable stress.

Several co-chaperones and genes involved in cellular response to heat stress were also DEs; these included *HMOX1* (heme oxygenase 1)*, TRPM2* (transient receptor potential melastatin 2), which were up-regulated in the male group, *IL1A* (interleukin 1 α); up-regulated in both male and female groups, *SOD1* (superoxide dismutase 1) down-regulated in male group and TFEC (Transcription factor EC), DAXX (death domain associated protein), and TRPM2 which were down-regulated in the female group. *HMOX1* is an oxidative stress marker [74], and is up-regulated in response to oxidative stress [75]. *SOD1* catalyzes the removal of superoxide radicals that are generated due to biological oxidation [76]. *SOD* expression was found to be down-regulated in porcine skeletal muscle in response to long-term (3 days) heat stress [58]. In this study SOD expression was down-regulated in the male group following heat stress; this along with the increased expression of *HMOX1* suggests that the male group animals may have experienced significant oxidative stress following heat stress. Several genes involved in response to oxidative stress were up-regulated in the male group, including *TRPM2* [77], *LRRK2* [78], *MAPK14* [79] and *VNN1* [80]. In the female group *PRDX2* [81], *SLC7A11* [82] and *FOXO3* [83] which protects cells from ROX were up-regulated, while *TREX1* involved in cellular response to oxidative stress [84] and *SETX* (Senataxin) involved in defense against DNA damage due to oxidative stress [85] were down-regulated in the female group.

Immune response is closely linked with response to heat stress. Immune stimulation could occur either due to the hyper thermic effect on immune cells or due to indirect effects such as the activation of heat shock factors, which are potent immune modulators and can stimulate both innate and adaptive immune response [86]. Heat stress triggered a massive innate immune response in the male group, with several genes including *CD14*, *BMX*, *S100A8*, *TLR4*, *MYD88*, *SRC*, *TLR2*, *SLEC5A*, *S100A9*, *PTK2B* and *VNN1* being up-regulated; interestingly in the female group, 13 genes involved in innate immunity were down-regulated; those included, *TRIM5*, *ISG20*, *RSAD2*, *DDX58*, *PML*, *TRIM25*, *TEC, IRF3*, *TRIM26*, *SH2D1A*, *ANXA1*, *FADD* and *JAK3.* In the male group, among the DEGs several inflammatory response genes were up-regulated; these included *CD14*, *CCR1*, *TSPAN2*, *TLR4*, *THBS1*, *PTGER3*,

*MYD88*, *C5AR2*, *SELP*, *TLR2* and *CHI3L1*, while inflammatory response was not enriched amongst the DEGs in the female group. Anti-tumor immune response were also markedly enhanced due to heat stress, genes involved in proteoglycans in cancer, choline metabolism in cancer, pathways in cancer were up-regulated in both male and female groups, while apoptotic genes were only triggered in the female group; these included genes *PIK3CA*, *BCL2L1*, *AKT3* and *ATM.* Several studies have shown that heat stress triggers anti-tumor and apoptotic pathways possibly due to protein aggregation due to heat stress induced denaturation [16].

#### *4.2. Metabolome Regulation in Response to Heat Stress*

Metabolites are building blocks for growth and development; they are also key regulators and markers of animal health [87]. Blood metabolites are highly sensitive to environmental stress [88] and HS alters protein metabolism in a number of species [19,89–91]. Metabolome analysis showed that several plasma metabolites were enriched in response to HS in both male and female groups, this included, creatinine, histidine, lysine, methionine, ornithine, serine, proline and pyruvate, while 4-aminobutryic acid, Creatine, Taurine and Urea concentration was lower post HS. Plasma creatinine concentration has been found to be increased due to heat stress in several species including cattle, pigs and sheep [92–94]. Creatinine along with methylhistidine is a known indicator of muscle breakdown [95], the concentration of 1-Methylhistidine was also significantly increased in the female group. This indicates that the animals might be experiencing tissue break-down due to heat stress, the reason for such breakdown is not clear, however it could probably be due to increased protein catabolism, which is required for utilizing the carbon in amino acids for gluconeogenesis [96]. Methionine is a sulfur containing amino acid, and along with creatinine is involved in cellular antioxidant mechanism and is found abundantly on the surface of proteins exposed to very high oxidant fluxes [97,98]. The concentration of Pyruvate, an anti-oxidant [99] has been found to be increased under hyperthermia in pigs [100]. Pyruvate is also the precursor to alanine via alanine aminotransferase, and the entry of pyruvate into the TCA cycle through the pyruvate dehydrogenase complex is impaired due to heat stress, possibly due to the inactivation of PDH (Pyruvate dehydrogenase complex) complex due to the HS induced oxidative chain reactions generated by intracellular ROS (reactive oxygen species) [101]. Along with increased plasma concentration of Pyruvate due to the inactivation of PDH, the increased expression of *PDK4* (Pyruvate Dehydrogenase Kinase 4) is suggested to serve as a mechanism to reduce substrate oxidation and mitochondrial ROS production to protect HS induced cellular damage [6]. The expression of *PDK4* was significantly increased in the male group. Glycine stimulates protein synthesis, and has been found to inhibit oxidative stress in pig small intestine, the plasma concentration of glycine in the female group increased post HS while it reduced in the male group [102]. Similarly Alanine and Citrulline have protective effect against oxidative damage [103]. The difference in the plasma concentration of several of the above discussed metabolites together with increased expression of *PDK4*, *HMOX1* and down-regulation of *SOD1* indicates that the male group could have experience heat stress induced oxidative stress.

Several plasma metabolites were also oppositely enriched in the male and female groups. The plasma concentration of Tryptophan, 1-Methylhistidine, Acetate, Glycine, Aminoadipic acid, Alanine, Arginine, Citrulline, Cystathionine, Glutamate, Threonine, and Valine were all increased in the female group, while their concentration was reduced in the male group. Tryptophan metabolites are key neurotransmitters, that regulate immune response, increased plasma concentration of tryptophan could indicate intracellular protein degradation, and were found to increase the expression of apoptosis initiators in pig [104,105]. The plasma concentration of metabolites like glutamate, cystathionine and threonine are all associated with apoptosis and autophagy [106–108]. Several anti and pro apoptotic genes were also amongst the DEGs. This along with previously known effects of HS on protein denaturation and the proteotoxic effect of HS [16], and the increased concentration of amino acids that are by products of protein degradation suggest that HS induced cellular damage must have occurred in female group.

#### **5. Conclusions**

In conclusion, heat stress triggered a dynamic response in pigs; however, the response to heat stress was sex-specific. The reason for this sexual dimorphic response is not clear; however, evidence from other studies suggests that it could be due to the anatomical, physiological and hormonal differences. Future transcriptome and metabolome studies, along with blood parameters and hormonal analyses could provide more understanding about the effect of sex on the heat stress response in pigs. The results of this study along with increasing our understanding of porcine heat stress response will serve as a good reference for future studies. However, there are some limitations in this study, such as only six animals being used in this study for both transcriptome and metabolome analyses, and since only NMR was used for metabolite quantification, it resulted in limited detection of metabolites (only 48 metabolites) resulting in only a limited understanding of the effect of HS on metabolite accumulation, so complementing NMR analysis with mass spectrometric analysis might result in a more comprehensive understanding of the effect of HS on metabolite accumulation.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/5/540/s1. File S1: List of primers used in the q-RT PCR analysis. File S2: Summary of the mapping and alignment statistics. File S3: List of DEGs identified in response to heat stress. File S4: List of metabolites found to be differentially enriched in response to heat stress. Figure S1: PCA using RNA-Seq and metabolome data from all samples. Table S1: Summary of the measurements of average daily feed and water intake.

**Author Contributions:** Conceptualization, K.S., S.D.L., G.-W.J. and J.-E.P.; methodology, K.S., S.Y.J., K.H.K., Y.K.L., M.K., Y.C.B. and S.D.L.; software, K.S.; validation, H.K. (Hana Kim); formal analysis, K.S., Y.C.B. and M.K.; investigation, K.S. and H.K. (Himansu Kumar); resources, J.-E.P. and S.D.L.; data curation, B.-H.C., G.-W.J. and S.Y.J.; writing—original draft preparation, K.S.; writing—review and editing, S.D.L. and J.-E.P.; visualization, K.S.; supervision, G.-W.J. and J.-E.P.; project administration, S.D.L.; funding acquisition, S.D.L. and J.-E.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was carried out with the support of the Cooperative Research Program for Agriculture Science and Technology Development (project number: PJ01277102), Rural Development Administration, Republic of Korea. K.S. and H.K.U. were supported by the 2020 RDA fellowship program of the National Institute of Animal Science, Rural Development Administration, Republic of Korea.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Review*
