**Prolonged E**ff**ect of Seminal Plasma on Global Gene Expression in Porcine Endometrium**

#### **Marek Bogacki \* , Beenu Moza Jalali, Anna Wieckowska and Monika M. Kaczmarek**

Institute of Animal Reproduction and Food Research of Polish Academy of Sciences in Olsztyn, 10-748 Olsztyn, Poland; beenu.jalali@pan.olsztyn.pl (B.M.J.); anna.kitewska@gmail.com (A.W.); m.kaczmarek@pan.olsztyn.pl (M.M.K.)

**\*** Correspondence: m.bogacki@pan.olsztyn.pl; Tel.: +48-89-5391-3131

Received: 8 October 2020; Accepted: 29 October 2020; Published: 3 November 2020

**Abstract:** Seminal plasma (SP) deposited in the porcine uterine tract at the time of mating is known to elicit an initial response that is beneficial for pregnancy outcome. However, whether SP has any long-term effect on alterations in endometrial molecular and cellular processes is not known. In this study, using microarray analyses, differential changes in endometrial transcriptome were evaluated after Day 6 of SP-infusion (6DPI) or Day 6 of pregnancy as compared to corresponding day of estrous cycle. Both, pregnancy and SP induced significant changes in the endometrial transcriptome and most of these changes were specific for a particular group. Functional analysis of differentially expressed genes (DEGs) using Ingenuity Pathway Analysis revealed that inhibition in immune response was affected by both pregnancy and SP infusion. Long-term effects of SP included differential expression of genes involved in inhibition of apoptosis, production of reactive oxygen species and steroid biosynthesis, and activation of processes such as proliferation of connective tissue cells and microvascular endothelial cells. Moreover, interleukin-2 and interferon-γ was identified to be responsible for regulating expression of many DEGs identified on 6DPI. The present study provides evidence for the long-term effects of SP on porcine endometrium that can be beneficial for pregnancy success.

**Keywords:** seminal plasma; endometrium; global gene expression; microarray; pig

#### **1. Introduction**

The high rate of pregnancy failure in human and livestock has been attributed mainly to the unsynchronized development of the embryos with the proper preparation of the female reproductive tract and the impaired communication between the developing embryos and uterus [1–3].

Understanding of the molecular embryo-maternal cross-talk is crucial for solving infertility problems, reducing pregnancy loss and identifying hormonal, paracrine, and autocrine factors regulating the developmental potential of the offspring. Effective recognition of the embryo in the maternal tract is crucial for the preparation of an appropriate environment in the uterus for the embryo's development, implantation, and final establishment of pregnancy [4]. However, exactly when the oviduct and uterus recognize the presence of embryos and how the maternal pathway changes its environment in response to embryos is not fully understood.

In pigs, transcriptomic profiling of pregnant and non-pregnant animals has been conducted and pointed major differences in endometrial genes activities in the post-conception period of pregnancy [5], pre-attachment phase [6,7] and during onset of implantation [8]. Identified alterations in uterine transcriptome lead to morphological, biochemical and immunological changes and are reflection of action of para- and autocrine signals released by maternal tract as well as developing embryos.

The application of artificial insemination (AI) and other reproductive technologies shows that pregnancy can be maintained without any semen being deposited in the uterus (embryo transfer outcome) or with highly diluted semen for AI [9]. However, studies conducted in pigs and other livestock species show reduction of early fetal loss due to exposure to seminal constituents. Seminal plasma (SP) is a mixture of various components and serves not only as a vehicle for sperm transport, protection, and nutrition but also affects gamete interaction and successful establishment of pregnancy. Biologically active molecules present in SP (estrogens, testosterone, prostaglandins, cytokines, and growth factors) interact with uterine epithelium to induce a series of changes in the maternal immune response and modify cellular composition, structure, and function of the reproductive tract, directly in local tissues or indirectly in tissues distal to the uterus, for example ovary [10]. In pigs, it was shown that introduction of SP before natural service or/and AI leads to the increase in farrowing rate and litter size [11]. Additionally, increased litter size was also reported after pre-sensitization of the uterus to sperm and seminal antigens in a previous estrous cycle [12] and an increased embryo survival was noted after addiction of leukocyte antigens to boar semen at breeding [13].

The immediate response to full semen insemination in pigs is a rapid influx of inflammatory cells such as neutrophils into the uterine lumen [14] and macrophages, granulocytes, and lymphocytes into the endometrial stroma [15]. The activation of maternal immune system does not cause the rejection of seminal antigens due to the presence of several immunoregulatory molecules in boar SP, such as prostaglandin E (PGE) and tumor growth factor β (TGF-β) [16]. It has been suggested that constituents of SP deposited in lower reproductive tract can easily access the upper reproductive tract and induce biologically relevant changes in the endometrium [17]. Concomitant with this observation, SP interacts with uterine cells in pigs to induce expression of several cytokines: granulocyte-macrophage colony-stimulating factor (GM-CSF), interleukin 6 (IL-6), and macrophage chemotactic protein-1 (MCP1) 34 h after infusion [18]. Moreover, components of SP can alter expression of prostaglandin synthesis enzymes in the porcine oviduct and uterus, acting in favor of PGE2 action as a critical element of early embryo transport and development but also can modify endometrial vascularity up to 10 days after infusion [19,20]. Additionally, in our previous study, prolonged modulatory effects of SP infusion at least for 6 days, demonstrated by induction of T helper (Th) and T regulatory (Treg) cells, increased interleukin 10 (IL10), and decreased expression granulocyte-macrophage colony stimulating factor (GM-CSF), was observed [21].

Though, the immune modulatory effects of SP on porcine endometrium are well documented, to our knowledge there is no published microarray data showing long term effect of SP on molecular changes in porcine endometrium that might be relevant for hatching blastocysts on Day 6 of pregnancy when they reach the uterine horn from the oviduct. That is why, we hypothesized that SP infusion can induce not only immediate but also prolonged transcriptome changes required for endometrial receptivity and modulation of immune response in the uterine environment.

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

#### *2.1. Animals and Treatments*

All procedures involving the use of animals were conducted in accordance with the national guidelines for agricultural animal care and approved by the Animal Ethics Committee, University of Warmia and Mazury in Olsztyn, Poland; Decision 86/2011. Estrous induction and synchronization, insemination, and SP infusion were performed as previously described [21]. Cross-bred gilts (*Sus scrofa domesticus*) of similar age and weight were subjected to a hormonal treatment with an intramuscular injection of 750 I.U equine chorionic gonadotropin (eCG) and 500 I.U. human chorionic gonadotropin (hCG) given after 72 h. Gilts randomly divided into three groups (*n* = 5), 24 h after hCG injection were treated as following: (1) artificially inseminated twice within an interval of 12 h, (2) received two intrauterine infusions of 100 mL SP with an interval of 12 h or (3) received two intrauterine infusions of 100 mL PBS within 12 h (reference group). All the treatments were given at two time points

within an interval of 12 h to mimic regular procedure of artificial insemination. Artificial insemination was performed with 100 mL of 2.5 × 109 spermatozoa diluted from the fresh semen collected from a boar with semen extender Safe Cell + (IMV technologies, L'Aigle, France). All the treatments were deposited into the uterus using post cervical artificial insemination methods (PC Blue, SafeBlue Foamtip® with PC Cannula-Minitube) to facilitate interaction of treatments with the endometrium. SP for intrauterine infusion was prepared from whole semen collected and pooled together from four fertile boars, which were used for AI. SP was separated by centrifugation of the whole semen at 1200× *g* at 4 ◦C for 20 min, divided into 100 mL aliquots and frozen at –20 ◦C until needed for intrauterine infusion. Animals were slaughtered in a local abattoir at day 6 of pregnancy or at day 6 after SP or PBS infusion. Uterine horns were flushed with PBS and opened longitudinally along the anti-mesometrial surface. Endometrial explants were collected from the upper part of uterine horns and snap-frozen in liquid nitrogen. In the group of artificially inseminated animals, only those animals were included in the study in which pregnancy was confirmed by the presence of blastocysts in the uterine flushings.

#### *2.2. RNA Isolation*

Total RNA was isolated from 30 mg of grinded in liquid nitrogen and homogenized endometrial tissue with the use of a Qiagen RNeasy Mini Kit (Qiagen, Valencia, CA, USA) and genomic DNA contamination was removed by DNAse treatment (RNase free DNAse Kit, Qiagen, Valencia, CA, USA). The initial RNA quality and quantity were determined spectrophotometry using NanoDrop ND-1000 (Thermo Scientific, Pittsburgh, PA, USA). Subsequently, RNA integrity was evaluated with microfluidic electrophoresis by 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and RNA integrity number (RIN) was calculated for each sample using Agilent 2100 Expert software (Agilent Technologies, Inc., Santa Clara, CA, USA). Only samples with a RIN above 8.0 were processed further.

#### *2.3. Microarrays*

The Porcine (V2) Gene Expression Microarray 4\_44 (Agilent Technologies, Santa Clara, CA, USA) was used for differential gene expression analysis. As positive internal controls of microarrays performance an RNA Agilent Spike-In Kit, One Color was used (Agilent Technologies, Santa Clara, CA, USA). Total RNA obtained from cyclic and pregnant gilts was amplified and labeled with Cy3 dye using Quick Amp Kit, One Color (Agilent Technologies, Santa Clara, CA, USA). After purification of labeled RNA (Qiagen RNeasy Kit), the yield (ng of cRNA) and specific activity (pmol of Cy3/mg of cRNA) were quantified using NanoDrop ND-1000. Labeled cRNA was fragmented, mixed with hybridization buffer, and added to the microarray slide. On each array (*n* = 4, one slide) a combination of samples from all three groups were hybridized for 17 h at 65 ◦C in an Agilent hybridization oven. Afterwards, arrays were dissociated from the hybridization chamber and washed twice in GE wash buffers. After washing, slides were scanned using Agilent G2565CA Microarray Scanner at settings recommended for the 4\_44 K array format. Images obtained after scanning were analyzed using Agilent Feature Extraction software v. 10.5.1.1 (Agilent Technologies Inc., Santa Clara, CA, USA). A detailed analysis including filtering of outlier spots, background subtraction from features, and dye normalization was performed.

#### *2.4. Data Analysis*

Data obtained after extraction was further analyzed using GeneSpring GX 11.0.2 (Agilent Technologies, Santa Clara, CA, USA). To determine differentially expressed genes (DEGs) data were normalized with quantile method and afterwards moderated *t*-test (Benjamini–Hochberg false discovery rate (FDR) < 0.05, absolute fold change |Fc| > 1.5) was performed to compare endometrial transcriptomes between: (1) pregnant and cyclic (*n* = 5) as well as (2) SP infused and cyclic animals (*n* = 4, data from one array were not correlating with other arrays after principal component analysis (PCA)). For identification of differentially expressed probe sets with unknown target sequence the annotation was done manually using NCBI blast algorithm [22]. When porcine sequence for particular mRNA was not available the annotation was performed for human, murine, and cattle transcript, and the date was included only if the query cover and percent identity was equal or higher than 70%. To identify biological processes, pathways and upstream regulators Ingenuity Pathway Analysis (IP Ingenuity Systems-Qiagen, Aarhus, Denmark). Biological processes, pathways, and upstream regulators were considered statistically significant if Fishers exact tests *p*-value ≥ 0.05 and each process associated with at least four DEGs. Biological processes and pathways connected to cancer, diseases and disorders, nervous, respiratory, digestive, renal and urological system, organismal survival and functions, drug metabolism, organ development and behavior were not taken into consideration while examining IPA results.

#### *2.5. Quantitative Real-Time PCR*

For validation of microarray results 2 µg of total RNA were transcribed to cDNA with the use of High Capacity RNA to cDNA kit (Applied Biosystems, Foster City, CA, USA). Real-time PCRs were performed using 7900 HT Real-Time PCR System (Applied Biosystems) using 15 ng cDNA, TaqMan Universal MasterMix II, no UNG. TaqMan assays are listed in Table 1. The initial denaturation was carried out at 95 ◦C


**Table 1.** TaqMan assays used for real-time PCR validation of microarray results.

\* Reference genes.

For 15 min and was followed by 40 cycles of denaturation at 95 ◦C (15 s) and primer annealing and elongation at 60 ◦C. Non template controls and non-enzyme controls were included in the experiment. Gene expression levels were calculated with the use of Real-Time PCR Miner (http: //ewindup.info/miner/) and normalized using the geometric mean of expression levels of two reference genes-hypoxanthine guanine phosphoribosyl transferase (HPRT) and β-actin (ACTB). The statistical differences in gene expression between the endometrium from pregnant and cyclic as well as seminal plasma infused and cyclic animals was analyzed with GraphPad PRISM v. 5.0 Software (GraphPad Software, Inc., San Diego, CA, USA) by Student's *t*-test. Confirmed differences in gene expression (*p* < 0.05) were expressed as fold changes.

#### **3. Results**

#### *3.1. Di*ff*erential Changes in Endometrial Transcriptome*

Pairwise comparisons of endometrial samples collected from pigs on Day 6 of pregnancy (6DP) and Day 6 after SP infusion (6DPI) with PBS-infused cycling control pigs on corresponding day of estrous cycle (6DC) were performed to identify endometrial transcriptome changes in response to pregnancy and SP constituents, respectively. Statistical analysis revealed a pregnancy-induced change (fold change > 1.5 or < −1.5; *p* < 0.05; false discovery rate (FDR) = 5%) in 444 probes representing 281 differentially regulated (DEGs) of which 225 were downregulated and 56 were upregulated, on Day 6

as compared to controls on 6DC (Table S1). Whereas, genes showing the most downregulation included *S1008*, *CGA*, and *HLA-DQA1*, genes with highest upregulation were *GPR116*, *COL4A1*, and *CADPS2*. On the other hand, SP-infusion resulted in statistically significant alterations in 342 probes representing 255 DEGs. A downregulation of 118 genes and upregulation of 137 genes was observed on Day 6 after SP infusion as compared to 6DC (Table S2). SP infusion decreased the expression of *ATP6V1C2, NMU, S100A8, S100A12, ANGPTL3, NOS1, CCR3* and upregulated *PCDHB15*, *KLHL5*, *RASGEF1A, NMB*, and *CAPZB*.

#### *3.2. Comparison between Pregnancy-Induced and Seminal Plasma-Induced Changes in Endometrial Transcriptome*

The list of pregnancy and SP-induced DEGs were uploaded to jvenn software (http://jvenn. toulouse.inra.fr/app/example.html) to visualize common DEGs across both presented comparisons and DEGs that were identified as a result of either pregnancy or SP-infusion (Figure 1). This comparison identified 19 common genes (Table S3) that were differentially regulated by pregnancy and SP-infusion. Whereas, in both the groups, 15 genes were found to be downregulated only three genes were upregulated. The expression of transcript coding for chloride channel, voltage sensitive 5 (CLCN5) was lower during Day 6 of pregnancy (Fc = −2.66) and higher on Day 6 after SP infusion (Fc = 1.53). Most of the DEGs common between two groups were involved with immune regulation (*IL15, IL18, LGALS1, FKBP3*) or were DEGs related with molecular transport (*S100A8, S100A12, CLCN5*) and structural organization (*FN1, COL7A1, TUBA1B*). −

**Figure 1.** Venn diagram showing number of shared and unique DEGs altered on Day 6 of pregnancy (6DP) and Day 6 after SP infusion (6DPI) in comparison to 6 Day of estrous cycle. Genes chosen for qPCR validation are listed in boxes.

#### *3.3. Analysis of Biological Processes, Pathways, and Upstream Regulators of Identified DEGs*

To classify identified DEGs altered on Day 6 of pregnancy and Day 6 after SP infusion under functional categories, tools available in IPA database were employed. Analysis of DEGs using core analysis module of IPA revealed many altered functions in porcine endometrium as a result of pregnancy or SP-infusion. Furthermore, a comparison analysis, comparing disease and biological functions, canonical pathways, and upstream regulators of DEGs was also performed to evaluate the differences in activated and inhibited functions between pregnant and SP-infused animals (Figure 2). Whereas, most of

the identified functions were specific to either pregnancy status or SP-infusion, inhibition of immune functions was common to both the groups. A pregnancy specific activation was observed in processes related to organization of cytoskeleton, organization of cytoplasm, and transmigration of leukocytes (*Z* > 2.0; Figure 3A). The molecular functions inhibited by pregnancy included chemoattraction of leukocytes, homing of leukocytes, and cytotoxicity of lymphocytes (*Z* < −2.0) (Figure 3B, Table S4). On the other hand, SP infusion, besides inhibiting cytotoxicity of lymphocytes, also inhibited molecular functions such as apoptosis, lipid metabolism, senescence of fibroblasts, and production of reactive oxygen species (*Z* < −2.0; Figure 3C, Table S5) and activated only function related to proliferation of microvascular endothelial cells and connective tissue cells (*Z* > 2.0; Figure 3D). Interestingly, alterations in molecular signaling pathways associated HIF1α signaling were specific to Day 6 of pregnancy and activation of Wnt/β-catenin signaling was specific to SP-infused group. Genes associated with PPAR signaling were differentially altered in both the groups (Figure 2B).

**Figure 2.** Comparison analysis of (**A**) diseases and bio functions, (**B**) canonical pathways, and (**C**) upstream regulators associated with DEGs identified on 6DP and 6DPI. A *Z*-score > 2 (orange color) is associated with activated functions, pathways, or upstream regulators and a *Z*-score < 2 (blue color) is associated with inhibited functions, pathways, or upstream regulators.

**Figure 3.** IPA depicting networks integrating DEGs identified on 6DP as compared to 6DC in (**A**) activated functions: organization of cytoskeleton, organization of cytoplasm, and transmigration of leukocytes, and (**B**) inhibited immunological functions: cytotoxicity of lymphocytes, chemoattraction of leukocytes, and homing of leukocytes. DEGs identified that 6DPI as compared to 6DC were found to be associated with (**C**) activated functions: proliferation of connective tissue cells and proliferation of microvascular endothelial cells and (**D**) inhibited functions: cytotoxicity of lymphocytes, apoptosis, production of reactive oxygen species, concentration and release of lipid and quantity of steroid, on Day 6 after SP infusion in the porcine endometrium. Red and green colors depict an increase or decrease, respectively, in the expression of a gene. The color intensity of nodes indicates a fold change increase or decrease associated with a particular DEG.

In our study, upstream analysis function of IPA was used to identify the molecules including cytokines, transcription factors, or hormones (2.0 < *Z* score < −2.0) that might be upstream regulators of altered gene expression in porcine endometrium as a result of pregnancy or SP-infusion. Many cytokines including interleukin (IL-1β), IL-2, TNF, transcription regulators such as FOXO1, NUPR1 and nuclear receptor ESR2, PPARA, and AHR were found to affect the gene expression in endometrium on Day 6 of pregnancy (Figure 4A, Table S6). Upstream analysis of SP-induced DEGs revealed that upstream regulators such as IL-2 and RICTOR were the only common regulators among two groups. Other molecules affecting the expression of genes in SP-infused endometrium on Day 6 included *IFN*γ*, GFI1, HSF1, TNFRSF1A, TLR2, and NR1H3* (Figure 4B, Table S7). Interestingly, we observed PGR to be an upstream regulator of SP-induced DEGs (Figure 4B).

**Figure 4.** Networks generated for identified upstream regulators of DEGs in porcine endometrium on (**A**) Day 6 of pregnancy or (**B**) Day 6 after SP infusion. Red and green colors depict an increase or decrease, respectively, in the expression of a gene. The color intensity of nodes indicates a fold change increase or decrease associated with a particular DEG.

#### *3.4. qRT-PCR Validation of Microarray Results*

For q PCR validation of microarray data, 10 DEGs, shown in Figure 1, were chosen. These genes were associated with immune function, molecular transport, and cell proliferation. Most of the assessed DEGs showed similar expression profiles when comparing microarray and qPCR data (Figure 5). However, qPCR data revealed a pregnancy induced upregulation of TGFA expression that was not observed in the microarray data. A comparison of fold change and *p* value of DEGs obtained after qPCR and microarray data analysis is presented in Table 2.


**Table 2.** Results of microarray experiment validation with qPCR. 6DP—6 Day of pregnancy vs. 6 Day of estrous cycle, 6DPI—6 Day after SP infusion vs. 6 Day of estrous cycle, Fc—fold change.

**Figure 5.** *Cont.*

**Figure 5.** *Cont.*

**Figure 5.** Validation of microarray results using qPCR. Expression of *CCR3* (*chemokine* (*C-C motif*) *receptor 3*), *CXCL11* (*chemokine* (*C-X-C motif*) *ligand 11*), *TGFA* (*transforming growth factor,* α), *LGALS1* (*galectin 1*), *IL-18* (*interleukin 18*), *LY96* (*lymphocyte antigen 96*), *PDCD10* (*programmed cell death 10*), *SLA-DQA1* (*MHC class II histocompatibility antigen SLA-DQA*), *S100A8* (*S100 calcium binding protein A8*), and *S100A12* (*S100 calcium binding protein A12*) in the Day 6 of pregnant or SP infused animals as compared to Day 6 of cycling control animals. Expression values were normalized to expression of *ACTB (Actin,* β*)* and *HPRT (Hypoxanthine phosphoribosyltransferase1)*. Data are presented as mean ± standard error. \* *p* ≤ 0.05, \*\* *p* ≤ 0.01, \*\*\* *p* ≤ 0.001 by *t*-test. The fold change and *p* values are shown in Table 2.

#### **4. Discussion**

Although, nowadays, pregnancy in pigs is a result of AI with diluted semen or the result of embryo transfer techniques during which only the residual SP enters the reproductive tract, there is documented evidence that SP affects the biological functions of the uterus and evidence that interaction between male SP and female tissues promotes fertility, pregnancy, and finally health of offspring [23]. Many transcriptomic studies in humans, cattle, and pigs have been carried out to evaluate the effects of SP on endometrium [17,24,25]. These reports support the results that SP itself modifies the transcriptome, although semen either after mating or AI results in the maximum changes in the molecular profiles in peri-ovulatory uterine tract of pigs [25,26]. However, in pigs, these studies either reported the immediate effects (after 24 h) of SP on uterine tract or effect of SP followed by AI [27]. Our previous studies have shown that SP can induce long term changes in the endometrium that can be observed at least till Day 6 of its infusion, therefore, in this study, we evaluated SP-induced long-term changes in endometrial transcriptome to identify significantly altered molecular and cellular processes that might prepare endometrium for a possible pregnancy. We also compared these changes with the list of DEGs obtained on Day 6 of pregnancy to evaluate distinct and shared pathways between the two treatments.

Our data demonstrated that as many as 255 and 281 genes are differentially regulated after 6 days of SP infusion and on Day 6 of pregnancy as compared to corresponding day of estrous cycle with only 19 being common to both the groups. A comparison of the biological, molecular, and cellular functions altered by SP-infusion or pregnancy revealed that most of these processes are specific to either SP-infused or pregnant groups of animals, highlighting specific actions of SP constituents. Many DEGs found in both the groups were responsible for inhibition of immune function. Processes such as

organization of cytoskeleton and transmigration of leukocytes were specific for pregnancy induced DEGs. Treatment with SP inhibited processes such as apoptosis, necrosis, production of reactive oxygen species and steroid transport. On the other hand, connective tissue cell and microvascular endothelial cell proliferation was activated by SP. Whereas, pathways affected by SP, such as endometrial immune response and steroid biosynthesis were inhibited after Day 6 of its infusion, these responses were activated immediately after SP infusion [25].

#### *4.1. Immune Regulation*

Consistent with the literature reports, modulation of immune responses was one of the topmost processes altered by both the treatments, i.e., AI and SP. Blastocysts and SP are known to differentially regulate genes involved in immune response on Day 6 [5,21]. It is well known that immediate effects of SP on endometrium include an inflammatory reaction, a response to paternal antigens and mainly to clear reproductive tract of any pathogen deposited at the time of mating [13]. A recent study reported minimal effect of SP treatments on the differential expression of genes in the porcine upper reproductive tract [25]. Consistent with these reports, a very small number of immune related genes were differentially regulated in our study. However, there was no difference between the genes regulated either as a result of pregnancy after AI (six DEGs) or SP treatment (eight DEGs) showing the similarities between both the treatments. In the present study genes involved with immune functioning such as *IL15*, *IL18*, and *LGALS1* were found to be downregulated in both the groups and additionally, *STAT5* and *GZMA* was downregulated by SP infusion. Interleukins 15 and 18 are pro-inflammatory cytokines, these cytokines are not only the regulators of innate immune response but also enhance the cytotoxicity of natural killer cells (NK cells) [28,29]. Furthermore, a decrease in cytotoxicity of lymphocytes was also evident from the downregulation of *GZMA*, a factor secreted by the cytotoxic T cells and natural killer cells that induces apoptosis [30]. A downregulation of these factors might result in a dampened innate immune response at the time of blastocyst hatching and in turn may affect endometrial immune tolerance to paternal antigens. Downregulation of many of the genes associated with immune regulation in SP-infused animals was found to be a result of inhibition of toll-like receptor (TLR) 2 and IL-2 signaling. The activation of TLR-signaling is indeed found to be detrimental to the success of pregnancy [5]. In the endometrium of pregnant pigs, negative regulation of immune responses could be a result of inhibition TNF or IL-1 signaling (*Z* > –2.0). Both SP and pregnancy status are able to generate moderate changes in endometrial immune response. Current data closely corresponded with our previous results suggesting that the immunomodulatory effects of SP, also observed at the protein level, last up to at least 6 days after its infusion at which time blastocysts enter the uterus from the oviduct [21,23].

#### *4.2. Cell Death and Survival*

A large number of endometrial genes associated with the category cell death and survival were downregulated in SP-infused animals as compared to PBS-infused controls. These DEGs resulted in an inhibition of the apoptotic signaling and promotion of microvascular and connective tissue cell proliferation. Corresponding to the proliferation of connective tissue, an inhibition in senescence of fibroblasts was also observed, confirming an earlier observation that SP constituents can also have an effect on stromal layer of endometrium [17]. Our data emphasize the utility of SP in suppressing apoptosis in endometrial cells during early pregnancy. Increase in proliferation and inhibition of apoptosis of endometrial cells during early pregnancy period is an important step for generation of receptivity to implanting embryos in pigs [31]. A recent report reveals pro-survival effect of SP on endometrial epithelial and stromal cells [17]. In vitro treatment of human endometrial cells resulted in inhibition of pathways promoting apoptosis [17]. In this study we observed downregulation of many apoptotic genes such as *TNFAIP3, CXCL11, ABCC4,* and *B2M* and upregulation of genes inhibiting apoptosis including *FN1, GLRX1, CDK2AP1*, and *LGALS3*. Though direct participation of all of these genes has not been evaluated in endometrial cell proliferation or in apoptosis, their role in apoptosis

is established. The tumor necrosis factor-α-induced-protein 3 (TNFAIP3) and chemokines including CXCL11 are inducers of epithelial apoptosis [32,33] which is important for species with invasive implantation, however, in pigs, conceptuses do not breach the epithelial layer during implantation. An inhibition in endometrial apoptosis resulting from downregulation of these genes will be favorable for generation of receptive endometrium in pigs.

Many genes that were upregulated in the endometrium of SP-infused pigs were found to participate in proliferation of microvascular endothelial cells and proliferation of connective tissue cells that consists mostly of fibroblast. The genes associated with these categories included *TGFA*, *MAPK1, ADM, LGALS3, DNMT3B* and *DKK*. Galectins, including LGALS3, are multifunctional proteins associated with immune regulation, cell proliferation, and angiogenesis. Growth factor, TGFA, that activates epidermal growth factor receptor (EGFR) and LGALS3 has been shown to mediate angiogenesis through VEGF and bFGF-mediated angiogenic response [34–36]. Increased endometrial angiogenesis is a hallmark of successful pregnancy, it ensures proper growth and development of the embryo. A SP-induced upregulation of these factors suggests its possible effect on endometrial vasculogenesis.

#### *4.3. Oxidation Stress*

A downregulation of genes associated with production of reactive oxygen species (ROS) such as *MAOA, NCF2, NOS3* and *CCR3* was observed in endometrial transcriptome of SP-infused pigs. A balance in production ROS is important as at moderate concentrations it has important signaling roles under physiological conditions, but sustained ROS production can have detrimental effects. The expression of *MAOA* which is induced by TNFα has been reported in human and rat endometrium where it is upregulated during the implantation period. The role of this gene in endometrium is not clear yet. However, as *MAOA* is identified as a source of ROS generation, we speculate that a decreased activity of this gene during early pregnancy might create a decrease in detrimental ROS species which in turn might act on downregulation of *CCR3* expression. We also observed an increase in the expression of gene coding antioxidative enzyme, *GLRX*. Glutaredoxins were reported to be expressed in human endometrium where they are associated with the regulation of the cellular redox state and antioxidation defense mechanisms [37].

#### *4.4. Lipid Metabolism*

Processes related to lipid metabolism, such as steroid and lipid quantity and lipid release were found to be inhibited in the endometrium of SP-infused animals (Table S5). Some of the DEGs related to this process included downregulated expression of STAR, NR5A1, and Cyp7A1. Steroid biosynthesis has previously been reported to be regulated by the SP-infusion. However, the process was shown to be activated after 24 h of the SP treatment [25]. In our study, the one probable explanation for downregulation of these genes could be SP-induced upregulation of COX-2 [18]. Though, regulation of steroid biosynthesis by COX-2 has not been evaluated in porcine endometrium, it was shown to be a negative modulator of steroid biosynthesis in Leydig and bovine luteal cells [38,39]. More importantly, it has been suggested that COX-2 inhibits steroidogenesis through MAPK-signaling [40] and MAPK was found to be upregulated after SP-infusion in our study. As we did not measure systemic progesterone concentrations, whether a result of downregulation of these genes was observed in systemic progesterone concentration needs to be evaluated. Interestingly, it has been reported that early progesterone treatment decreases uterine capacity at 105 days of gestation due to accelerated fetal growth [41]. Moreover, a slower rate of conceptus development is attributed to greater fertility as in Meishan breed of pigs. Consistent with these findings, SP has also been shown to increase embryo viability and at the same time slow the growth (size) of embryos on day 9 of gestation [18]. Our observations concerning a possible long-term effect of SP on steroid biosynthesis are worth further exploration in terms of effect of SP on fertility rate in *Sus scrofa domesticus* through moderate changes in progesterone synthesis.

In conclusion, our results clearly show that SP can induce long term effects on the gene expression in the porcine endometrium. Long-term effects of SP on endometrium include inhibition of processes related to immune response, apoptosis, and steroid biosynthesis and activation of processes such as proliferation of cells. In modern pig breeding, use of SP has been neglected. Our study paves the way for further research on the benefits of addition of SP or its constituents to the semen extenders during AI. The effects of SP on endometrium might prove to be advantageous for blastocyst development, preparing uterus for the conceptus attachment and finally for improving the fertility rate in pigs.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/11/1302/s1, Table S1: Microarray analysis data listing differentially regulated genes that were up- or downregulated on day 6 of pregnancy as compared to day 6 of the estrous cycle, Table S2: Microarray analysis data listing differentially regulated genes that were up- or downregulated on day 6 after seminal plasma infusion as compared to day 6 of the estrous cycle, Table S3: List of 19 common differentially regulated genes between day 6 of pregnancy and day 6 after seminal plasma infusion, Table S4: Disease and function analysis (IPA analysis) of DEGs identified on day 6 of pregnancy as compared to day 6 of the estrous cycle, Table S5: Disease and function analysis (IPA analysis) of DEGs identified on day 6 after seminal plasma infusion as compared to day 6 of the estrous cycle, Table S6: List of upstream regulators that might regulate the expression of DEGs identified on day 6 of pregnancy, Table S7: List of upstream regulators that might regulate the expression of DEGs identified on Tabelday 6 after seminal plasma infusion.

**Author Contributions:** Conceptualization, M.B.; methodology, A.W., M.B. and M.M.K.; software, A.W., M.M.K.; validation, A.W.; investigation, A.W., M.B.; resources, M.B.; data curation, B.M.J.; writing—original draft preparation, A.W., M.B., B.M.J.; supervision, M.B.; project administration, M.B.; funding acquisition, M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Science Center Grant No 2011/01/B/NZ4/03603.

**Acknowledgments:** We would like to thank Michal Blitek and Piotr Grezlikowski for helping with the animal preparation and handling as well as Paula Wojnicz and Marta Romaniewicz for their technical assistance.

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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*

## **Expression Profile of Porcine TRIM26 and Its Inhibitory E**ff**ect on Interferon-**β **Production and Antiviral Response**

**Hui Huang 1,2,**† **, Mona Sharma 1,**† **, Yanbing Zhang 1 , Chenxi Li 1 , Ke Liu 1 , Jianchao Wei 1 , Donghua Shao 1 , Beibei Li 1 , Zhiyong Ma 1 , Ruibing Cao 2, \* and Yafeng Qiu 1, \***


Received: 29 August 2020; Accepted: 15 October 2020; Published: 19 October 2020

**Abstract:** TRIM26, a member of the tripartite motif (TRIM) family has been shown to be involved in modulation of innate antiviral response. However, the functional characteristics of porcine TRIM26 (porTRIM26) are unclear. In this study, we used a synthesized antigen peptide to generate a polyclonal antibody against porTRIM26 with which to study the expression and function of porTRIM26. We demonstrated that polyinosinic:polycytidylic acid (poly (I:C)) stimulation and viral infection (vesicular stomatitis (VSV) or porcine reproductive and respiratory syndrome virus (PRRSV)) induce expression of porTRIM26, whereas knock-down expression of porTRIM26 promotes interferon (IFN)-β production after poly (I:C) stimulation and virus infection (VSV or PRRSV). The importance of the porTRIM26-mediated modulation of the antiviral response was also shown in VSV- or PRRSV-infected cells. In summary, these findings show that porTRIM26 has an inhibitory role in IFN-β expression and the antiviral response.

**Keywords:** TRIM26; antiviral response; IFN-β; poly (I:C); VSV; PRRSV

#### **1. Introduction**

Tripartite motif (TRIM) proteins, a large family of ubiquitin E3 ligase, include more than 80 and 60 members in human and mouse, respectively [1]. Moreover, more than 50 members of porcine *TRIM* genes have been annotated in the GenBank database, although most of them were predicted with computational analyses [2,3]. Recent studies have shown that many members of the TRIM family are expressed in response to interferons (IFNs) and are involved in the processes of innate immune response, especially during viral infections [4–6]. These reports strongly suggest that understanding the molecular functions of porcine TRIM proteins could offer insights into the regulation of innate immune response in swine species.

TRIM26 is a member of the TRIM family, and has a structure similar to that of many other members of the family. They are all structurally characterized by a RING finger domain (E3 ligase with ubiquitin ligase activity) with two B-box domains, followed by a coiled-coil (CC) region and a C-terminal protein binding domain [7–9]. In human, *TRIM26* gene is located in the MHC class I region [10,11]. Likewise, the porcine MHC (swine leukocyte antigen (SLA)) region contains *TRIM26* gene according to a sequencing analysis. However, although porcine *TRIM26* (*porTRIM26*) gene has been identified [2,12,13], the biological functions of the protein, especially in the immune response, remain to be determined.

Previous studies in human and mouse cell lines have shown that TRIM26 plays a controversial role in the regulation of IFN-β production and innate antiviral response, which may be attributable to the different experimental systems used [14,15]. However, how porTRIM26 affects the regulation of IFN-β expression and the innate antiviral response is unknown. In this study, for the first time, we determined the expression profile of porTRIM26 in different pig tissues and clarified its effect on the modulation of IFN-β expression and innate antiviral response.

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

#### *2.1. Tissue Sample Collection*

Piglets (Shanghai great white pig strain; ~30 days old) were purchased from the Shanghai Academy of Agricultural Sciences (Shanghai, China), euthanized, and dissected to obtain various tissue samples (liver, spleen, lung, kidney, submandibular lymph node, hilar lymph node, mesenteric lymph node, inguinal lymph node, and thymus). The tissue samples were immediately snap-frozen in liquid nitrogen and stored at −80 ◦C. All animal experiments were approved by the Institutional Animal Care and Use Committee of Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Shanghai, China (IACUC No: Shvri-po-2016060501) and were performed in compliance with the Guidelines on the Humane Treatment of Laboratory Animals (Ministry of Science and Technology of the People's Republic of China, Policy No. 2006398).

#### *2.2. Cells, Viruses and Infections*

Porcine alveolar macrophages (PAM) were generated as shown in our previous study [16]. HEK 293T cells, porcine iliac artery endothelial cells (PIEC), ST cells, PK-15 cells, BHK-21 cells, and Marc-145 cells were maintained in Dulbecco modified Eagle medium (Thermo Fisher Scientific, Shanghai, China) supplemented with 10% fetal bovine serum (Thermo Fisher Scientific, Shanghai, China) at 37 ◦C in a 5% CO<sup>2</sup> atmosphere. Vesicular stomatitis virus (VSV) and highly pathogenic porcine reproductive and respiratory syndrome virus (HP-PRRSV) were maintained in our lab and were propagated on BHK-21 cells and Marc-145 cells, respectively. The viruses were titrated with the median tissue culture infective dose (TCID50) methods, as described in a previous study [17]. According to the progress of virus infection and expression of porTRIM26, we chose 24 h after infection for collecting the viral infected samples: PIEC cells were infected with VSV at a multiplicity of infection (MOI) of 1 for 24 h, at which point peak titer was reached with the induction of porTRIM26; PAM were infected with PRRSV at a multiplicity of infection (MOI) of 1 for 24 h, at which point peak titer was reached with the induction of porTRIM26.

#### *2.3. Cloning and Sequence Analysis of porTRIM26*

The primers for cloning *porTRIM26* gene (shown in Table S1) were designed based on the *porTRIM26* gene sequence (GenBank accession number, NM\_001123209.1). The amplified sequence was confirmed with DNA sequencing and inserted into the p3×Flag-CMV-14 vector (Sigma, St. Louis, MO, USA) to generate a recombinant plasmid expressing FLAG-tagged porTRIM26 (pFlag-porTRIM26). An amino acid sequence alignment of the deduced protein sequence by this construct and the TRIM26 proteins of another three species, *Homo sapiens* (NP\_001229712.1), *Mus musculus* (NP\_001020770.2), and *Rattus norvegicus* (XP\_008770914.1), was conducted using the software Lasergene version 7.1 (Madison, WI, USA). A phylogenetic tree based on the sequences of different species was constructed by the neighbor-joining method using the MEGA software (version 6.06).

#### *2.4. Generation of Polyclonal Antibody Against porTRIM26*

A polyclonal antibody directed against porTRIM26 was generated as described in a previous study [18]. Briefly, a peptide of 20 amino acids corresponding to residues 364–383 of the porTRIM26 sequence was synthesized chemically and conjugated with keyhole limpet hemocyanin (KLH) as the carrier protein. Rabbits were immunized five times with the peptide-KLH conjugate combined with complete or incomplete Freund's adjuvants. HEK 293T cells were then transfected with the porTRIM26-expressing recombinant plasmid described in Section 2.3 and a Western blotting analysis was performed to confirm the specificity of the polyclonal antibody. All of the animal experiments were approved by IACUC in Shanghai Veterinary Research Institute, CAAS (No: Shvri-po-201606 0501) and followed the guidelines described in Section 2.1.

#### *2.5. Plasmid Transfection, Small Interfering RNA (siRNA), and Polyinosinic:polycytidylic Acid (poly (I:C)) Stimulation*

Cells were grown to 70–80% confluence and transfected with the plasmids expressing FLAG-porTRIM26 or the empty vector (p3×Flag-CMV-14 vector, as a negative control) with Lipofectamine 2000 (Thermo Fisher Scientific, Shanghai, China), according to the manufacturer's instructions. After 24 h of transfection, PIEC cells were transfected with poly (I:C) at a final concentration of 3 µg/mL (InvivoGen, Toulouse, France) and stimulated for 6 h, at which *porTRIM26* and *IFN-*β was obviously induced according to a multiple time-point sample analysis by qPCR. At 6 h after stimulation, the samples were collected for further analysis.

To investigate the role of porTRIM26, one porTRIM26 specific siRNA out of 4 (the targeted sequence: GCCTGTACCAGAGCTCTTA) was selected and transfected into PIEC cells or PAM by Lipofectamine RNAiMAX (Thermo Fisher Scientific, Shanghai, China) following the manufacture's instructions. Likewise, a scrambled siRNA (sequence: TTCTCCGAACGTGTCACGT) was transfected as the negative control (NC). After 72 h of transfection, cells were treated with poly (I:C): the PIEC cells were treated as described above, and PAM were stimulated with poly (I:C) at a final concentration of 3 µg/mL without transfection. At 6 h after stimulation, the samples were collected for further analysis.

#### *2.6. Western Blottting Analysis*

The protein samples were prepared as previously described [19]. In brief, the membranes transferred with the protein samples were blocked with 5% skim milk for 1 h at room temperature. The membrane was incubated overnight at 4 ◦C with the individual primary antibodies (anti-Flag (1:1000, M2, Sigma), anti-porTrim26 (1:1000, generated in this study), anti-VSV G (1:1000, Abcam), anti-PRRSV N (1:1000, generated by a synthetic peptide of N), and anti-actin (1:10,000, Sigma)). The membrane was then incubated for 1 h at room temperature with the secondary antibodies: horseradish-peroxidase-conjugated goat anti-mouse IgG (1:5000, Abcam) or goat anti-rabbit IgG (1:10,000, Abcam).

#### *2.7. Enzyme-Linked Immunosorbent Assay (ELISA)*

The concentrations of porcine IFN-β in supernatants from PIEC cell culture or PAM culture were measured by ELISA kit (Lengton, Shanghai, China).

#### *2.8. Reverse Transcription (RT)-Quantitative PCR (qPCR) Analysis*

Total RNA was extracted with RNAiso Reagent (Takara, Dalian, China) and the cDNA was prepared with PrimeScript RT Reagent Kit (Takara, Dalian, China). Gene expression was analyzed by RT-qPCR using SYBR Green qPCR Master Mix (Takara, Dalian, China). The specific primers are shown in Table S1. The expression of the glyceraldehyde 3-phosphate dehydrogenase gene (*GAPDH*) was used as the reference. Expression was calculated relative to that of *GAPDH (2*−∆*Ct)*.

#### *2.9. Statistical Analysis*

All data were analyzed with GraphPad Prism software (Graphpad Software, Inc, La Jolla, CA, USA). An unpaired Student's *t*-test was used to determine significant differences. Values were considered statistically significant when *p* < 0.05. Data were given as mean ± SEM as indicated; 'n' refers to the sample size.

#### **3. Results**

### *3.1. Sequence Analysis and Generation of Polyclonal Antibody against porTRIM26*

Based on the porcine *TRIM26* gene sequence (GenBank accession number NM\_001123209.1), the porcine *TRIM26* gene was amplified with RT-PCR from the total RNA extracted from pig lungs (Shanghai great white pig strain). A sequence analysis with BLAST showed that the cloned gene sequence from Shanghai great white pig strain was identical to the porcine *TRIM26* gene sequence in GenBank. The full-length cDNA of porcine *TRIM26* contained 1638 bp and encoded a protein of 546 amino acid residues including one Ring domain, two B-box domains, one coiled-coil domain, and one C-terminal domain, according to the protein BLAST information. A BLASTp analysis in the National Center for Biotechnology Information (NCBI) nonredundant database detected more than 200 TRIM26 orthologues (>50% identity) from more than 100 species. Notably, among the mammalian sequences, human and mouse TRIM26 shared 91% identity and 84% identity with porTRIM26, respectively, which is consistent with the data generated with the Clustal W method (Figure 1A). A phylogenetic analysis was performed to investigate the difference in the identities with TRIM26 of the other species. The *porTRIM26* gene clustered on a separate branch of the dendrogram within the sequences of Mammalia and was phylogenetically closest to the human sequence than to the mouse sequence (Figure 1B).

**Figure 1.** *Cont*.

**Figure 1.** Multiple alignment and phylogenetic analysis of TRIM26. (**A**) Polypeptide of porcine TRIM26 was aligned with that of human, mouse, and rat. The RING finger domain (RING), the B-box domains (B-box), and the coiled-coil region (CC) are labeled in reference to the human TRIM26. The antigenic peptide is shown in the red box. (**B**) The phylogenetic tree of TRIM26, constructed by neighbor-joining method using MEGA 6.06.

To study porTRIM26 expression in different tissues and cells of pigs, we developed an anti-porTRIM26 antibody after synthesizing an antigen peptide (Figure 1A). We tested the specificity of the antibody with a Western blotting analysis. Our result indicated that the anti-porTRIM26 antibody specifically detected the expression of porTRIM26 in HEK 293T cells (Figure 2A), and that the signal was equal to that detected with an anti-Flag antibody.

**Figure 2.** *Cont*.

**Figure 2.** (**A**) Specificity of polyclonal antibody against porcine TRIM26 was analyzed with Western blotting. (**B**) Expression profiles of TRIM26 in different porcine cell lines, analyzed with Western blotting with anti-TRIM26 antibody. (**C**) Expression profiles of TRIM26 in different pig tissues (*n* = 3), analyzed with Western blotting and an anti-TRIM26 polyclonal antibody. Approximately 100 mg of each tissue was homogenized with lysis buffer as per the Materials and Method section for protein sample preparation. (**D**) Expression profiles of *TRIM26* in different pig tissues (*n* = 3) analyzed by reverse transcription-quantitative PCR.

#### *3.2. Expression Profiles of porTRIM26 in Di*ff*erent Cell Lines and Tissues*

The polyclonal antibody was used to determine the expression profiles of porTRIM26 in different porcine cell lines and tissues in a Western blotting analysis. First, the expression of porTRIM26 in different porcine cell lines was determined by Western blotting. The PK15 cell line had a lower level of porTRIM26 in comparison to that in other cell lines (Figure 2B). The expression of porTRIM26 in different tissues was also determined with Western blotting. In almost all examined tissues porTRIM26 was differentially expressed, including in liver, spleen, lungs, kidneys, lymph nodes, and thymus (Figure 2C). The lowest expression of *porTRIM26* was observed in the kidneys in RNA level (Figure 2D), consistent with its expression in protein level (Figure 2C).

#### *3.3. porTRIM26 Negatively Regulates Expression of IFN-*β

The role of porTRIM26 in the regulation of IFN-β expression is controversial, which may be attributable to the different experimental systems, such as different cell types used in different studies. Whether porTRIM26 affects IFN-β production, either positively or negatively, remains unclear. We first investigated its effect on poly (I:C)-induced IFN-β expression using two kinds of porcine cells, PIEC cells and PAM. Notable, poly (I:C) stimulation induced the expression of porTRIM26 in both porcine cell lines (Figure 3A), consistent with the previous report in human and mouse cell lines [15]. To investigate the effect of porTRIM26 on IFN-β expression, a porTRIM26-specific siRNA was designed and used to transfect PIEC cells and PAM. A Western blotting analysis showed that the transfection of the siRNA reduced the expression of porTRIM26 (Figure 3B,C). The transfection of porTRIM26 siRNA promoted poly (I:C)-induced IFN-β expression in both PAM and PIEC cells (Figure 3B,C). In contrast, the overexpression of porTRIM26 in PIEC attenuated the poly (I:C)-induced IFN-β expression (Figure 3D). Collectively, these results reveal that porTRIM26 modulates IFN-β expression downstream of poly (I:C)-stimulated innate signaling.

β *β* ≥ **Figure 3.** Effect of porTRIM26 on poly (I:C)-induced interferon (IFN)-β. (**A**) Expression of TRIM26 in porcine alveolar macrophages (PAM) or porcine iliac artery endothelial cells (PIEC) with or without poly(I:C) treatment was measured with qPCR and Western blotting. Poly (I:C) treatment induced the expression of TRIM26 in both PAM and PIEC cells. (**B**) PAM were transfected with porTRIM26 siRNA or negative control (NC) siRNA for 72 h and then stimulated with poly (I:C) for 6 h. Supernatant and cells were harvested and analyzed by enzyme-linked immunosorbent assay (ELISA) and qPCR, respectively. (**C**) PIEC cells were transfected with porTRIM26 siRNA or NC siRNA for 72 h, and then stimulated with poly (I:C) for 6 h. Supernatant and cells were harvested and analyzed by ELISA and qPCR, respectively. (**D**) PIEC cells were transfected with plasmid encoding porTRIM26 or the empty vector. After 24 h, the cells were treated with poly (I·C) or vehicle and incubated for another 6 h. Level of *IFN-*β mRNA was determined with RT-qPCR. Data present are mean ± SEM pooled from one independent experiment; *n* ≥ 3 for each of the analyzed parameters. \*, *p* < 0.05; \*\*, *p* < 0.01; \*\*\*, *p* < 0.001 in comparison between mock and poly (I:C) group (**A**); between NC and RNAi with poly (I:C) stimulation (**B**,**C**); and between the empty vector and pFlag-porTRIM26 with poly (I:C) stimulation (**D**).

#### *3.4. porTRIM26 Negatively Regulates Antiviral Response to VSV Infection*

Although we had shown that porTRIM26 plays a role in poly (I:C)-stimulated innate immune signaling, whether it plays a role in virus-triggered signaling remained to be determined. In previous studies, VSV has been used to study the effect of TRIM26 on IFN-β expression and antiviral response. Because PAM are resistant to VSV infection according to our preliminary experiments, we next investigated the effect of porTRIM26 (either positive or negative) on VSV infection in PIEC cells. Similar

to poly (I:C) stimulation of PIEC cells, VSV infection also induced the expression of porTRIM26 in PIEC cells (Figure 4A). To identify the role of porTRIM26 in VSV infection we transfected PIEC cells with porTRIM26 siRNA, as described above. The knock-down of porTRIM26 expression significantly increased the expression of IFN-β during VSV infection compared with that in cells infected with the negative control. The viral titers were higher in the VSV infected-negative control cells than in VSV infected cells in which the expression of TRIM26 was knocked down (Figure 4B). To confirm these results, we transfected PIEC cells with the plasmid, as described above, to overexpress porTRIM26. ELISA data showed that the overexpression of porTRIM26 significantly reduced the expression of IFN-β relative to that in the VSV-infected empty-vector-transfected group. Notable, VSV infection was significantly greater in the porTRIM26-overexpressed cells than in the VSV-infected empty-vector-transfected cells. Taken together, these data suggest that porTRIM26 promotes viral infection by inhibiting IFN-β expression.

β ≥ **Figure 4.** Effect of porTRIM26 on IFN-β expression and vesicular stomatitis virus (VSV) infection in PIEC cells. (**A**) PIEC cells were infected with VSV at a multiplicity of infection (MOI) of 1 for 24 h. The porTRIM26 mRNA and protein level was measured by RT-qPCR and Western blotting, respectively. (**B**) PIEC cells were transfected with porTRIM26 siRNA or NC siRNA. After 72 h, the cells were mock infected or infected with VSV at a MOI of 1. Supernatants were collected at 24 h post infection (hpi) for tissue culture infective dose (TCID50) assay. (**C**) PIEC cells were transfected with plasmid encoding porTRIM26 or empty vector for 24 h and then infected with VSV at an MOI of 1. Supernatants were collected at 24 hpi for either a TCID<sup>50</sup> assay or ELISA. Data are mean ± SEM pooled from one independent experiment; *n* ≥ 3 for each of the analyzed parameters. ND, no detected. \*\*, *p* < 0.01; \*\*\*, *p* < 0.001 in comparison: mock-infected vs. VSV-infected group (**A**); NC-treated vs. siRNA-treated groups after VSV infection (**B**); empty-vector-transfected vs. pFlag-pTRIM26-transfected after VSV infection (**C**).

#### *3.5. porTRIM26 Negatively Regulates Antiviral Response to PRRSV Infection*

Our data showed that porTRIM26 inhibits IFN-β production and innate antiviral response. Therefore, we examined whether this also occurs during other viral infections. A previous study showed with an RNA-sequencing (RNA-seq) analysis that PRRSV infection could induce the expression of porTRIM26 [20], so we examined the biological functions of porTRIM26 using PRRSV, to which PAM are susceptible. Our data show that PRRSV infection induced the expression of TRIM26 in PAM (Figure 5A), consistent with the previous RNA-seq data. We also investigated the role of porTRIM26 in IFN-β production and the antiviral response in PRRSV-infected PAM. Our ELISA data showed that knocking down the expression of porTRIM26 increased the expression of IFN-β compared with that in the PRRSV-infected negative control (NC) cells (Figure 5B). We also noted that the viral titer was significantly lower in the PRRSV-infected porTRIM26-siRNA-transfected cells than in the PRRSV-infected negative control (NC) cells. These data suggest that PRRSV may inhibit IFN-β production and the antiviral response by inducing porTRIM26. These results confirm that porTRIM26 plays an inhibitory role in IFN-β expression and the antiviral response.

β ≥ **Figure 5.** Effect of porTRIM26 on IFN-β expression and porcine reproductive and respiratory syndrome virus (PRRSV) infection in PAM. (**A**) PAM were infected with PRRSV at a MOI of 1 for 24 h. The mRNA levels of *TRIM26* and PRRSV *ORF7* were measured with RT-qPCR. The protein levels of TRIM26 and PRRSV N were measured with Western blotting. PRRSV infection induced expression of TRIM26. (**B**) PAM were transfected with the porTRIM26 siRNA or NC siRNA. After 72 h, cells were mock-infected or infected with PRRSV at an MOI of 1. Supernatants were collected at 24 hpi for either a TCID<sup>50</sup> assay or ELISA. The knock-down of porTRIM26 expression increased the production of IFN-b and inhibited PRRSV infection. Data are mean ± SEM pooled from one independent experiment; *n* ≥ 3 for each of the analyzed parameters. ND: not detected. \*\*, *p* < 0.01; \*\*\*, *p* < 0.001 in comparison: mock-infected vs. PRRSV-infected group (**A**); NC-siRNA-transfected vs. *porTRIM26*-siRNA-transfected cells after infection with PRRSV (**B**).

#### **4. Discussion**

Several studies have reported the gene sequence of *porTRIM26* based on genomic sequencing analyses [12,13]. However, the expression profiles of porTRIM26 in tissues and its biological functions have been unclear until now. In this study, we confirmed that porTRIM26 inhibits IFN-βproduction after poly (I:C) stimulation or viral infection. Using a rabbit anti-pTRIM26 polyclonal antibody generated with synthesized antigen peptide, we showed that poly (I:C) stimulation or viral infection (VSV or PRRSV) induces the expression of porTRIM26. By overexpressing or knocking down the expression of

porTRIM26, we demonstrated that the induction of porTRIM26 modulates the IFN-βexpression induced by poly (I:C) stimulation. Our data demonstrate that the induction of porTRIM26 negatively regulates IFN-β production and the antiviral response to VSV or PRRSV infection. Collectively, these data provide novel evidence that porTRIM26 modulates the innate antiviral response by inhibiting IFN-β production.

Although more than 50 *TRIM* genes of the pig have been annotated, only TRIM21 has yet been functionally analyzed. Swine TRIM21 restricts foot-and-mouth disease virus (FMDV) infection with an intracellular neutralization mechanism [21]. Our data demonstrate that porTRIM26 negatively regulates IFN-β production and the antiviral response to PRRSV infection. These findings not only verify the role of porTRIM26 in IFN-β expression and the antiviral response but also extend our understanding of how PRRSV uses host proteins, such as porTRIM26, to interfere with the innate antiviral response.

A previous study has shown that TRIM26 negatively regulates IFN-β production and the innate antiviral response by inhibiting activation of interferon regulatory factor 3 (IRF3) [15]. Because PRRSV infection interferes with the activation of IRF3, it cannot be excluded that PRRSV infection inhibits the activation of IRF3 by inducing porTRIM26. The exact mechanism by which PRRSV inhibits IFN-β production via induction of porTRIM26 needs to be clarified in future research.

Although TRIM26 negatively regulates IFN-β production and the innate antiviral response, one study has reported that TRIM26 has the opposite function, promoting IFN-β production and the innate antiviral response [14]. It is possible that the different experimental systems used, including different methods and viruses, have generated different results. For example, TRIM21 is reported to degrade IRF3, thus limiting the type I IFN response after Japanese encephalitis virus (JEV) infection [22]. In contrast, another report suggested that TRIM21 acts as a positive regulator of the IRF3 pathway during Sendai virus infection [23]. Whether or not porTRIM26 positively affects IFN-β production and the innate antiviral response needs to be investigated in other swine viruses. However, our results suggest that porTRIM26 modulates IFN-β production and the innate antiviral response.

In summary, this is the first study to identify the biological functions of porTRIM26. We have demonstrated that porTRIM26 plays an inhibitory role in IFN-β expression and the innate antiviral response. These findings extend our understanding of how some swine viruses, such as PRRSV, inhibit IFN-β production to evade the host's innate immune response.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/10/1226/s1, Table S1: PCR primer used in the study.

**Author Contributions:** H.H. and M.S.: Conceptualization, methodology, writing–original daft, visualization. Y.Z. and C.L.: Conceptualization, visualization, methodology. K.L., J.W., D.S. and B.L.: Methodology, visualization. R.C. and Z.M.: Supervision. Y.Q.: Conceptualization, supervision, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported, in part, by the national key Research and Development program of China (2018YFD0500100), the National Natural Science Foundation of China (31972693), the Chinese Special Fund for Ago-scientific Research in the Public Interest (No. 2020JB04), and the Elite program of CAAS (to YQ).

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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*

### **3** ′**quant mRNA-Seq of Porcine Liver Reveals Alterations in UPR, Acute Phase Response, and Cholesterol and Bile Acid Metabolism in Response to Di**ff**erent Dietary Fats**

#### **Maria Oczkowicz 1, \* , Tomasz Szmatoła 1,2 , Małgorzata Swi ˛atkiewicz ´ 3 , Anna Koseniuk 1 , Grzegorz Smołucha 1 , Wojciech Witarski <sup>1</sup> and Alicja Wierzbicka 1**


Received: 20 August 2020; Accepted: 16 September 2020; Published: 18 September 2020

**Abstract:** Animal fats are considered to be unhealthy, in contrast to vegetable fats, which are rich in unsaturated fatty acids. However, the use of some fats, such as coconut oil, is still controversial. In our experiment, we divided experimental animals (domestic pigs) into three groups differing only in the type of fat used in the diet: group R: rapeseed oil (*n* = 5); group B: beef tallow (*n* = 5); group C: coconut oil (*n* = 6). After transcriptomic analysis of liver samples, we identified 188, 93, and 53 DEGs (differentially expressed genes) in R vs. B, R vs. C, and B vs. C comparisons, respectively. Next, we performed a functional analysis of identified DEGs with String and IPA software. We observed the enrichment of genes engaged in the unfolded protein response (UPR) and the acute phase response among genes upregulated in B compared to R. In contrast, cholesterol biosynthesis and cholesterol efflux enrichments were observed among genes downregulated in B when compared to R. Moreover, activation of the UPR and inhibition of the sirtuin signaling pathway were noted in C when compared to R. The most striking difference in liver transcriptomic response between C and B was the activation of the acute phase response and inhibition of bile acid synthesis in the latest group. Our results suggest that excessive consumption of animal fats leads to the activation of a cascade of mutually propelling processes harmful to the liver: inflammation, UPR, and imbalances in the biosynthesis of cholesterol and bile acids via altered organelle membrane composition. Nevertheless, these studies should be extended with analysis at the level of proteins and their function.

**Keywords:** pigs; fatty acids; 3 ′quant mRNA-seq; nutrigenomics

#### **1. Introduction**

The study of lipid metabolism is becoming increasingly important in the context of the growing incidence of human metabolic diseases like obesity, NAFLD (Non-Alcoholic Fatty Liver Disease), and type 2 diabetes, as well as neurodegenerative disorders and cancer [1,2]. Dietary fats are one of the most critical determinants of the vulnerability of organisms to the development of diseases. Fatty acids are essential components of cell membranes and function as signaling molecules, regulating enzyme activity and preserving homoeostasis [3].

The recent development of advanced sequencing techniques accelerated research on the interactions between nutrients and diet. The RNA-seq results have given insight into the lipid metabolism of many species, including mice, humans, and pigs. The domestic pig is widely used as a model animal in biomedical research. It seems to be perfect for nutrigenomic research since, as an omnivore, it allows researchers to test various diets. Additionally, genetically, it is much more diverse than laboratory rats or mice, and thus better reflects human biological processes. Thus, this species is commonly used as a model in studies of the molecular background of some human diseases, including cancer and atherogenesis [4,5].

Recently, several studies investigated the effect of omega-3 and omega-6 in the diet of pigs on the transcriptome in different tissues [6–8], showing that different fatty acid compositions have a significant effect on the expression of genes engaged in fatty acid synthesis and metabolism. We also performed a preliminary study on the impact of different sources of fat in the diet on the liver transcriptome in pigs [9]; however, due to the limited number of samples and identified differentially expressed genes, we were not able to recognize all aspects of the observed effects. In the present paper, we described the comprehensive functional analysis of Quant-seq 3′mRNA-seq quantitative profiling (Lexogen, Vienna, Austria) of liver samples collected from animals receiving rapeseed oil, beef tallow, or coconut oil in the diet.

Although dietary recommendations point to the harmful effects of excessive consumption of saturated fat, there is still a lot of doubt about the health effects of consuming fats from various sources, e.g., coconut or rapeseed oil. Our investigation aimed to identify the possible molecular mechanism of metabolic changes that occur after receiving different fats in the diet and how they contribute to the pathogenesis of diseases.

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

#### *2.1. Animals and Diets*

Animals for the study were kept at the Testing Station of the National Research Institute of Animal Production in Grodziec Sl ˛aski. The experiment described in this manuscript was conducted on ´ animals used in previously published studies [10–13]. All procedures included in this study relating to the use of live animals were in agreement with the local Ethics Committee for Experiments with Animals in Cracow (Resolution No. 912 dated 26 April 2012). In this study, we used 16 samples of liver collected from crossbred fatteners divided into three dietary groups, in which the diets differed among each other in terms of fodder fat: 3% rapeseed oil (group R−*n* = 5, 3 gilts and 2 barrows), 3% beef tallow (group B−*n* = 5, 2 gilts and 3 barrows), and 3% coconut oil (group C−*n* = 6, 3 gilts and 3 barrows). Two samples were excluded from the initial number of 18 samples, due to the low quality of RNA-sequencing. All animals were kept in individual straw-bedded pens in uniform conditions. The animals were healthy, and as equal as possible in regards to body weight. The diets of all of the groups were isonitrogenous and isoenergetic (metabolized energy (ME): *R* = 13.4, *B* = 13.3, *C* = 13.4 MJ/kg feed), and were formulated to cover the nutritional requirements of the pigs. The ingredient composition and nutritive value of the diets are presented elsewhere [11], while the fatty acid compositions of the feed mixtures are shown in Figure 1. Briefly, the group I feed mix contained 80% UFA content (44% MUFA and 36% PUFA), group B contained 67% UFA (32% MUFA and 35% PUFA), and group C contained 45% UFA (16% MUFA and 29% PUFA). The experimental fattening lasted from 60 to 118 kg of live weight of the animals. At the end of the experiment, all the pigs were slaughtered by stunning with high-voltage electric tongs (voltage 240–400 V), and samples of subcutaneous adipose tissue from the area between the last thoracic and the first lumbar vertebrae were collected for transcriptome analysis. All samples were stored in a freezer (−85 ◦C) until analysis.

**Figure 1.** Percentage of the most abundant fatty acids in the diets.

#### *′ 2.2. RNA Isolation and 3*′*Quant mRNA Library Construction and Sequencing*

′ RNA isolation was performed using the SPLIT RNA Extraction Kit (Lexogen, Vienna, Austria) according to the manufacturer's recommendations. Quality of RNA was assessed using a Tapestation 2200 (Agilent, Santa Clara, CA, USA), while quantity was evaluated by nanodrop. A sample of 200 ng of RNA was used for library preparation using the QuantSeq 3′ mRNA-Seq Library Prep Kit FWD for Illumina (Lexogen, Vienna, Austria) according to the manufacturer's protocol. Assessment of library quantity and quality was performed using a Qubit (Thermofisher Scientific, Foster City, CA, USA) and Tapestation 2200 (Agilent, Santa Clara, CA, USA), respectively. The sequencing of the pooled libraries (50 bp single-read) was performed on an Illumina Hiseq 2500 instrument (Illumina, San Diego, CA, USA), using the High-Output v4—SR 50 Cycle kit (Illumina, San Diego, CA, USA) at the DNA Sequencing Center at Brigham Young University (Provo, UT, USA).

#### *2.3. Bioinformatic Analysis*

Demultiplexed fastq files were downloaded from the sequencing provider server. Next, the quality check, trimming of reads, and mapping of reads were conducted with FastQC 11.8, FLEXBAR 3.5.0, and TopHat 2.1.1 software, respectively. Samtools 1.9, RSeQC, and HTSeq-count 0.11.1 software, and Gtf-Ensembl annotation 96 were used for evaluation of the mapping statistics and read counts. Differential expression analysis was performed using DEseq 2 software. Genes with *p*-adjusted < 0.05 (FDR-False Discovery Rate) Benjamini–Hochberg (BH) adjustment and no fold-change threshold were regarded as differentially expressed. Functional analysis with STRING software was performed separately for upregulated and downregulated genes using the *Sus scrofa* database. Functional analysis with IPA software was performed using porcine gene names with databases for all available species (human, rat, mouse). Figures were made with the Biorender.com and Paint programs.

#### *2.4. Quantitative PCR*

RNA was reverse transcribed using a high capacity cDNA archive kit (Thermo Fisher Scientific). Next, qPCR was performed using TaqMan gene expression assays and TaqMan gene expression master mix on a QuantStudio 7-flex instrument. Relative quantity data were analyzed on the Thermo Fisher cloud. Statistical significance was assessed using the Relative Quantification application on Thermo Fisher Connect.

#### **3. Results**

#### *3.1. Fatty Acids Profiles of Diets Used in the Experiment*

In our experiment, we compared the expression of genes in the liver of pigs fed with three isoenergetic diets differing only in the type of fat. The detailed information about the diets are presented elsewhere [10]. The proportions of the most important fatty acids in each diet are presented in Figure 1. Group R (rapeseed oil) obtained feed with the highest amounts of unsaturated fatty acids (UFA) and the lowest quantity of saturated fatty acids (SFA). The predominant fatty acids in this group were behenic acid (C22:0), eicosapentaenoic acid (EPA) (C20:5), docosahexaenoic acid (DHA) (C22:6), and α-linolenic acid (C18:3). The second group had an intermediate ratio of SFA/UFA and high content of palmitic acid (C16:0), stearic acid (C18:0), arachidonic acid (C20:4), and palmitoleic acid (C16:1). Group C displayed the highest SFA/UFA ratio and a high content of short/medium-chain saturated fatty acids: myristic acid (C14:0), capric acid (C10:0), and lauric acid (C12:0).

Diets used in the study did not change the phenotype of animals dramatically [10]. We observed no differences in the weight of animals, the weight of the liver and main muscles, or backfat thickness. Different diets did not affect average feed utilization or growth of animals. However, substantial differences were observed in the fatty acid composition of the adipose tissue. We found a high correlation between the dietary fatty acid composition of the diet and adipose tissue fatty acid composition [10–12].

### *3.2. 3*′*Quant mRNA Statistics and DEGs Identified in the Liver after Di*ff*erent Dietary Treatments*

After 3′quant mRNA-sequencing, we obtained, on average, 2,776,234 raw reads per sample. More than 79% of them were mapped to *Sus scrofa* 11.1. (Table 1). On average, 11,500 genes with read counts >1 were identified in each sample. All gene expression data were submitted to the GEO NCBI database (accession number: GSE144247). Using DESeq2 software, we performed comparisons of gene expression between all three dietary groups (rapeseed oil vs. beef tallow, rapeseed vs. coconut oil, beef tallow vs. coconut oil, the first group in the comparison is the reference group). MA plots and PCA obtained after DESeq2 analysis are presented in Supplementary Figure S1. We identified 188 DEGs in the group R vs. group B comparison (77 downregulated in group B and 111 upregulated in group B (adjusted *p*-value < 0.05). In the comparison of group R with group C, 93 DEGs were identified (45 downregulated in group C and 48 upregulated in group C). The lowest number of DEGs (53) was noted in the group B versus group C comparison (32 downregulated in group C and 21 upregulated in group C) (Supplementary Table S1. Twenty genes with the highest adjusted *p*-value from each comparison are listed in Table 2.


**Table 1.** RNA-seq statistics of Quant-seq 3′mRNA-sequencing.

**Table 2.** Top 20 DEGs with the highest adjusted *p*-value for each comparison.





**Table 2.** *Cont.*

Bold genes were annotated with the Uniprot database.

Among the identified DEGs, there were 30 genes in common between comparison groups R vs. B and groups R vs. group C, 25 common genes between comparison groups R vs. B and groups B vs. C and eight genes in common between the group R vs. C comparison and the B vs. C comparison (Figure 2). Two differentially expressed genes (EGFR, HSPA5) were common for all three comparisons. Interestingly, expression of all common genes was altered in the same direction in both groups: in B and C when compared to R, and in R and C when compared to B. All identified genes in common between specific comparisons are listed in Supplementary Table S2.

**Figure 2.** Venn diagram illustrating differentially expressed genes that are common between comparisons.

Since there was a sex imbalance in some groups in our study, we performed a comparison between females and males with DEseq2 software, regardless of dietary group. We observed only 11 differentially expressed genes between female and males, and all of them were located on sex chromosomes (Supplementary Table S3). What is more, none of the genes differentially expressed between males and females were present among the genes identified as differentially expressed between dietary groups. We concluded that the effect of sex is negligible in our study.

#### *3.3. Functional Analysis of Identified DEGs with STRING*

To get insights into the biological processes that are activated by specific dietary fats, we performed functional analysis of identified DEGs with String software using *Sus scrofa* genes as a background. We separately analyzed upregulated and downregulated genes in each comparison. We observed 86, 82, and 18 enrichments in rapeseed oil vs. beef tallow, rapeseed oil vs. coconut oil, and the beef tallow vs. coconut oil comparisons, respectively (Supplementary Table S4). In R vs. B and R vs. C comparisons, we observed a higher number of enrichments among downregulated genes (Supplementary Table S4). Among genes downregulated by beef tallow in the R vs. B comparison, we observed overrepresentation of genes engaged in oxidoreductase activity (*CAT*, *MSMO*, *GLUD1*), cholesterol biosynthesis (e.g., *HSD3B1*, *HMGCS1*, *MVK*), metabolism of steroids (*FDFT1*, *MVK*, *HMGCS1*), and bile acid secretion (*CYP7A1*, *HMGCR*) (Figure 3, Supplementary Table S4). Furthermore, among genes upregulated by beef tallow in this comparison, enrichment of genes associated with the metabolism of RNA (*GNL3*, *NOP58*, *NOPB*), cellular response to stress (*HSPA5, HSPA8*, *BAG3*), the urea cycle (*ENSSSCG00000016159*, *SLC25A15*), and metabolism of polyamines (*AMD1*, *ENSSSCG00000016159*, *SLC25A15*) were observed (Figure 4, Supplementary Table S4).

**Figure 3.** The network of DEGs downregulated in the beef tallow group when compared to the rapeseed oil group according to STRING software analysis.

**Figure 4.** The network of DEGs upregulated in the beef tallow group when compared to the rapeseed oil group according to STRING software analysis.

In the second comparison (R vs. C), we observed that genes upregulated by coconut oil are overrepresented in pathways connected to stress response and cellular signaling: the regulation of HSF1-mediated heat shock response (*BAG3*, *ENSSSCG00000015140*, *HSPA9*), and NOTCH3 activation and transmission of signal to the nucleus (*EGFR, PSEN2*). Moreover, downregulated genes were enriched in many metabolism-related pathways and terms (e.g., primary metabolic process, organic substance metabolic process, cellular metabolic process, *ANPEP*, *C1QA*, *ERH*, *FTCD*, *HSD3B1*, *MAN1A1*, *RPS3*, *SDHD*, *TMEM86B*, *UOX*).

In the comparison between beef tallow and coconut oil, we observed downregulation of genes engaged in immunity (defense response (*C1QA, HP, IL4R, ITIH4*), and genes engaged in protein processing in the endoplasmic reticulum (*HSPA5, HYOU1, PDIA3*). Genes upregulated by coconut oil were associated with PPAR signaling pathways (*CYP7A1, FABP1)* and bile secretion and cholesterol metabolism (*ABCB11, CYP7A1)*. All enrichments with adequate FDR values and matching proteins are presented in Supplementary Table S4.

#### *3.4. Functional Analysis of Identified DEGs with IPA*

To further analyze differentially expressed gene sets for their association with human diseases and to assess to what extent the relationships observed in pigs can be translated into humans, we performed functional analysis of identified DEGs from all comparisons with ingenuity pathway analysis (IPA) (Qiagen). Sixty-six significant canonical pathways were observed in the comparison between group R vs. B (log2 *p*-value > 1.5) (Supplementary Figure S2). Among them, the superpathway of cholesterol biosynthesis (Z-score = −2.646, *p*-value < 2.34 × 10<sup>9</sup> , mevalonate pathway I (Z-score = −2, *p*-value < 3.57 × 10<sup>6</sup> ), the superpathway of geranylgeranyldiphosphate biosynthesis I (Z-score = −2, *p*-value < 1.07 × 10<sup>5</sup> ), and LXR/RXR activation (Z-score = −0.447, *p*-value < 5.12 × 10<sup>5</sup> ) were inhibited, while LPS/IL-1-mediated inhibition of RXR function (Z-score = 1.342, *p*-value < 1.25 × 10<sup>5</sup> ), acute phase response signaling (Z-score = 2.236, *p*-value < 5.74 × 10<sup>4</sup> ), Huntington's disease signaling (Z-score = 1, *p*-value < 2.88 × 10<sup>3</sup> ), and several pathways of inositol metabolism were activated by the beef tallow diet when compared to the rapeseed oil diet (Figure 5) (Supplementary Figure S2). In the second comparison (rapeseed oil vs. coconut oil), 30 canonical pathways were noted, with the sirtulin signaling pathway (Z-score = −2, *p*-value < 1.1 × 10<sup>4</sup> ) being inhibited and the BAG2 signaling pathway (Z-score = 1, *p*-value < 2.01 × 10<sup>5</sup> ) and unfolded protein response (Z-score = 2, *p*-value < 5.76 × 10<sup>5</sup> ) being activated by coconut oil (Figure 5). The beef tallow vs. coconut oil comparison revealed 43 enriched canonical pathways. Beef tallow activated acute phase response signaling (Z-score = −2.236, *p*-value < 4.32 × 10<sup>9</sup> ), while coconut oil stimulated RXR/LXR activation (Z-score = 0.816, *p*-value < 7.24 × 10<sup>9</sup> ) in this comparison (Figure 5). All identified canonical pathways and connections between them are presented as networks in Supplementary Figure S2.

#### *3.5. Identification of Hub Genes with Cytohubba*

Our next step was to identify critical genes responsible for changes observed in the transcriptome under the influence of different types of fat. For this purpose, we used the cytoHubba–Cytoscape plugin for ranking nodes in a network by their network features. Among eleven available methods, we chose MCC, which is the most precise in predicting essential proteins from the yeast PPI network [13]. In the first comparison, R vs. B, *GNL3* was ranked number one, followed by *RSL1D1, UTP3*, and *DDX24*. The essential genes in the B vs. C comparison were *ORM1* and *TTR*, while in the rapeseed oil vs. coconut oil they were *RPS3* and *HSPA9* (Figure 6).

### *3.6. Validation of Quant 3*′*mRNA Profiling by qPCR*

To validate the 3′quant mRNA-seq results, we performed qPCR analysis of several DEGs identified in the study. We observed high concordance between RNA-seq and qPCR results. As expected, we found strong overexpression of genes engaged in the unfolded protein response (*BAG3*, *HSPA5*) and inflammation (*PID1*, *IHIT4*, *ALPL*, *LITAF*) in the beef tallow group compared to the rapeseed oil and coconut oil groups. In contrast, genes involved in bile acid secretion (*CYP7A1*) and cholesterol biosynthesis *(HSD3B1)* were downregulated in the beef tallow group when compared to both groups or the rapeseed oil group only (Figure 7). Furthermore, none of the genes analyzed by qPCR showed significantly different expression between gilts and barrows.

**Figure 6.** Top 10 Hub genes identified in the R vs. B comparison, R vs. C comparison, and B vs. C comparison using Cytohubba, ranked by MCC (Maximal Clique Centrality). The more intense the red color, the higher the position in the rank.

**Figure 7.** Results of the qPCR analysis of selected DEGs. RQ—Relative Quantity of mRNA, \*\*\* *p* < 0.01, \*\* *p* < 0.05, \* *p* < 0.1, ns—not significant.

#### **4. Discussion**

The progress of civilization has initiated many changes in the human environment. This applies not only to nature that surrounds us but also to our way of life, especially our nutrition. First of all, fat consumption and the amount of saturated fatty acids consumed increased significantly during the last

′

′

century. Moreover, the ratio of omega-6 to omega-3 in food decreased, resulting in worsening of the dietary fatty acid profile [14]. At the same time, there was a sharp increase in the incidence of diseases related to metabolism, cardiovascular diseases, neurodegenerative diseases, and cancer. These diseases are also called civilization diseases because their occurrence is closely related to the modern lifestyle.

In our experiment, we used a domestic pig (*Sus scrofa*) to illustrate the changes that occur in the liver transcriptome as a result of using various fats in the diet. Using isoenergetic and isonitrogenous diets made it possible to observe transcriptome changes at a specific stage, before drastic changes at the phenotype level occurred. For quantitative analysis of the transcriptome, we chose the 3′quant mRNA method, which allows transcriptome profiling based on the 3′UTR ends of the gene and is a cost-efficient alternative to whole transcriptome RNA-seq. By this method, we identified 11,500 genes with >1 read number, of which 308 (2.67%) were differentially expressed between samples of animals fed different diets, despite a very low sequencing coverage—approximately two million reads per sample. We observed that many of these genes are engaged in the pathogenesis of human civilization diseases and that depending on the source of dietary fat, different pathways are activated or inhibited in liver at the gene expression level.

#### *4.1. Biosynthesis and Catabolism of Cholesterol are Inhibited in Animals Obtaining Beef Tallow in the Diet*

The opinion that the consumption of saturated fats increases blood cholesterol levels is generally approved and well documented [15]. Maintaining adequate blood cholesterol levels is crucial for the body since its excess results in the development of atherosclerosis, heart disease, and neurodegenerative diseases. On the other hand, cholesterol is an essential component of cell membranes and has many important functions in the organism. Cholesterol homeostasis in the body is directed by the interaction between absorption, synthesis, and excretion or conversion of cholesterol into bile acids. A reciprocal relationship between these processes is known to regulate circulating cholesterol levels in response to dietary or therapeutic interventions. Cholesterol biosynthesis is self-regulating; a high cholesterol level in the blood forces the inhibition of its synthesis in the liver. One of the mechanisms involved in this process is inhibition of the expression of the gene encoding HMGCR, a rate-limiting enzyme in cholesterol biosynthesis. This is exactly what we observed in animals fed beef tallow—the expression of *HMGCR* was ~3 fold lower than in animals fed rapeseed oil. We observed inhibition of several other genes engaged in cholesterol biosynthesis (*FDFT1*, *MVD*, *EBP*, *HSDB31*) in the group receiving beef tallow in the diet when compared to the group receiving rapeseed oil. Two of these genes (*FDFT1*, *HSD3B1*) were also downregulated in the coconut oil group compared to the rapeseed oil group but to a lesser extent, suggesting a different response of the mechanisms responsible for cholesterol homeostasis to medium-chain saturated fatty acids. The results of the functional analysis using the STRING and IPA programs indicate the existence of an additional mechanism affecting the level of cholesterol synthesis. According to the IPA results (Figure 5), RXR function was inhibited through LPS/IL-1 mediation, and simultaneously the super pathway of cholesterol biosynthesis was repressed in the beef tallow group. We suppose that it may have occurred as a result of a series of changes triggered by gut microbiome misbalance under the influence of dietary beef tallow. A diet rich in saturated fatty acid affects gut microbiota composition by enhancing overflow of dietary fat to the distal intestine in mice [15]. In the pigs fed beef tallow, an increase in pathogenic LPS-secreting bacteria appeared, which resulted in an increase in expression of *LBP* and other acute-phase response genes in the liver. Moreover, beef tallow contains substantial amounts of arachidonic acid, which is known for its pro-inflammatory properties. As a consequence of inflammation, the RXR function was inhibited, which affected the level of expression of genes engaged in cholesterol metabolism (Figure 8).

**Figure 8.** Graphical illustration of selected biological processes and canonical pathways altered in the beef tallow group when compared to rapeseed oil. Beef tallow contains pro-inflammatory ingredients (SFA, AA—arachidonic acid) that could change gut microbiota and promote inflammation. SFA and cholesterol from beef tallow decrease cholesterol biosynthesis directly in the liver and indirectly through inhibition of LXR/RXR by LPS/IL-1.

The use of coconut oil in the diet still causes a lot of controversies. As a fat containing many saturated fatty acids, and thus causing an increase in blood cholesterol levels, it is not recommended for people at risk of cardiovascular disease. On the other hand, it was shown recently that it may improve intestinal microbiota, antioxidant status, and immunity of growing rabbits [16]. Moreover, the latest meta-analysis showed that MCFA (medium-chain fatty acids), which predominate in coconut oil, increase HDLcholesterol—responsible for cholesterol efflux—content in comparison with long-chain fatty acids. *CYP7A1* is a gene coding for a rate limiting enzyme in the cholesterol catabolic pathway in the liver, which converts cholesterol to bile acids. This reaction is the major site of regulation of bile acid synthesis, which is the primary mechanism for the removal of cholesterol from the body. We observed a 6-fold decrease in *CYP7A1* gene expression in the beef tallow group compared to coconut oil, emphasizing the advantage of coconut oil over beef tallow in cholesterol efflux. It was shown in vitro that arachidonic acid is a potent inhibitor of *CYP7A1* expression [17], which is in accordance with our results.

#### *4.2. Acute-Phase Response Signaling is Activated in the Beef Tallow Group when Compared to both Rapeseed Oil and Coconut Oil*

Low-grade inflammation is one of the leading causes of NAFLD, cardiovascular diseases, diabetes, and neurodegenerative diseases. We observed a significant increase in the expression level of many genes connected to the immune system response and inflammation markers in liver tissue collected from animals fed with beef tallow (*LITAF*, *ALPL*, *PID1*, *IHIT4*, *LBP*, *HP*) (Figure 7, Supplementary

Table S4). Strikingly, the *LBP* gene, which codes for lipopolysaccharide binding protein, and *CD14* were both upregulated in the beef tallow group. LBP and CD14 drive ternary complex formation and TLR activation and as a consequence trigger the whole cascade of immune response stimulated by *NFKBiB* [18]. Interestingly, the expression of some genes (*LBP*, *IHIT4*) was also upregulated in rapeseed oil when compared to coconut oil, supporting information about the antibacterial properties of coconut oil [16]. Ingenuity pathway analysis revealed that the acute phase response canonical pathway is highly significantly activated in animals obtaining beef tallow when compared to both rapeseed oil and coconut oil. The pro-inflammatory effects of long-chain saturated fatty acids have been known for some time [19]. It has been observed that long-chain saturated fatty acids increase haptoglobin gene expression (inflammation marker) in mice adipose tissue [19]. Another study showed that the composition and metabolic activity of the gut microbiota change as a result of a steatohepatitis-inducing high-fat diet in mice [20]. Moreover, in these animals, the level of saturated fatty acids (palmitic acid) in the gut increased significantly, activated macrophages in the liver, and promoted TNF-α expression. The consequence of these changes was the development of NASH, which was reversible under the influence of antibiotics [20]. On the other hand, lauric acid—the main component of coconut oil—was shown previously to alleviate neuroinflammatory responses by LPS-activated microglia, supporting its beneficial effect on neurodegenerative diseases [21]. Our results, in agreement with previous studies, underline the difference in response of the immune system to dietary long-chain SFA (beef tallow) and medium-chain SFA (coconut oil), (Figure 9).

**Figure 9.** Graphical illustration of selected biological processes and canonical pathways altered in the coconut oil group when compared to the rapeseed oil (UPR, sirtuin signaling pathway) and beef tallow groups (LXR/RXR activation, cholesterol metabolism, bile secretion). Coconut oil contains antibacterial and anti-inflammatory ingredients (lauric acids, polyphenols), and cholesterol metabolism and bile acid secretion are not reduced as a result of increased inflammation, which is observed in the beef tallow group.

#### *4.3. Unfolded Protein Response is Activated in both the Beef Tallow and Coconut Oil Groups when Compared to the Rapeseed Oil Group*

The unfolded protein response (UPR) is a highly conserved pathway that allows the cell to manage endoplasmic reticulum (ER) stress that is imposed by the secretory demands associated with environmental forces [22]. Our results show that upon dietary beef tallow or coconut oil uptake, expression of genes engaged in UPR—*BAG3*, *HSP5*, *HSP7*—increases when compared to rapeseed oil. Consequently, functional analysis with IPA and STRING revealed that UPR is strongly activated in the coconut oil and beef tallow groups when compared to rapeseed oil. qPCR analysis confirmed overexpression of *HSPA5* and *BAG3* in beef tallow and both saturated fatty acids rich groups, respectively, when compared to rapeseed oil. All these results support the opinion that UPR is regulated by lipid-dependent mechanisms [23,24]. It is assumed that in an environment rich in saturated fatty acids, the composition of membranes in the endoplasmic reticulum changes, leads to disturbances in protein folding and results in activation of the unfolded protein response pathway [24]. Recent in vitro studies [25] indicate that palmitic acid induces ER stress and simultaneously increases inflammatory indices, while oleic acid ameliorates this action in exocrine pancreas cells. In our experiment, the ratio of palmitic to oleic acid was about twice as high in coconut oil and beef tallow as compared to rapeseed oil (Figure 1), supporting in vivo associations previously observed in vitro in pancreas cells [25].

When the proteins present in mitochondria are damaged, and their accumulation threatens the maintenance of homeostasis, mitochondrial specific UPR (UPRmt) is activated. The central UPRmt coordinator is *SIRT3*, which encodes a member of the sirtuin family of class III histone deacetylases. In our study, we observed lower expression of *SIRT3* in animals fed coconut oil in the diet when compared to the rapeseed oil group. Moreover, according to IPA analysis, the sirtuin signaling pathway (with *SIRT3*, *SDHD*, *UOX* genes) was downregulated in the coconut oil group. This result may have an important clinical significance since the loss of *SIRT3* leads to deregulation of several mitochondrial pathways, which contributes to the accelerated development of the disease of ageing [26]. In general, lowering the *SIRT3* level is associated with adverse health effects. It is considered a mitochondria-localized tumor suppressor, which opposes reprogramming of cancer cell metabolism through HIF1α destabilization [27]. Moreover, it was shown that *SIRT3* deficiency accelerates the development of metabolic syndrome. On the other hand, in vivo experiments show that a chronic high-fat diet decreases expression of *SIRT3* in liver tissue [28]. In contrast, in vitro experiments suggest that palmitic acids increase *SIRT3* expression, contrary to oleic acids [29]. Thus, the effect of coconut oil consumption on the sirtuin signaling pathway should be further analyzed on the protein and protein function levels, especially in the context of using it to prevent neurodegenerative diseases.

During UPR activation, a decrease in RNA synthesis is often observed to protect cells against excessive accumulation of misfolded proteins [30]. According to our results, dietary beef tallow-activated genes (*DDX21*, *NOLC1*) suppress pre-rRNA transcription by forming ring-shaped structures surrounding multiple Pol I complexes [31,32]. Additionally, we observed upregulation of *NOB1*, which blocks the recruitment of mRNAs to the nascent ribosome [33]. Thus, it suggests the launching of repair mechanisms to stop the excessive production of misfolded proteins in the beef tallow group. On the other hand, several genes (*EIF3A*, *EIF5B*) engaged in translation initiation were upregulated in the beef tallow and coconut groups. It was recently found that eIF3a regulates HIF1α protein synthesis through the internal ribosomal entry site (IRES)-dependent translation. Therefore, it was concluded that eIF3a might be a potential therapeutic target for hepatic carcinoma (HCC) since it acts as a regulator for glycolysis—a process that is central to cancerous reprogramming of metabolism [34].

#### *4.4. UPR, Inflammation, and Cholesterol and Bile Acid Metabolism are the Main Processes A*ff*ected by Dietary Fatty Acids*

In our experiment, we observed that UPR, inflammation, and cholesterol and bile acid metabolism are the most altered processes under different dietary fat treatments. Interaction between endoplasmic reticulum stress—a trigger for UPR and inflammation—is involved in a variety of human pathologies [22,35]. Our results show that cholesterol and bile acid metabolism are additional components of this complex. They also support the hypothesis that a key point in these interactions is the composition of lipid membranes, since "organelle membrane" is the common cellular component enriched in all comparisons of our study (Supplementary Table S2). Considering the results of functional analysis of the identified DEGs as well as the qPCR results, it can be stated that the UPR level is low in the rapeseed oil group, while it is high in the beef tallow and coconut oil groups. The level of inflammation and cholesterol efflux is high only in the group receiving beef tallow. In contrast, cholesterol biosynthesis is low in the group receiving beef tallow and to a lesser extent in the coconut oil group. Our next step was the analysis of the causal relations between these processes and the identification of key genes regulating these interactions. The effect of bile acid and cholesterol metabolism on inflammatory processes has been well known for a long time. Numerous studies indicate that reduced flow of bile acids leads to their accumulation in the liver, which in turn causes inflammation. It was shown that bile acids act as an inflammagen, and directly activate signaling pathways in hepatocytes that stimulate the production of the proinflammatory mediators. However, other mechanisms were considered; the inflammatory response is triggered by activation of Toll-like receptor 4 (TLR4), either by bacterial lipopolysaccharide or by damage-associated molecular pattern molecules released from dead hepatocytes [36]. When we compared the effect of beef tallow and coconut oil in relation to rapeseed oil, we observed the downregulation of genes responsible for the synthesis and transport of bile acids (*CYP7A1*, *ABCB11*) only in the group receiving beef tallow. Similarly, only in this group, we observed an increase in expression of genes coding for acute-phase proteins. The exception here was the *ORM1* gene (identified as a hub gene by cytohubba software), which was upregulated in the coconut oil group. The product of this gene is classified as an acute phase protein, but it also has immunosuppressive activity. More than 25 years ago, it was found that this protein protects mice against lethality shock induced by tumor necrosis factor (TNF) or endotoxins [37]. More recently, a decrease of ORMDL protein without decreases in ORMDL mRNA levels was observed in HepG2 liver cells treated with the pro-inflammatory stimulus, and this observation was extended to in vivo models of inflammation [38]. In contrast, we observed a decrease in ORM expression accompanied by an increase in acute phase response in the beef tallow group at the mRNA level. The situation is further complicated by the fact that the *ORM1* gene is activated by bile acid through FXR—the nuclear bile acid receptor in mice [39]. This may indicate that the main difference in the action of coconut oil compared to beef tallow is that coconut oil does not reduce cholesterol catabolism and its disposal with bile acids. Thanks to this, FXR is activated, and thus ORM1′ s protective effect is maintained.

The interaction between inflammation and cholesterol and bile acid metabolism also acts in the opposite direction; the activation of immune response proteins affects the level of gene expression associated with cholesterol and bile acid metabolism [40]. In the liver, LPS markedly decreases the mRNA expression and activity of CYP7A1 and cholesterol transporters ABCG5 and ABCG8, which mediate cholesterol excretion into the bile in Syrian hamsters [41]. The relationship discussed here has a significant clinical implication. Until recently, hypercholesterolemia was considered the main cause of heart disease, but some scientists indicate that inflammation may play a more important role in the pathogenesis of CVD. Our results underline the significance of inflammation in this context. In essence, findings from epidemiological studies report low rates of cardiovascular disease among populations who consume coconut oil as part of their traditional diets (in India, the Philippines, and Polynesia) even though this fat is cholesterol-raising. Considering that in our study the animals receiving coconut oil did not show an increase in activation of acute-phase response, the limited influence of coconut oil consumption for the occurrence of heart disease may be due to the anti-inflammatory properties of this oil (Figure 9).

Inflammatory processes are strictly connected to endoplasmic reticulum stress, and UPR and the *ORM* genes link these two processes. Yeast cells lacking the *ORM1* and *ORM2* genes show a constitutive unfolded protein response, sensitivity to stress, and slow ER-to-Golgi transport [42]. ORM proteins function as the main regulators of the enzyme serine palmitoyltransferase—necessary for the process of sphingolipid biosynthesis, which is one of the main components of the cytoplasmic membranes of the endoplasmic reticulum, and disturbances in its homeostasis lead to ER stress [43]. The unfolded protein response can also alter sphingolipid metabolism. Bi-directional interactions between sphingolipids and the UPR have now been observed in a range of diseases, including cancer, diabetes, and liver disease [43]. Other studies have shown that ER stress can directly initiate the proinflammatory pathways, while proinflammatory agents such as ROS, TLR ligands, and cytokines can induce ER stress. As a result, activated UPR may further enhance the pro-inflammatory response [44].

As in the case of relationships discussed earlier, the connection between UPR and biosynthesis of cholesterol and bile acids is well documented and mutual. It was shown that within the ER, there are numerous membrane receptors detecting changes in cholesterol levels and cholesterol overload causes severe dysregulation of the ER [45]. To counteract these abnormalities, transcription factors (e.g., *SREBF2*, *LXR*) regulating cholesterol biosynthesis and immunological responses are triggered by the NRF2 protein [46]. In our experiment, we did not observe changes in the level of expression of genes coding for these molecules; however, they are largely regulated at the level of translation or post-translational modifications [47]. Among the DEGs identified in our report, the two genes *S1PR1* (upregulated in the beef tallow and coconut oil groups when compared to the rapeseed oil group) and *ORAI3* (downregulated only in the beef tallow group when compared to the rapeseed oil group) are potential links between cholesterol metabolism and ER stress and UPR. The protein encoded by *S1PR1* is structurally similar to G protein-coupled receptors and binds the ligand sphingosine-1-phosphate (S1P) with high affinity and high specificity. S1P, in turn, is a signaling sphingolipid acting as a bioactive lipid mediator. It is transported mainly by HDL and activates one of the S1PR1-mediated biological functions: calcium flux [48]. It was recently demonstrated that S1PR1-mediated calcium efflux is achieved through ORAI-membrane calcium channels. Thus, our data suggest that both the dietary fats beef tallow and coconut oil activate *S1PR1*, but only beef tallow downregulates *ORAI3* transcript expression. Interestingly, store-operated calcium (Ca2+) entry (SOCE) is mediated by *Orai3* only in breast cancer cells that express the estrogen receptor, contrary to estrogen receptor-negative cancer cells, which suggests a relationship between estrogen concentration and *ORAI3* expression [49]. In our experiment, probably as a consequence of reduced cholesterol catabolism in the beef tallow group, we observed a decrease in steroid hormone biosynthesis, which theoretically could lead to a decrease in *ORAI3* expression and dysfunctional calcium channels.

Due to the limited amount of space, we are not able to discuss all the interesting relationships observed in our experiment. To mention only a few: the "inositol metabolism", "urea cycle", and "glutathione metabolism" pathways deserve additional detailed analysis in terms of the effect of dietary fats on cancer and the development of metabolic, neurodegenerative, and cardiovascular diseases. Moreover, expression of a number genes—potential therapeutic targets in liver diseases—including *EGFR*, *GLUD1*, *GNL3*, and *RGS5*, has been shown to be modified by dietary fats and should be further investigated in this regard. Although our work provides a large amount of new information on the impact of consuming different sources of fat on gene expression in the liver, we are aware that these studies should be extended with analysis at the level of proteins and their function. Furthermore, analysis of bile acid species content in the intestine and liver and histological examination of these organs would give a full view of changes introduced by consumption of different sources of fat. Even though the material was collected from all animals from the same part of the liver, our samples most likely contained a mixture of different cells (hepatocytes, parenchymal cells, immune cells). Therefore, differences in expression observed between the groups of animals tested may partly result from the proportion of the content of these cells in the samples. To accurately assess what processes occur in the specific cells of the liver under the influence of dietary fat, experiments using single-cell methods are necessary.

### **5. Conclusions**

To sum up, our experiment showed that the type of consumed fat has a very significant impact on changes in gene expression in the liver of pigs. We observed that these changes are most intensively visible in three related processes: cholesterol and bile acid biosynthesis, UPR, and inflammation, playing a key role in the pathogenesis of civilization diseases. If one of these processes is dysregulated, repair mechanisms are triggered by activation of connected pathways. Therefore, when the body is excessively stimulated by improper nutrition, a vicious circle starts, in which the dysregulation of one process results in the dysregulation of the next. In our experiment, this situation most likely arose as a result of a diet containing beef tallow. In this group of animals, we observed deregulation of cholesterol and bile acid metabolism, activation of genes coding for acute phase proteins, and activation of the UPR when compared to animals fed with rapeseed oil. In contrast, in the coconut oil group, no activation of inflammation genes was observed, suggesting that some ingredients of coconut oil (lauric acid, polyphenols) can stop this vicious cycle and prevent the development of civilization diseases.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/9/1087/s1, Supplementary Table S1. List of DEGs identified in each comparison. Supplementary Table S2. List of DEGs and Gene ontology terms, KEGG and Reactome pathways, common for different comparisons. Supplementary Table S3. Results of the DEseq2 comparisons between males and females. Supplementary Table S4. List of Gene ontology terms, KEGG and Reactome pathways overrepresented in each comparison after STRING functional analysis. Supplementary Figure S1. MA and PCA plots after analysis of RNA-seq results with DESeq2 software. Supplementary Figure S2. Networks of canonical pathways identified after IPA software analysis in each comparison.

**Author Contributions:** Conceptualization, M.O. and M.S.; methodology, A.W., A.K., M.O.; software, T.S.; ´ validation, G.S., M.O.; formal analysis, M.O.; investigation M.O.; resources, M.S.; data curation, T.S.; ´ writing—Original draft preparation, M.O.; writing—Review and editing A.W., A.K., W.W., M.O., M.S. and T.S.; ´ visualization, M.O.; W.W., supervision, M.O.; project administration, M.O.; funding acquisition, M.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Science Centre, Poland (grant no 2014/13/B/NZ9/02134).

**Acknowledgments:** We thank Justyna Mrozowicz for technical assistance.

**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/).

*Article*

## **Identification of DNMT3B2 as the Predominant Isoform of DNMT3B in Porcine Alveolar Macrophages and Its Involvement in LPS-Stimulated TNF-**α **Expression**

### **Yanbing Zhang, Hui Li, Xiao Xiang, Yan Lu, Mona Sharma, Zongjie Li, Ke Liu, Jianchao Wei , Donghua Shao, Beibei Li, Zhiyong Ma and Yafeng Qiu \***

Shanghai Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Shanghai 200241, China; zhangyanbing129@outlook.com (Y.Z.); lihui022715@outlook.com (H.L.); xiao.xiang@wur.nl (X.X.); yanlu013194@outlook.com (Y.L.); monasharma1990@yahoo.com (M.S.); lizongjie@shvri.ac.cn (Z.L.); liuke@shvri.ac.cn (K.L.); jianchaowei@shvri.ac.cn (J.W.); shaodonghua@shvri.ac.cn (D.S.); lbb@shvri.ac.cn (B.L.); zhiyongma@shvri.ac.cn (Z.M.)

**\*** Correspondence: yafengq@shvri.ac.cn

Received: 20 August 2020; Accepted: 8 September 2020; Published: 10 September 2020

**Abstract:** DNA methyltransferase 3B (DNMT3B) as one member of the DNMT family functions as a de novo methyltransferase, characterized as more than 30 splice variants in humans and mice. However, the expression patterns of DNMT3B in pig as well as the biological function of porcine DNMT3B remain to be determined. In this study, we first examined the expression patterns of DNMT3B in porcine alveolar macrophages (PAM). We demonstrated that only DNMT3B2 and DNMT3B3 were the detectable isoforms in PAM. Furthermore, we revealed that DNTM3B2 was the predominant isoform in PAM. Next, in the model of LPS (lipopolysaccharide)-activated PAM, we showed that in comparison to the unstimulated PAM, (1) expression of DNTM3B is reduced; (2) the methylation level of TNF-α gene promoter is decreased. We further establish that DNMT3B2-mediated methylation of TNF-α gene promoter restricts induction of TNF-α in the LPS-stimulated PAM. In summary, these findings reveal that DNMT3B2 is the predominant isoform in PAM and its downregulation contributes to expression of TNF-α via hypomethylation of TNF-α gene promoter in the LPS-stimulated PAM.

**Keywords:** porcine alveolar macrophages; DNMT3B; DNA methylation; isoform; TNF-α

#### **1. Introduction**

DNA methyltransferase 3B (DNMT3B) is one member of the DNMT family which comprises DNMT1, DNMT3A, DNMT3B, as well as DNMT3L (DNMT3-like) in mammals. In mammalian systems, DNMT3B similar with DNMT3A serves as de novo methyltransferase for the establishment of DNA methylation [1]; in comparison, DNMT1 acts as maintenance of DNA methylation [2]. Moreover, DNMT3L has no DNA methyltransferase activity, but it can act as an accessory factor of the other DNMTs to regulate DNA methylation [3]. Although the functional characteristics between DNMT3A and DNMT3B are similar, the expression patterns of them are very different, i.e., in comparison to two isoforms of mouse DNMT3A [4], more than 30 isoforms of DNMT3B are identified in humans and mice [5,6].

Although there are more than 30 isoforms of DNMT3B, the expression patterns of DNMT3B appear to be highly conserved, at least in humans and mice. For example, human DNMT3B2 has a 60 bp-deletion (representative of exon 10) in comparison to the canonical isoform DNMT3B1 [5,7]; furthermore, DNMT3B3 of humans and mice has two deletions including a 60-bp deletion and a 189-bp deletion (representative of exon 21 and 22) in comparison to DNMT3B1 [5,8]. Notable, the alternative splicing of DNMT3B could influence DNMT3B functions. Thus, to clarify the expression patterns of DNMT3B is important to understand DNMT3B functions in the other species like pig.

Given the reported data, many spliced variants of DNMT3B are expressed in a tissue, cell, and/or developmental stage-specific manner [6,9], in this study, we focus on the expression pattern of DNMT3B in porcine alveolar macrophages (PAM), which are not only the major immune cells in pig lungs, but also the important resources of the inflammatory cytokines in pneumonia of pigs [10,11]. Though it has been shown that DNA methylation plays role in modulation of lipopolysaccharide (LPS)-induced inflammation in PAM [12], little is known about how DNMT3B is involved in this process. Furthermore, we determined effects of DNA methylation regulated by DNMT3B on TNF-α expression in the LPS-activated PAM. In this study, for the first time we determine the expression pattern of DNMT3B and clarify its effect on TNF-α expression in the LPS-stimulated PAM.

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

#### *2.1. Piglets and Porcine Alveolar Macrophages (PAM)*

Clinical healthy 35-day-old piglets (Shanghai great white pig strain) were purchased from the Shanghai Academy of Agricultural Sciences (Shanghai, China). All animal experiments were approved by the Institutional Animal Care and Use Committee of Shanghai Veterinary Research Institute (IACUC No: Shvri-po-201606 0501) and were performed in compliance with the Guidelines on the Humane Treatment of Laboratory Animals (Ministry of Science and Technology of the People's Republic of China, Policy No.2006 398). Porcine alveolar macrophages (PAM) were isolated from piglets as previously described [13] and cultured in RPMI 1640 containing 10% FBS, penicillin, streptomycin, GlutaMAX (all purchased from Thermo Fisher Scientific, Shanghai, China).

#### *2.2. Cloning of Porcine DNMT3B Isoforms and Sequence Analysis*

Total RNA was extracted from PAM (at least 1.0 × 10<sup>6</sup> cells) by using Trizol method [14] (Takara Biotechnology, Dalian, China), and cDNA was prepared by using Super ScriptII Reverse Transcriptase (Thermo Fisher Scientific, Shanghai, China). Three pairs of primers (Table S1) used to amplify full-length porcine DNTM3B were designed based on the porcine DNMT3B1 sequence deposited in GenBank (XM\_013985274.2). All PCR products were cloned using pMD19-T Vector Cloning Kit (Takara Biotechnology, Dalian, China), then positive clones were selected and sequenced by using M13 (Bacteriophage M13) forward and reverse primers. The sequences of DNMT3B2 and DNMT3B3 were obtained and deposited in NCBI GenBank with accession numbers MN873575 and MN207312, respectively. Subsequently, porcine DNMT3B2 or DNMT3B3 was inserted into the p3 × Flag-CMV-14 vector (Sigma, St. Louis, MO, USA) by homologous recombination with the ClonExpress MultiS One Step Cloning Kit (Vazyme Biotech, Nanjing, China) (Table S1), named Flag-DNMT3B2 or DNMT3B3, respectively.

Furthermore, the expression of exon 10 was analyzed by RT-PCR (reverse transcription-PCR) with an exon 9 forward primer and an exon 11 reverse primer (Table S1). Moreover, PCR was performed to detect the presence of exon 10 in genomic DNA of PAM with an intron 9 forward primer and an intron 10 reverse primer (Table S1). Furthermore, semi-quantitative RT-PCR was performed with an exon 20 forward primer and an exon 23 reverse prime (Table S1) to analyze the expression profile of DNMT3B2 and DNMT3B3 in PAM. All PCR products were sequenced using gene-specific primers. All of the images for agarose gel electrophoresis were captured by image lab version 5.1 (Bio-Rad Laboratories, Hercules, CA, USA).

The amino acid sequence alignment of human (NP\_008823.1), mouse (NP\_001003961.2), and three porcine DNMT3B isoforms was performed using the Clustal V method and edited using Genedoc. The phylogenetic tree was constructed using the available DNMT3B proteins by the neighbor-joining method with 1000 bootstrap replicates in MEGA version 6.06 [15,16]

#### *2.3. Generation of Polyclonal Antibodies against Porcine DNMT3B*

Rabbit experiments were approved by the Institutional Animal Care and Use Committee of Shanghai Veterinary Research Institute (IACUC No: Shvri-po-201606 0501) and were performed in compliance with the Guidelines on the Humane Treatment of Laboratory Animals (Ministry of Science and Technology of the People's Republic of China, Policy No.2006 398). Rabbit against porcine DNMT3B antibodies were generated according to a previous publication [17]. Briefly, a peptide of 15 amino acids (SYTQDLTGDGDGEGE) residues of the porcine DNMT3B was synthesized. Rabbits were immunized five times with the peptide in combination with complete or incomplete Freund's adjuvants every 14 days.After 7 days of the fifth immunization, rabbits were euthanized and serum was generated. The specificity of the generated DNMT3B antibodies was tested by Western blot.

#### *2.4. Western Blot*

The protein samples were prepared from the cell pellets as previously described [18,19]. Then the samples were separated on SDS-PAGE gel and transferred to NC (nitrocellulose, NC) membrane. After blocking in 5% nonfat milk, primary antibodies were added and incubated with membrane overnight at 4 ◦C (anti-Flag (1:5000, M2, Sigma), anti-DNMT3B (1:1000), anti-actin (1:10,000, clone C4, Sigma)). Secondary antibody was incubated for 1 h at room temperature (goat anti-mouse HRP (1:5000, Abcam), goat anti-rabbit (1:10,000, Abcam)). The images of Western blot were captured by image lab version 5.1 (Bio-Rad Laboratories, Hercules, California, USA).

#### *2.5. Quantitative Real-Time Reverse Transcription-PCR (qRT-PCR)*

cDNA was prepared using PrimeScript RT Reagent Kit including gDNA Eraser (Takara, Dalian, China). Then RT-PCR was conducted using a SYBR Premix Ex Taq kit (Takara, Dalian, China). Specific primers are shown in Table S1. The calculation was performed as described in our previous study [20].

#### *2.6. Bisulfite Sequencing PCR (BSP)*

Genomic DNA was extracted from PAM treated with LPS (1 µg/mL) and vehicle for 6 h. Then genomic DNA (0.8 µg) was subjected to bisulfite treatment using EZ DNA Methylation-Gold Kit (ZymoResearch, Los Angeles, CA, USA). Primers for Bisulfite sequencing PCR (BSP) (Table S1) were designed based on porcine *TNF*-α gene promoter sequence using online software (http: //www.urogene.org/cgi-bin/methprimer/methprimer.cgi) [21]. Subsequently, BSP [22] was performed using EpiMark Hot Start Taq DNA Polymerase (New England Biolabs, Ipswich, MA, USA) following the manufacturer's protocol. Briefly, using BSP primers amplified the region of TNF-α promoter, running an agarose gel to recover the PCR products. PCR products were cloned into the pMD19-T vector (Takara, Dalian, China). More than 10 positive clones were randomly selected for DNA sequencing [23,24]. The sequencing data and non-CpG-C-T conversion rates were analyzed using online QUMA software (http://quma.cdb.riken.jp/top/index.html) [25]. The total percentage of methylated CpG was calculated in each group including vehicle-treated, LPS-treated, vector-transfected, and DNMT3B2-transfected groups. Additionally, the difference of methylation level between certain groups was analyzed using Fisher's exact test of the online QUMA software.

#### *2.7. Lentivirus Production*

HEK293T cells were cultured in Dulbecco's modified Eagle's medium (DMEM) containing 10% FBS, penicillin, streptomycin (Thermo Fisher Scientific, Shanghai, China). The pLenO-DCE-DNMT3B2 or pLenO-DCE-Vector (Invabio, Shanghai, China) was co-transfected with pRsv-REV, pMDlg-pRRE, pMD2G (Addgene) into HEK293T cells using Lipofectamine 2000 reagent (Invitrogen, Carlsbad, CA, USA). The supernatants were collected at 72 h post-transfection and concentrated through ultra-centrifugation (25,000 rpm, 4 ◦C, 2 h, L7 Ultracentrifuge, Beckman, Duarete, CA, USA) after filtering through a 0.45 µm syringe filter [26,27].

### *2.8. Statistical Analysis*

All data shown are arithmetic means ± standard deviations. Statistical significance was assessed using unpaired Student's *t*-test by GraphPad Prism software version 5.01 (GraphPad Software, Inc., La Jolla, CA, USA). 'n' refers to the sample size.

### **3. Results**

#### *3.1. Identification of DNMT3B2 and DNMT3B3 as the Detectable Isoforms in Porcine Alveolar Macrophages (PAM)*

According to the predicted porcine *DNMT3B1* gene sequence (GenBank accession number: XM\_013985274.2), we designed the primers to amplify the DNMT3B ORF in PAM cDNA. Interestingly, the full-length sequencing results showed that only DNMT3B2 (GenBank accession number, MN873575) and DNMT3B3 (GenBank accession, MN207312) were identified in PAM (Figure 1A,B). Given that alternative splicing of DNMT3B exon 10 [5,7] distinguished DNMT3B1 (exon 10-included isoforms) with DNMT3B2 and DNMT3B3 (the exon 10-excluded isoforms), we further investigated the expression of DNMT3B exon 10 in PAM. In comparison to the expected fragment (about 160 bp) containing exon 10, a short fragment (about 100 bp) was obtained (Figure 1C). Moreover, the sequencing analysis confirmed that DNMT3B exon 10 was absent in the PCR product. We also investigated the presence of exon 10 in *DNMT3B* gene in PAM by PCR. Successfully, we obtained the fragment containing exon 10 (Figure 1D) which was further confirmed by the sequencing analysis. Taken together, these data reveal that expression of DNMT3B exon 10 is lost in PAM. Consistently, DNMT3B2 and DNMT3B3, the exon 10-excluded isoforms, are detectable in PAM.

**Figure 1.** *Cont.*

**Figure 1.** Identification of two splice variants of DNA methyltransferase 3B (DNMT3B) in porcine alveolar macrophages (PAM). (**A**) Schematic representation of the predicted DNMT3B1 and two splice variants identified in PAM. (**B**) Sequencing chromatogram showing the deleted region in the short splice variant named DNMT3B3 in comparison with the long splice variant named DNMT3B2. (**C**) Expression analysis of exon 10 of DNTM3B and GAPDH (about 90 bp) in PAM cDNA. Schematic representation of the exon 10-including fragment (about 160 bp) and the exon 10-excluded fragment (about 100 bp) from DNMT3B1 and DNMT3B2/3, respectively. Note that only a short fragment was obtained, indicating that expression of exon 10 of DNMT3B was absent in PAM. (**D**) As shown in the schematic, PCR was performed with an intron 9-forward primer and an intron 10-reverse primer to identify the presence of exon 10 (the positive fragment with 121 bp) in the genomic DNA from PAM. Note that a unique fragment was obtained, which was confirmed by sequencing analysis to show the presence of exon 10 in the *DNMT3B* gene of PAM.

#### *3.2. Identification of DNMT3B2 as the Predominant Isoform in Porcine Alveolar Macrophages*

In comparison to DNMT3B2, we observed a 189-bp deletion in DNTM3B3 (Figure 1A,B), which is attributed to lack expression of exon 21 and exon 22 according to the previous reports [5,8]. Based on this expression pattern between DNMT3B2 and DNMT3B3, we set up the RT-PCR with an exon 20 forward primer and an exon 23 reverse prime to analyze the expression profile of DNMT3B2 and DNMT3B3 in PAM. As expected, we totally obtained two fragments by analysis of the RNA samples extracted from PAM: the long fragment (268 bp) is representative of expression level of DNMT3B2; the short fragment (79 bp) is representative of expression level of DNMT3B3 (Figure 2A). Notable, the expression abundance of DNMT3B2 looked greater than that of DNMT3B3 according to the results by agarose gel electrophoresis. The density calculation further confirmed that expression of DNMT3B2 was much higher than that of DNMT3B3 in PAM (Figure 2B). Furthermore, in order to determine the protein level of DNMT3B2 and DNMT3B3, we used an antigen peptide of DNMT3B to generate polyclonal antibodies against DNMT3B. Western blot showed that the polyclonal antibodies were able to specifically detect expression of porcine DNMT3B2 and DNMT3B3 (Figure 2C). Then, we used the polyclonal antibodies to detect the expression profile of DNMT3B in PAM. Our result showed that only DNMT3B2 was detected in PAM (Figure 2D). Collectively, our data reveal that DNMT3B2 is the predominant isoform in PAM.

**Figure 2.** (**A**) Total RNA was extracted from PAM and analyzed by RT-PCR with exon 20 forward primer and an exon 23 reverse primer for the expression profile of DNMT3B2 and DNMT3B3 in PAM. Note that two fragments were obtained: one is more than 250 bp, representative expression of DNMT3B2; another is less than 100 bp, representative of DNMT3B3. (**B**) Furthermore, the bar graph represents relative expression of DNMT3B2 or DNMT3B3 by calculating the DNMT3B2/GAPDH ratio. The data shown are mean (SD) (*n* = 3), \*\*\*, *p* < 0.001. (**C**) HEK 293T cells were transiently transfected with recombinant plasmids expressing Flag-pDNMT3B2 and Flag-pDNMT3B3, respectively. Cell lysates were analyzed by Western blot with anti-*p*DNMT3B or anti-Flag antibodies. (**D**) The expression profile of DNMT3B in PAM was analyzed by Western blot. Note that only DNMT3B2 was detected in PAM. \* is representative of nonspecific reaction.

#### *3.3. Sequence Analysis of Porcine DNMT3B2 and DNMT3B3*

Having demonstrated that porcine DNMT3B2 and DNMT3B3 share the same expression pattern with that in humans and mice, we next performed sequence analysis to determine the sequence similarity and evolutional relationship of porcine DNMT3B with the other species. Multiple sequence alignment illustrated that the amino acid sequence and function domain of DNMT3B were conserved among human, mouse, and pig (Figure 3A). Notable, both PWWP domain (named as the well-conserved residues, Pro-Trp-Trp-Pro) and PHD domain (the plant homeodomain) in porcine DNMT3B are more than 94% identical to that of human and mouse (data not shown). Moreover, two deletion regions in the porcine DNMT3B3 were founded in comparison to DNMT3B1 of human, mouse and pig, which were attributed to alternative splicing of exon 10, exon 21, and exon 22. Furthermore, phylogenetical analysis of 18 protein sequences of DNMT3B classified them into three branches: mammals, birds, and fishes (Figure 3B), indicating that porcine DNMT3B is clustered with the other mammals.

**Figure 3.** (**A**) Comparative sequence analysis of DNMT3B protein from different hosts. The black and gray shading highlight identity in all of the selected sequences and similarity, respectively. The PWWP (yellow line) and PHD (blue line) domains are labeled in reference to the human DNMT3B1 (NP\_008823.1), respectively. Moreover, the antigen peptide is shown in the red box. (**B**) Phylogeny analysis of DNMT3B protein from different hosts. The scale bar indicates the genetic distance.

*α*

*α*

α

*α*

*α* α

α

*α*

*α* − −

#### *3.4. Downregulation of DNMT3B Associates with the Demethylation of TNF-*α *Gene Promoter in the LPS-Activated PAM*

Furthermore, these results by sequence analysis motivated us to determine the biological function of porcine DNMT3B. Given the reported data that expression of DNMT3B is downregulated in LPS-activated mouse macrophages, which is associated with the global demethylation of DNA [28], thus we chose the LPS-activated PAM as our model for the functional analysis of porcine DNMT3B. Since the alteration of DNMT3B in the LPS-stimulated porcine alveolar macrophages remains unknown, we first investigated the expression profile of DNMT3B in the LPS-activated PAM. In comparison to the untreated PAM, the mRNA level and protein level of DNMT3B was significantly downregulated in LPS-stimulation PAM. Besides, the expression profile of DNMT1 and DNMT3A was also determined in the untreated and LPS-activated PAM (Figure 4A,D). Interestingly, expression of DNMT1 and DNMT3A were not affected by LPS stimulation compared with those in the unstimulated PAM (Figure 4B,C). Thus, these data demonstrated that expression of DNMT3B is reduced in LPS-activated PAM, which might involve in the demethylation of DNA during LPS stimulation.

**Figure 4.** Total RNA was extracted from the treated and untreated PAM by LPS and expression of DNMT3B, DNMT1, and DNMT3A was analyzed by qRT-PCR analysis. The relative expression of DNMT3B (**A**), DNMT1 (**B**), and DNMT3A (**C**) was calculated as per material and method. *n* = 3, the data shown are mean (SD), ns is representative of no significant difference between control and LPS stimulated groups, \*, *p* < 0.05. (**D**) PAM with or without LPS treatment was harvested and analyzed by Western blot with anti-pDNMT3B antibody. Note that LPS stimulation decreased expression of DNMT3B2 in comparison to control group. The representative bar graph was calculated by DNMT3B2/actin ration. The data shown are mean (SD) (*n* = 3), \*, *p* < 0.05.

Given that the methylation status of *TNF-*α gene promoter plays a role in modulation of TNF-α expression [29–31], we next used the model of LPS-activated PAM to determine whether reduction of DNMT3B contributes to the demethylation of the *TNF-*α gene promoter. We performed Bisulfite sequencing PCR to examine the methylation profile of the *TNF-*α gene promoter region (−397 to −151) in the untreated and LPS-activated PAM, respectively. In comparison to the untreated PAM, the methylation level of *TNF-*α gene promoters in the LPS-activated PAM was significantly decreased (Figure 5A), which is accompanied with the robust induction of TNF-α in the stimulated PAM (Figure 5B). Thus, our data reveal that the demethylation of *TNF-*α gene promoter is associated with induction of TNF-α in the LPS-activated PAM.

*α α* α **Figure 5.** (**A**) The DNA methylation level of *TNF-*α promoter region was analyzed by Bisulfite sequencing PCR (BSP). Each line represents one individual sequenced clone. One circle represents one CpG site: the open circle shows unmethylated CpG; the black circle shows methylated CpG. Note that LPS stimulation significantly reduced the DNA methylation level of *TNF-*α promoter region in comparison to control groups. The difference between LPS group (*n* = 3) and control group (*n* = 3) was calculated by Fisher's exact test, \*\*, *p* < 0.01. (**B**) The relative expression of TNF-α in PAM with or without LPS treatment was analyzed by qRT-PCR. Note that LPS stimulation promoted TNF-α expression in comparison to control group. *n* = 3, the data shown are mean (SD), \*, *p* < 0.05.

α

#### *α α 3.5. DNMT3B2-Mediated Methylation of TNF-*α *Gene Promoter Restricts Induction of TNF-*α *in the LPS-Stimulated PAM*

*α α* α α *α* − − *α* Given that reduction of DNMT3B is associated with the demethylation of *TNF-*α gene promoter and that DNMT3B2 is the predominant isoform in PAM, therefore, our final goal was to determine the biological effect of DNMT3B2 on methylation of *TNF-*α gene promoter as well as expression of TNF-α in PAM. In order to do so, PAM were treated with lentivirus expressing porcine DNTMT3B2 and lentivirus vector (negative control, NC), respectively. After 24 h, the cells were stimulated with LPS and then harvested for the certain analysis. As shown in Figure 6A, in comparison to NC group, PAM treated with lentivirus expressing porcine DNMT3B2 showed a high level of DNMT3B even under the LPS-stimulated conditions. In contrast, the induced TNF-α was restricted by the increased DNMT3B2 in comparison to NC-treated PAM (Figure 6B). Furthermore, the methylation level of *TNF-*α gene promoter (−397 to −151) was analyzed in these two groups. Notable, PAM treated with lentivirus expressing porcine DNMT3B2 showed the higher methylation level of the *TNF-*α gene promoter than that in NC-treated PAM (Figure 6C). Thus, these results reveal that DNMT3B2-mediated methylation of *TNF-*α gene promoter modulates expression of TNF-α in the LPS-stimulated PAM.

*α* α

α *α* μ α *α α* **Figure 6.** Effect of overexpression of DNMT3B2 by lentivirus on TNF-α transcription and methylation profile of *TNF-*α promoter. PAM were infected with lentivirus overexpressing DNMT3B2 or containing control vector. After 48 h of infection, cells were stimulated with LPS (1 µg/mL) for 6 h. Subsequently, cells were harvested, mRNA was extracted and analyzed for expression of DNMT3B2 (**A**) and TNF-α (**B**) by qRT-PCR, respectively. Data are shown as mean (SD) (*n* = 3). \*, *p* < 0.05; \*\*, *p* < 0.01. In the meanwhile, cells were treated as shown above, and the total DNA was collected and analyzed for the DNA methylation level of the *TNF-*α promoter region by BSP (**C**). Each line represents one individual sequenced clone. One circle represents one CpG site: the open circle shows unmethylated CpG; the black circle shows methylated CpG. Note that overexpression of DNMT3B2 significantly increased the DNA methylation level of the *TNF-*α promoter region in comparison to control group. The difference between DNMT3B2 group (*n* = 3) and control group (*n* = 3) was calculated by Fisher's exact test, \*\*\*, *p* < 0.001.

#### **4. Discussion**

α This study provides novel information regarding the expression pattern of DNMT3B in porcine alveolar macrophages and its effect on TNF-α expression in the LPS-activated PAM. Although there are more than 30 isoforms of DNMT3B identified in humans and mice, yet the expression pattern of DNMT3B in pigs remains to be determined. According to the predicted porcine DNMT3B1 gene sequence (GenBank accession number: XM\_013985274.2), we sought to clone the different isoforms of DNMT3B in porcine alveolar macrophages (PAM). In this study, we demonstrate that DNMT3B2 and DNMT3B3 are the detectable isoforms in PAM. In contrast, the canonical isoform, DNMT3B1 and the other splice variants were not founded in the analysis. Since expression of DNMT3B exon 10 was not detected in PAM, alternative splicing of DNMT3B exon 10 [5,7] is important to result in the exon 10-included isoforms (DNMT3B1) and the exon 10-excluded isoforms (DNMT3B2 and DNMT3B3), respectively. Furthermore, we identified DNTM3B2 as the predominant isoform in PAM.

In the model of LPS-activated PAM, we show that in comparison to the unstimulated PAM, (1) expression of DNTM3B is reduced; (2) the methylation level of *TNF-*α gene promoter is decreased. We further establish that DNMT3B2-mediated methylation of *TNF-*α gene promoter restricts induction of TNF-αin the LPS-stimulated PAM. Collectively, our data identify DNMT3B2 as the predominant isoform to have a role in modulation of expression of TNF-α via DNA methylation in the LPS-stimulated PAM.

It has been reported that more than 30 DNMT3B isoforms are identified in humans and mice [5,6]. However, the expression patterns of porcine DNMT3B remains unknown. Here, we demonstrate that DNMT3B2 and DNMT3B3 are detectable in PAM; in comparison, DNMT3B2 is the predominant isoform. Furthermore, we demonstrate that the expression pattern of porcine DNMT3B2 and DNMT3B3 is same as what is shown in humans and mice [5,7,8], i.e., DNMT3B2 is lacking exon10; DNMT3B3 is lacking exon 10, exon 21, and exon 22. Thus, these results provide the evidence that the expression pattern of DNMT3B is conserved in different species [32].

Interestingly, we could not detect expression of DNMT3B1 in PAM despite the presence of exon 10 in the genome. In fact, many spliced variants of DNMT3B are expressed in a tissue, cell, and/or developmental stage-specific manner [6,9]. For example, DNMT3B1 is highly expressed in human ES cells; in contrast, its expression is decreased in the somatic cells which is accompanied with the predominant expression of DNMT3B3 [33]. Therefore, it is possible that porcine DNMT3B1 might be highly expressed in certain cells but not PAM during the development of pigs.

DNMT3B as one member of DNMTs catalytically regulates de novo DNA methylation. However, not all of DNMT3B isoforms have active catalytic domains. Notable, DNM3B1 and DNMT3B2 have the active catalytic domains; in comparison, DNMT3B3 is considered as an inactive catalytic isoform with some deletions in the catalytic domain. Although the role of DNMT3B3 in catalytic activity is controversial, the related studies have shown that DNMT3B3 could regulate DNA methylation by acting as the accessory protein. In this study, we demonstrated porcine DNTM3B2 regulates methylation of the *TNF-*α gene promoter, which has a role in modulation of expression of TNF-α. Of course, it is not excluded that porcine DNMT3B3 still functions as a positive regulator of DNA methylation. Thus, future studies are needed to answer this question.

#### **5. Conclusions**

In summary, this study for the first time identified DNMT3B2 as the predominant isoform in PAM. We demonstrated an important role of DNMT3B2-mediated DNA methylation in expression of TNF-α in the LPS-stimulated PAM. In this context, this study could provide evidence that DNMT3B2-mediated DNA methylation is the important mechanism responsible for understanding the inflammatory response in the LPS-stimulated PAM, even possible during some bacterial infections in the lungs of pigs.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/9/1065/s1, Table S1: Sequence of primers.

**Author Contributions:** Y.Z.: Conceptualization, Methodology, Writing—original daft, Visualization. H.L., X.X., and Y.L.: Conceptualization, Visualization, Methodology. M.S., Z.L., K.L., J.W., D.S., and B.L.: Methodology, Visualization. Z.M.: Supervision. Y.Q.: Conceptualization, Supervision, Writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was in part supported by the national key R&D program of China (2018YFD0500101), the National Natural Science Foundation of China (31,972,693), Central Public-interest Scientific Institution Basal Research Fund (No. 2020JB04), and Elite program of CAAS (to Y.Q.).

**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/).

*Article*
