**Insight into the Possible Formation Mechanism of the Intersex Phenotype of Lanzhou Fat-Tailed Sheep Using Whole-Genome Resequencing**

#### **Jie Li 1,**†**, Han Xu 1,2,**†**, Xinfeng Liu 1, Hongwei Xu 3,4,\*, Yong Cai 5,6 and Xianyong Lan 1,\***


Received: 7 May 2020; Accepted: 25 May 2020; Published: 29 May 2020

**Simple Summary:** Individuals with hermaphroditism are a serious hazard to animal husbandry and production due to their abnormal fertility. The molecular mechanism of sheep intersex formation was unclear. This study was the first to locate the homologous sequence of the goat polled intersex syndrome (PIS) region in sheep and found that the intersex traits of Lanzhou fat-tailed sheep were not caused by the lack of this region. By detecting the selective sweep regions, vital genes associated with androgen biosynthesis and the follicle stimulating hormone response entry were found, including *steroid 5 alpha-reductase 2* (*SRD5A2*), and *pro-apoptotic WT1 regulator (PAWR*). Additionally, the copy number variations of the four regions on chr9, chr1, chr4, and chr16 may affect the expression of the gonadal development genes, *zinc finger protein*, *FOG family member 2* (*ZFPM2*), *LIM homeobox 8* (*LHX8*), *inner mitochondrial membrane peptidase subunit 2* (*IMMP2L*) and *slit guidance ligand 3* (*SLIT3*), respectively, and further affect the formation of intersex traits.

**Abstract:** Intersex, also known as hermaphroditism, is a serious hazard to animal husbandry and production. The mechanism of ovine intersex formation is not clear. Therefore, genome-wide resequencing on the only two intersex and two normal Lanzhou fat-tailed (LFT) sheep, an excellent but endangered Chinese indigenous sheep breed, was performed. Herein, the deletion of homologous sequences of the goat polled intersex syndrome (PIS) region (8787 bp, 247747059–247755846) on chromosome 1 of the LFT sheep was not the cause of the ovine intersex trait. By detecting the selective sweep regions, we found that the genes related to androgen biosynthesis and follicle stimulating hormone response items, such as *steroid 5 alpha-reductase 2* (*SRD5A2*), *steroid 5 alpha-reductase 3* (*SRD5A3*), and *pro-apoptotic WT1 regulator (PAWR*), may be involved in the formation of intersex traits. Furthermore, the copy number variations of the four regions, chr9: 71660801–71662800, chr1: 50776001–50778000, chr4: 58119201–58121600, and chr16: 778801–780800, may affect the expression of the *zinc finger protein*, *FOG family member 2* (*ZFPM2*), *LIM homeobox 8* (*LHX8*), *inner mitochondrial membrane peptidase subunit 2* (*IMMP2L*) and *slit guidance ligand 3* (*SLIT3*) genes, respectively, which contribute to the appearance of intersex traits. These results may supply a theoretical basis for the timely detection and elimination of intersex individuals in sheep, which could accelerate the healthy development of animal husbandry.

**Keywords:** sheep; intersex; whole-genome resequencing; copy number variation; forming mechanism

#### **1. Introduction**

Intersexuality, also known as hermaphroditism, refers to the phenomenon that a dioecious animal is characterized by female-to-male sex reversal or abnormal gonad development. Intersex individuals are unable to reproduce, which poses certain obstacles to the protection and breeding of endangered species, and causes production loss to animal husbandry. In vertebrates, the mechanisms of sex determination are mainly divided into two types, genetic sex determination and environmental sex determination [1]. Intersex mostly occurs in the goat population, with high occurrence frequencies (about 3%–10%) [2], while related reports on horses, donkeys, pigs, and sheep are few. In goats, it is also named polled intersex syndrome (PIS) for the phenomenon of intersex individuals often found in hornless goat populations [2]. In the previous study, a 11.7 kb deletion fragment containing a repeat sequence (AF404302) was cloned by PCR, and the complete absence of this fragment resulted in goat polled syndrome [3,4]. A recent study about long-read whole-genome sequencing of a PIS-affected goat and a horned control goat revealed the presence of a more complex structural variant consisting of a 10,159 bp deletion and an inversely inserted 480 kb duplicated segment containing two genes, *potassium inwardly rectifying channel subfamily J member 15* (*KCNJ15*) and *ETS transcription factor ERG* (*ERG*) [5]. The deletion of the PIS region was identified affecting the development of germ cell support cells, and it can also affect the expression of genes, including *Forkhead box L2* (*FOXL2*), *PIS-regulated transcript 1* (*PISRT1*), and *promoter FOXL2 inverse complementary* (*PFOXic*) [6], indicating that the lack of the PIS region is closely related to goat intersex traits.

Compared to goats, reports of intersex sheep are quite rare. Domestic sheep and domestic goats, diverging about 4 to 5 million years ago and evolving into two different branches, are relatively close in genetic distance, and they have many similar genetic targets during domestication [7]. Therefore, it was suspected that the cause of sexual traits in sheep may be similar to that of goat sex. The location of homologous sequences of goat PIS regions should be detected in sheep.

Genetic variations or regulatory regions may affect an individual's phenotypic traits by affecting the transcription or translation of key genes [8]. In bovines, a 110 kb deletion in the *MER1 repeat containing imprinted transcript 1* (*MIMT1*) gene was a prominent cause of bovine abortion and stillbirth [9]. Bovine osteosclerosis may be associated with a deletion of approximately 2.8 kb in exon 2 and part of exon 3 of the *solute carrier family 4 member 2* (*SLC4A2*) gene encoding an anion exchanger [10]. An increased copy number of the *prolactin receptor* (*PRLR*) and *sperm flagellar 2* (*SPEF2*) genes in the K locus on the Z chromosome in chickens were closely related to the slow feathering trait of the chicken [11,12]. Additionally, the copy number variation of the *sperm flagellar 2* (*ASIP*) gene in sheep was closely related to the coat color [13]. It was speculated that some mutations may lead to the abnormal expression of certain genes that affect sex formation, leading to the generation of intersex traits.

As an excellent sheep variety of meat and wool, Lanzhou fat-tailed (LFT) sheep are famous for their large and fat tail and a number of excellent characteristics, such as the crude feed tolerance, higher disease resistance, and higher resilience than other domestic sheep. In recent years, only two surviving intersexual sheep were found in the LFT sheep population. In this study, the genome-wide resequencing of two intersex LFT sheep and two normal LFT sheep was performed, and combined with the resequencing data of four normal Tan sheep to find potential genes or regions related with the formation of ovine intersex traits, and thus provide an important basis for the timely detection and elimination of intersex individuals in sheep conservation and expansion.

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

#### *2.1. Ethics Statement*

All implemented experiments were approved by the Institutional Animal Care and Use Committee and were in strict accordance with good animal practices as defined by the Northwest A&F University (protocol number NWAFAC1008).

#### *2.2. Animal and Sequencing*

In this study, the ear tissues of four sheep, including two surviving intersex (LZ1 and LZ2) and six normal (LZ3 to LZ8) Lanzhou fat-tail sheep from Lanzhou City were collected and stored in 70% alcohol at −80 ◦C. Genomic DNA was extracted from sheep ear tissues using the phenol-chloroform method according to a previously reported protocol [14]. The DNA samples were quantified using a Nanodrop 1000 (Thermo Scientific, Waltham, MA, USA). Four DNA libraries, two intersex (LZ1 and LZ2) and two normal LFT sheep (LZ3 and LZ4), with insert sizes of approximately 350 bp, were constructed following the manufacturer's instructions, and 150 bp paired-end reads were generated using the Illumina HiSeq X10 platform.

#### *2.3. Sequence Quality Checking and Mapping*

Considering the relatively close genetic relationship between LFT sheep and Tan sheep, the previous four resequencing data of Tan sheep were also used as a normal control group for further analysis [15,16]. Before alignment, the FastQC software was used (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) to detect the joint information, the length information, and the quality information for each base on each read of all the raw data. Based on the results of the above quality control, the adapter and low quality raw paired reads were filtered using Trimmomatic (v0.36) [17]. The high-quality reads were mapped to the sheep version 4.0 reference genome (GCF\_000298735.2\_Oar\_v4.0\_genomic.fna) using the 'mem' algorithm of the Burrows–Wheeler Alignment Tool (BWA) software [18]. SAMtools (http://samtools.sourceforge.net/) was used to remove replicate sequences [19].

#### *2.4. Calling and Validation of SNPs and CNVs*

The calling of single nucleotide polymorphisms (SNPs) was performed with Genome Analysis Toolkit (GATK, version 2.4–9) UnifiedGenotyper [20], and was annotated by ANNOVAR [21]. The called SNPs with reads >4 and quality ≥20 were used for further analysis. CNVcaller software (https://github.com/JiangYuLab/CNVcaller) was used for the detection of whole-genome copy number variation (CNV) [22]. The specific steps were in accordance with a reported study [23], including segmenting the reference genome into a window of a specified size (800 bp) with a step size of 400 bp. Herein, Vst was used to measure the difference in the size of each copy number variation between different groups [24,25]. The mean log2 ratio across all probes falling within a specific CNV region was calculated. The variance of the means was for the entire set (Vt). The average variance within populations was then calculated (Vs) by taking the mean between populations (i.e., the V intersex sheep and Vnormal sheep). Vst values were finally calculated using the standard formula Vst = (Vt − Vs)/Vt.

#### *2.5. Regional Localization and Depth Statistics for Ovine Homologous Sequences of PIS*

The sequence of the PIS region was extracted from the goat reference genome (GCF\_001704415.1\_ASM170441v1\_genomic.fna) by using SAMtools and then aligned to the sheep reference genome using BLAT v. 36 × 1 software (BLAST-Like Alignment Tool) to obtain the ovine homologous regions [26].

After that, we used the SAMtools-depth to count the reads depth for each locus in the candidate region. Furthermore, we corrected it based on the total depth of sequencing to obtain the reads depth in the homologous regions of the eight samples, including two intersex LFT sheep, two LFT sheep, and four reported Tan sheep [15,16]. If the reads depth was essentially 0, the candidate region was completely missing on the ovine genome.

#### *2.6. PCR Amplification of PIS Candidate Region Sequences in Sheep*

In order to further verify whether the PIS candidate region of the intersex individual was really missing, primers, namely PIS-1, PIS-2, PIS-3, and PIS-4, for the homologous sequence in the ovine PIS region were designed to reference the second generation genome of the sheep chr1: 247747059–247755846 sequence (Supplemental Table S1). We performed assays of amplification on two intersex and six normal LFT sheep according to previous reaction volume and amplification procedure [27]. The products were separated on 2.5% agarose gels.

#### *2.7. Sweep Analysis of SNPs in Selected Regions*

In order to identify the selection signatures in the genomes of sheep, two sequenced pools based on genetic differentiation (Fst) of each 150 KB genome window were separately performed. The specific formula of Fst was consistent with the previous studies [28,29].

After Z-transformation of Fst, the candidate selection windows were used by selecting the top 1% in Fst score intersections [30]. Finally, the geneview module of python and R were utilized for data visualization.

#### *2.8. Gene Annotation and Functional Enrichment Analysis of Selected Signal Regions*

Gene functional enrichment analysis was performed primarily on the KOBAS 3.0 website (http://kobas.cbi.pku.edu.cn/index.php) [31]. Herein, considering that the sheep database in KOBAS was not complete, a Perl script was used to compare the longest protein sequences of each gene in sheep and human by invoking Blastp, and the gene extracted in the previous step (2.4 and 2.7) was converted into a gene homologous to "human", and then human was selected as a species for annotation.

#### **3. Results**

#### *3.1. Sequencing, Filtering, and Mapping*

Whole-genome sequencing of two normal Tan sheep, as well as two intersex sheep and four normal LFT sheep, was performed on an Illumina HiSeq X10 platform using genomic DNA and 198.5 Gb of high quality paired-end reads were generated. After getting the raw data, we utilized FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality testing. According to the test results, the clean reads were obtained using trimmomatic (http://www.usadellab.org/cms/index.php/ page=trimmomatic) to filter the low-quality reads (Supplemental Table S2). The clean reads of eight sheep were mapped to the sheep reference genome (GCF\_000298735.2\_Oar\_v4.0\_genomic.fna) using BWA-MEM [17,31], with an average mapping rate above 94% (Supplemental Table S3).

#### *3.2. Single Nucleotide Polymorphisms (SNPs) Calling and Annotation*

A total of 29,732,629 SNPs were obtained among eight sheep, among which, the most distributed SNPs were on chromosome 1 (n = 3,186,739), and the least distributed, lowest number of SNPs were on chromosome 24 (n = 459,707) (Supplemental Table S4). According to the results of the annotation, the proportions of transition (ts) mutations (A/G 10333679 and T/C 10308562) and transversion (tv) mutations (A/T 2050141, A/C 2457813, G/T 2450228, and G/C 2132206) were 69.4% and 30.6%, respectively, which met the 3 to 1 ratio (Figure 1).

**Figure 1.** The statistics of mutation type of single nucleotide polymorphisms (SNPs). Note: The proportion of translation (ts) mutations and transversion (tv) mutations meet the 3 to 1 ratio.

#### *3.3. The Localization Results of Homologous Sequences of Goat PIS Region in Sheep*

According to the specific information of the goat PIS region (AF404302.1: 27015–38775) in the NCBI, the obtained sequence of this 11.76 kb segment [2] was aligned to the sheep genome using BLAT, and then the score was calculated using the pslScore.pl script provided by the BLAT software to evaluate the results of the comparison.

The BLAT results show that there were 79 alignment regions, and the alignment region with the highest score was located in the region of the ovine chromosome1 (start to end position: 247747059–250146105, 2399.046 kb) (Supplemental Table S5). Then, the detailed analysis of the 79 region revealed that the first 58 of them were a closely connected 8787 bp region (chr1: 247747059–247755846), corresponding to a segment of the 11.7 kb region of the goat (AF404302.1:30003–38775) with 100% coverage (Figure 2a). Further localization revealed that this region was located approximately 340 kb upstream of the *FOXL2* gene (chr1: 248088730–248095868) in the sheep genome, which was consistent with the position of the goat *FOXL2* gene. Therefore, the 8787 bp sequence (chr1: 247747059–247755846) was served as a candidate region for ovine intersex traits.

**Figure 2.** *Cont*.

(**b**)

**Figure 2.** Study on the homologous sequences of the polled intersex syndrome (PIS) regions in sheep and goats. (**a**) The comparison results between the sheep chromosome 1: 247747059–250146105 and goat PIS area (the color line represents the matching relationship: goat 8773 bp matched sheep 8787 bp, and goat 2987 kb matched sheep 1813.4 kb). (**b**) The statistics of the read depth in the sheep 8787 bp candidate area. (**c**) Electrophoresis results of the four pairs of primer amplification products, indicating that the four regions of intersex sheep were not missing.

#### *3.4. Statistics of Reads Depth and PCR Amplification of Sheep 8787 bp Candidate Sequences*

In the two intersex sheep and normal individuals, the read depth was between 5 and 20, and there was no significant difference between them (Figure 2b).

According to the results of agarose gel electrophoresis, in two intersex LFT sheep (LZ1 and LZ2) and six normal LFT sheep (LZ3, LZ4, LZ5, LZ6, LZ7, and LZ8), four pairs of primers (PIS-1, PIS-2, PIS-3, and PIS-4) were able to amplify the target band in the candidate region, and the product lengths were 118, 376, 380, and 454 bp (Figure 2c), respectively, indicating that the candidate region of intersex sheep was not missing. It is further explained that the intersex traits of LFT sheep was not caused by the lack of the region.

#### *3.5. Genome-Wide Selection Sweeping Analysis in Intersex and Normal Populations*

*ZFst* score were calculated for the only two intersex LFT sheep and six normal sheep populations (including two LFT sheep and four Tan sheep). The top 1% percent of the windows with the highest *ZFst* score were defined as candidate selective sweep regions. The results show that most chromosomes contained windows with a higher differentiation coefficient, and chromosomes 2, 6, 7, 10, and 11 were strongly selected (Figure 3a). The region with the largest *ZFst* value was located on the X chromosome (chrX: 99675001–99825000, *ZFst* = 10.69), and the second one was located on chromosome 3 (chr3: 105000001–105150000, *ZFst* = 10.38) (Figure 3a; Supplemental Table S6).

**Figure 3.** The detection of genome-wide selection signals in intersex and normal populations. (**a**) The genome-wide distribution of *ZFst* between intersex sheep and normal sheep (including normal Lanzhou large-tailed sheep and Tan sheep). (**b**) The genome-wide distribution of *ZFst* between intersex and normal individuals of Lanzhou large-tailed sheep. Note: The blue line represents the top 1% of the ZFst value.

A total of 466 genes were obtained by gene annotation (Supplemental Table S7), and then those genes were functionally enriched by KOBAS3.0. A total of 1757 significant entries were found (*p* < 0.05), and 1038 significant entries remained after the false discovery rate (FDR)-correction (corrected *p*-Value < 0.05), mainly including items, such as muscle development, fat development, and immunity. Furthermore, the pathways were screened and two more significant entries related to female gonadal development (GO: 0008585) and the development of primary female sexual characteristics (GO: 0046545) were found (*p* < 0.05) (Supplemental Table S7). Nevertheless, neither entry was significant after correction.

#### *3.6. Genome-Wide Selection Sweeping Analysis in Intersex LFT Sheep and Normal LFT Sheep*

To exclude the effects of SNP locus frequencies between different breeds, the study then compared two intersex individuals and two normal LFT sheep. Consistent with the previous section, the region with the largest *ZFst* value was located on the X chromosome (Figure 3b). By contrast, the second region was also in the X chromosome, rather than chromosome 3 (Figure 3b). Moreover, a total of 451 genes were obtained by annotation (Supplemental Table S8). After functional enrichment analysis, 1680 significant entries were found (*p* < 0.05) and 969 significant terms remained after false discovery rate (FDR)-correction (corrected *p*-value < 0.05) (Supplemental Table S9). Screening of the pathway revealed five significant entries related to the synthesis and response of the sex hormone (*p* < 0.05). After correction (corrected *p*-value < 0.05), the genes, such as *SRD5A2*, *SRD5A3*, and *PAWR*, were significantly rich in the androgen biosynthesis process and their responses to follicle stimulating hormones (Supplemental Table S9). Additionally, the genes involved in androgen receptor signaling pathways, androgen metabolism, and the regulation of intracellular estrogen receptor signaling pathways, and gonadotropin response processes, such as *UFM1 specific ligase 1* (*UFL1*), *GTP-binding nuclear protein Ran-like* (*LOC105609617*), *mediator complex subunit 14* (*MED14*), *DEAD-box helicase 5* (*DDX5*), etc. (Supplemental Table S9).

#### *3.7. Detection of Genome-Wide Copy Number Variation (CNV) in Intersex and Normal Populations*

As copy number variation regions (CNVRs) can be separated by gaps or poorly assembled regions, the adjacent initial calls were merged if their reads depth were highly correlated. The default parameters were as follows: the distance between the two initial calls was less than 20% of their combined length, and the Pearson's correlation index of the two CNVRs was significant at the *p* = 0.01 level [21]. After calculation and combination through the CNVcaller software, 87,729 CNVRs were obtained, of which 1817 CNVRs were located on the scaffold sequence (not assembled into chromosomes) and 11,170 CNVRs were located on the X chromosome (Supplemental Table S10).

As the number of X chromosomes is different in females and males, the copy number variation on the scaffold and X chromosomes was not considered in the results. Therefore, 74,302 CNVRs located on autosomes were used for the subsequent analysis. As shown in Supplemental Table S10, the smallest proportion of CNVRs was chromosome 26 (5.8%), and the largest was chromosome 11 (9.92%). In addition, the number of CNVRs distributed on chromosome 1 (*n* = 8267) was the highest with a total length of 19,584,800 bp, while the smallest was on chromosome 26 (*n* = 1116) with 2,556,800 bp. The longest CNVR was 380,800 bp and located on chromosome 13 (Supplemental Table S10).

The copy number of the normal and intersex populations was screened with variation length >2000 bp and Vst value >0.25, and a total of 238 candidate regions were obtained. Gene annotation of the above 238 candidate CNVRs revealed that 140 were located in the intergenic region, and the remaining 98 overlapped with the genes (Supplemental Table S11). Functional enrichment analysis was performed on the annotated genes, revealing a total of 1838 significant entries (*p* < 0.05) and 1039 significant entries after FDR-correction (corrected *p*-value < 0.05) (Supplemental Table S12).

Through further GO enrichment analysis of CNVRs, four GO items related to female gonad development were found. Based on functional enrichment analysis, the enriched genes were mainly *ZFPM2, LHX8* (located at 43,507 nt downstream of its corresponding copy number region), IMMP2L (located at 400835 nt upstream of its corresponding copy number region), and *SLIT3* (Supplemental Table S12). The above four genes corresponded to four CNVRs, respectively. Furthermore, in the intersex individuals, the copy number region associated with the *LHX8*, *IMMP2L*, and *ZFPM2* genes were gained, while the copy number region associated with the *SLIT3* gene was lost (Table 1).


**Table 1.** The information of the four candidates copy number variation region and the genotype of eight samples.

Note: AA means two copies; AB means three copies; Ad means just one copy. *LHX8*, *LIM homeobox 8* gene; *IMMP2L*, *inner mitochondrial membrane peptidase subunit 2* gene; *ZFPM2*, *zinc finger protein FOG family member 2* gene; *SLIT3*, *slit guidance ligand 3* gene. The *V*st-values measure the difference in the size of each copy number variation between different groups.

#### **4. Discussion**

It is well-known that the mutation of the goat PIS region contributes to the absence of horns and sex-reversal [4]. Herein, through the reads depth and the PCR amplification experiments, the PIS homologous segments of two intersex sheep and normal sheep were not in a missing state, which indicated that the intersex trait of LFT sheep may not be caused by the absence of this region. The intersexuality of goats was always accompanied by hornlessness, manifested as non-interval syndrome, and the probability of the intersex appearance in the hornless goat population was 3% to 10% [32].

Genetic variations or regulatory regions may affect an individual's phenotypic traits by affecting the transcription or transition of key genes. We speculated that some mutations may lead to the abnormal expression of certain genes that affect sex formation, leading to the appearance of intersex. Furthermore, a previous study reported that the replication of the *SRY-box transcription factor 3* (*SOX3*) gene may also cause human developmental disorders [33]. An increased copy number of the *SRY-box transcription factor 9* (*SOX9*) gene may trigger the probability of intersex individuals [34], so it is suspected that the intersexuality of the sheep may be due to other reasons, such as genetic variations.

Functional enrichment analysis revealed significant entries for the androgen biosynthesis processes and follicle stimulating hormone responses, including *SRD5A2*, *SRD5A3*, and *PAWR* genes. Among them, mutations of the *SRD5A2*, which were closely related to testicular decline, could affect the formation of the urethra and external genitalia, leading to hypoplasia of the male reproductive organs [35]. Therefore, polymorphisms of the *SRD5A2* gene may be a key point leading to intersexuality in sheep.

Functional enrichment of the sequencing data for intersex sheep and normal sheep populations (including Tan sheep) did not find significant entries related to gonadal development or sex hormone metabolisms. A previous study reported that testicular tissue dysplasia was a key cause of human gender developmental disorders [35]. Additionally, the excessive synthesis of estrogen in vivo also affected the formation and development of female reproductive organs. Therefore, it was speculated that various mutations of genes related to the synthesis and secretion of androgen led to the fact that the testis and its accessory reproductive organs could not be maintained.

As hermaphroditism is not conducive to animal reproduction, it is speculated that intersexual individuals are selected to be eliminated in evolution, so the occurrence in the current group is low. A major limitation of this study is that there were only two intersexual individuals sequenced, which may result in a certain number of false positives in the sequencing results. If the limitation of the sample size is eliminated, the study may obtain more mutation regions or copy number variation regions that are more reliable than in this paper, and may even lock in the major genes that lead to intersex traits. It is undeniable that this study provides a direction and reference basis for further in-depth exploration of the molecular mechanism of sheep intersex traits.

#### **5. Conclusions**

This study was the first to locate the homologous sequence of the goat PIS region in sheep and found that the intersex traits of LFT sheep were not caused by the lack of this region. Through detecting the selective sweep regions, the vital genes associated with androgen biosynthesis and the follicle stimulating hormone response entry were found, including *SRD5A2*, *SRD5A3*, and *PAWR*. Additionally, the copy number variations of the four regions on chr9, chr1, chr4, and chr16 may affect the expression of the gonadal development genes, *ZFPM2*, *LHX8*, *IMMP2L*, and *SLIT3*, respectively, which contribute to the appearance of intersex traits.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/6/ 944/s1, Table S1: PCR amplification primer information for candidate region on sheep chromosome 1; Table S2: The statistics of resequencing samples from 4 Lanzhou large-tailed sheep; Table S3: The statistics of resequencing samples from 8 sheep; Table S4: Chromosome distribution of SNPs in sheep; Table S5: Goat PIS regional positioning on sheep genome (the top 10 score); Table S6: List of genes in the selected overlapping regions by top 1% highest Z(Fst) between intersex sheep and normal sheep (including normal Lanzhou large-tailed sheep and Tan sheep); Table S7: GO term of the candidate genes associated with the intersex; Table S8: List of genes in the selected overlapping regions by top 1% highest Z(Fst) between intersex and normal individuals of Lanzhou large-tailed sheep; Table S9: GO term of the candidate genes associated with the intersex; Table S10: Autosomal distribution of CNVRs in sheep; Table S11: Gene annotation of the candidate CNVRs; Table S12: Functional enrichment analysis of annotated genes.

**Author Contributions:** Methodology, H.X.; software, H.X. and X.L. (Xinfeng Liu); validation, H.X.; formal analysis, H.X. and X.L. (Xinfeng Liu); resources, H.X. and Y.C.; writing—original draft preparation, J.L.; writing—review and editing, X.L. (Xianyong Lan); supervision, X.L. (Xianyong Lan) and H.X.; funding acquisition, X.L. (Xianyong Lan) and H.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 31660642; 31360529; 31760649), Fundamentai Research Funds for the Central Universities (31920190020, 31920190004) and Science-technology Support Plan Project of Gansu Province (18JR3RA373,18YF1FA121), The Program for Changjiang Scholars and Innovative Research Team in the University (IRT\_17R88).

**Acknowledgments:** We greatly thank the staffs of Ruilin Sci-Tech Cluture and Breeding Limit Company Yongjing county, Gansu Province, for collecting samples of Tan sheep and Lanzhou fat-tail sheep. Moreover, we sincerely thank Wang XL and Jiang Y from Northwest A&F University for their data and technical support.

**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* **Influence of Age and Immunostimulation on the Level of Toll-Like Receptor Gene (***TLR3***,** *4***, and** *7***) Expression in Foals**

#### **Anna Migdał 1,\*, Łukasz Migdał 1, Maria Oczkowicz 2, Adam Okólski <sup>3</sup> and Anna Chełmo ´nska-Soyta 4,5**


Received: 19 August 2020; Accepted: 20 October 2020; Published: 26 October 2020

**Simple Summary:** Detailed knowledge of the molecular mechanisms of immunoglobulin synthesis appears necessary for a better understanding of foal immunity maturity and its influencing factors. At the same time, it encourages studies regarding the influence of the signaling cascade's proteins on the primary immunological response, which provides an opportunity to develop extremely precise methods of regulating acquired immunity. The results revealed that the expression of the*TLR3* and *TLR4* genes, as well as the levels of immunoglobulins and interleukins, can be modulated by stimulation with the pharmacological agent, and that the expression of the *TLR3* and *TLR4*genes in peripheral blood cells is dependent on age.

**Abstract:** The aim of this study was to investigate the molecular mechanisms leading to the identification of pathogens by congenital immune receptors in foals up to 60 days of age. The study was conducted on 16 foal Polish Pony Horses (Polish Konik) divided into two study groups: control (*n* = 9) and experimental (*n* = 7). Foals from the experimental group received an intramuscular duplicate injection of 5 mL of Biotropina (Biowet) at 35 and 40 days of age. The RNA isolated from venous blood was used to evaluate the expression of the*TLR3*, *TLR4*, and *TLR7* genes using RT-PCR. The results of the experiment demonstrated a statistically significant increase in the level of *TLR3* gene expression and a decrease in the level of*TLR4* gene expression with foal aging. The level of *TLR7* gene expression did not show age dependence. Immunostimulation with Biotropina had a significant impact on the level of the genes' expression for Toll-like receptors. It increased the level of *TLR4* expression and decreased *TLR3* expression. Thus, it was concluded that the expression of the*TLR3* and *TLR4*genes in peripheral blood cells is dependent on age. This experiment demonstrated a strong negative correlation between *TLR3* and *TLR4* gene expression.

**Keywords:** *TLR3*; *TLR4*; *TLR7*; foals; immunostimulation; gene expression

#### **1. Introduction**

Immune response differs in newborn and adult horses. Despite involving similar components, the regulation of immunity and the response to antigens vary. Foals are born with small, non-protective amounts of endogenous serum immunoglobulins (i.e., IgM and IgG) [1]. Immunoglobulin transport is limited due to the structure of the horse placenta (placenta spuria). That is why the suckling of colostrum is essential in the first hours of a foal's life. Immunological outcomes in newborn foals differ as compared to adults and are distinguished by modified cytokine profiles, as well as reduced antibody and T-cell responses [2]. Moreover, foals have a very low level of immunoglobulins in their blood plasma. An innate immune system composed of pre-existing or rapidly induced defenses is critical for newborn foals, while an antigen-specific response requires exposure to pathogens and time for development after birth [3]. The effectiveness of the immunological reaction is controlled by a complexity of direct and indirect mechanisms involving interactions among various cells and cytokine-induced actions. Many different immune cells are involved in maintaining a balanced immune response [4]. Receptors called pathogen recognition receptors (PRRs) represent very important elements found in immune cells. Thanks to their conservative structure, these receptors can recognize signals associated with a variety of pathogens. Pathogen-associated molecular patterns (PAMPs), PRR ligands, are molecules specific to viruses, bacteria, and other microorganisms with evolutionarily conserved structures [5]. Toll-like receptors (TLRs) are the most closely investigated PRRs and are one of the most essential components of immune responses [6,7]. Toll-like receptors play a key role in activating and stimulating innate as well as acquired immunity [8]. Identification of the threat and activation of TLRs triggers an immune response, leading to the elimination of this threat from the organism, which involves two basic reactions, namely, inflammatory, and antiviral. Cells with activated receptors release large amounts of proinflammatory cytokines, chemokines, and defensins, and these released factors initiate the migration and aggregation of immune cells (e.g., leukocytes, macrophages, mast cells, and dendritic cells) at the site of the pathogen invasion [7]. Activated TLRs present on the surface of macrophages lead to increased synthesis of the proinflammatory cytokines IL-1, -6, -8, -12, and TNF-α. In addition, complexes of ligands and TLR4 receptors increase the phagocytic activity of macrophages and stimulate the production of reactive oxygen species (ROIs) and the synthesis of nitric oxide (NO). TLR-activated macrophages enhance the expression of major histocompatibility complexes I and II (MHCI and MHCII), CD80, CD86, and co-stimulators that make immune cells more efficient in displaying T-cell antigens that induce specific immune responses [9,10].

Approximately 20% of foals die before the end of the second month, which brings significant economic and breeding losses [11,12]. This justifies undertaking research into new and novel techniques for stimulation of the immune system. Gaining insight into the behavior of TLRs seems to be necessary for expanding our understanding of the mechanism responsible for the development of foal resistance/immunity, and the identification of determinants for further building it up. In addition, it is also an important element contributing to finding innovative solutions in the fight against infections and new ways to improve the prevention and treatment of infections in animals. The paradox of neonatal vaccination is the need of immediate protection during early days, the perceived limitations of the immune system of neonate foals, and the theory of maternal antibody interference [13]. Studies have shown that the immune system of neonatal foals is also naive and immature relative to juvenile and adult horses [14–17]. Several studies have suggested that basal TLR expression in full-term neonatal blood monocytes is similar to that of adults [18,19]. The TLR-mediated production of cytokines by neonatal monocytes, however, is very different in newborns compared to that of adults [19]. Thus far, little is known about the development of the horse immune system during pre- and postnatal periods, which negatively affects the ability to devise strategies for maintaining and improving foal health. Based on its biological properties, as well as the influence of Toll-like receptors on the immune response traits of farm animals and humans, we hypothesized that gene expression for Toll-like receptors TLR-3, TLR-4, and TLR-7 in foals is dependent on factors such as age and immunostimulation. The aim of this work was to investigate the molecular mechanisms leading to the identification of pathogens by

congenital immune receptors in foals up to 60 days of age, including the verification of the hypothesis concerning age-related expression of the*TLR3*, *TLR4*, and *TLR7* genes.

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

This experiment was granted permission from the Local Ethics Committee in Kraków (no 37, 30 May 2016).

#### *2.1. Animals and Feeding*

Studies were carried out on 16 foals representing Polish Pony horses (Polish Konik). This primitive horse breed is genetically and phenotypically closely related to its wild ancestor, the Tarpan Horse (Eurasian wild horse) [20].

All foals with mares were kept in the same stable in individual boxes (size 2.15 × 3.50 m) on permanent straw bedding at the Experimental Station of the University of Agriculture in Krakow. All animals were clinically healthy throughout the experimental period. Mares of 5–17 years of age and 270–340 kg live body weight were not vaccinated during pregnancy. Foal birth weight was 27–35 kg, and weight loss on the first day of life was <1.5%. The horses had all been used by university students in the teaching program. No horses were used for equestrian purposes. Inclusion criteria consisted of foals born from healthy mares with no placentitis, a normal gestational period, an uneventful birth, and normal physical and neurological examination findings. The foals had to successfully stand and nurse within 2 h of birth and remain clinically healthy during the study period.

Mares were fed ad libitum with hay (*Lolium* 40% and *Trifolium L.* 20%) with the addition of oats in the amount of 1.5 kg/mare/day [21]. Foals were fed only with colostrum and mother's milk ad libitum, without additional supplementation. Water was offered from automatic water drinkers (flow ~ 10 L/min).

#### *2.2. Experimental Design*

Two weeks before delivery, birth alarms (Abfohlsystem, Jan Wolters, Steinfeld, Germany) were placed in the labia, and mares were moved to box stalls inside a stable lit with natural light (Figure S1). During the experiment, foals were kept with their mothers in individual boxes, and when leaving the stalls with their mothers for the pasture, they were randomly assigned into the following groups:


For the immunostimulation, a commercially available immunostimulator was used in the present study, namely, Biotropine (Biowet Drwalew S.A., Drwalew, Poland), which consists of a mixture of inactivated Gram-positive bacteria, e.g., *Staphylococcus aureus* (74 mg/mL), *Streptococcus zooepidermicus* (24.6 mg/mL), *Streptococcus equi* (24.6 mg/mL), *Streptococcus equisimilis* (24.6 mg/mL), *Streptococcus agalactiae* (24.6 mg/mL), *Streptococcus dysgalactiae* (24.6 mg/mL), *Erysipelothrix insidiosa* (49 mg/mL),and Gram-negative bacteria, e.g., *Escherichia coli* (123 mg/mL) and *Pasteurella multocida* (123 mg/mL) as well as pork spleen extract (10 mg/mL). On days 35 and 40 after birth, the foals from the experimental group received an intramuscular (*m. pectoralis descendens*) injection of 5 mL of Biotropine.

#### *2.3. Blood Sampling and Blood Analysis*

Blood samples were collected from foals by jugular venipuncture. Blood samples were obtained from foals up until 60 days of age according to the following scheme: After birth before the first suckling and then on the 1st, 3rd, 5th, 10th, 20th, 30th, 40th, 50th, and 60th days of age. Three milliliters of blood were collected into TEMPUS tubes (Applied Biosystems, Foster City, CA, USA) with RNA stabilizing factor. Samples were stored at −20 ◦C until further processing. Isolation of RNA was carried out using TEMPUS SPIN (Ambion, Waltham, MA, USA) according to the manufacturer's protocol (Supplementary File 1). One microgram of RNA was transcribed into cDNA using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's protocol.

A "No-RT" (non-reverse transcriptase) control was used for selected RNA samples to analyzed contamination in samples.

Gene expression analyses (Table S2 presents the reaction efficiency of each gene) were performed on an Illumina Eco system (Illumina, San Diego, CA, USA, Country)using TaqMan®MGB (Applied Biosystems, Foster City, CA, USA) probes (Table 1). Every sample was analyzed in triplicate in a final volume of 10 μL (Table S1). Amplification was performed according to the following protocol: polymerase activation at 95 ◦C (2 min) and 40 cycles at 95 ◦C for 15 s and 60 ◦C for 1 min. The *SDHA* and *HPRT* genes were used as housekeeping genes (Table 1).

**Table 1.** Probes used for amplification of Toll-like receptor (TLR) genes and housekeeping genes.


In addition, an analysis of the blood morphotic parameters was performed (Supplementary File 1).

#### *2.4. Statistical Analysis*

Data are presented as means ± standard error. The data were analyzed using SAS 9.4 software (SAS Institute Inc., Cary, NC, USA). The Shapiro–Wilk test was considered the best test to check the normality of the distribution of random variables. Because the data did not have a normal distribution, the Kruskal–Wallis test was used with immunostimulation and age as the effects. The degree of association between the parameters was examined using a non-parametric Spearman's rank correlation coefficient. Values ranging from 0.0 to 0.5, from 0.5 to 1.0, from −0.5 to 0.0, and from −1.0 to −0.5 indicate weak positive, strong positive, weak negative, and strong negative correlations, respectively.

#### **3. Results**

#### *3.1. Influence of Age on the Expression of TLR3, 4, and 7 mRNA*

The lowest expression of *TLR3* was observed during delivery (6.20 ± 0.89) (Figure 1—data presented from control group). After delivery, the level of *TLR3* mRNA increased. In the period between delivery and 60 days of age, the level of *TLR3* expression increased by 94.34%. We found a highly statistically significant difference between the age and expression of *TLR3* mRNA (*p* > 0.01), as presented in Table 2. From the 5th day of age, we found statistical differences between samples, with the highest expression level on the 60th day after delivery.

The highest level of *TLR4* mRNA expression was observed during the delivery (18.3 ± 2.6) of newborn foals. From the day of the delivery to 20th day of age, we observed a steady significant decrease of *TLR4* mRNA expression (Figure 1) in the blood of the examined foals. During the subsequent days of observation, the expression of the mRNA of this receptor remained at a similar level. The lowest expression value was observed at the 60th day of age (5.82 ± 0.96) (Table 2). Between the delivery day and the 60th day of age, the expression of mRNA for *TLR4* decreased by 76.89%. Statistical analysis showed highly statistically significant differences between the age of foals and the expression of *TLR4* mRNA (*p* < 0.01).

**Figure 1.** Trends in the change of TLR3 (**A**), TLR4 (**B**), and TLR7 (**C**) expression during foals' growth. Delivery, sample collected at delivery; 24 h, sample collected 24 h after delivery; 3 days, sample collected every 3rd day after delivery; 5 days, sample collected every 5 days after delivery; 10 days, sample collected every 10 days after delivery; 20 days, sample collected every 20 days after delivery; 30 days, sample collected 30 days after delivery; 40 days, sample collected every 40 days after delivery; 50 days, sample collected every 50 days after delivery; 60 days, sample collected 60 days after delivery. The means are reported with their standard errors.

The expression ofthe*TLR7* gene remained statistically unchanged throughout the experiment (Figure 1). The highest values were observed during the day of delivery (9.20 ± 1.20), while the lowest was observed 20 days after delivery (7.20 ± 0.61). There was no statistically significant correlation between the age and the expression of *TLR7* mRNA (*p* > 0.2366) (Table 2).


**Table 2.** Expression of *TLR3*, *TLR4*, and *TLR7* mRNA over the foals' subsequent days of age from the control group (mean ± standard error (SE).

<sup>1</sup> Delivery, sample collected at delivery; 24 h, sample collected 24 h after delivery; 3 days, sample collected 3 days after delivery; 5 days, sample collected 5 days after delivery; 10 days, sample collected 10 days after delivery; 20 days, sample collected 20 days after delivery; 30 days, sample collected 30 days after delivery; 40 days, sample collected 40 days after delivery; 50 days, sample collected 50 days after delivery; 60 days, sample collected 60 days after delivery. <sup>2</sup> Means are reported with their standard errors. Means with same letter in column show highly statistically significant differences (*p* < 0.01).

#### *3.2. Influence of Stimulation with Biotropina on the Expression of TLR3, 4, and 7 mRNA*

Analysis of changes in the expression of *TLR3* mRNA after the injection of the immunostimulant (Group E) showed a decrease by 41.65% (Table 3), while in the control group (Group C), a dynamic increase in the expression was observed until the last day of observation (116.22 ± 13.93). Highly statistically significant differences were found between immunostimulation and the expression of *TLR3* mRNA.


**Table 3.** Influence of stimulation with Biotropina on the expression of mRNA for selected Toll-like receptors (*TLR3*, *TRL4*, *TLR7*) (mean ± SE).

<sup>1</sup> Delivery, sample collected at delivery; 24 h, sample collected 24 h after delivery; 3 days, sample collected 3 days after delivery; 5 days, sample collected 5 days after delivery; 10 days, sample collected 10 days after delivery; 20 days, sample collected 20 days after delivery; 30 days, sample collected 30 days after delivery; 40 days, sample collected 40 days after delivery; 50 days, sample collected 50 days after delivery; 60 days, sample collected 60 days after delivery. <sup>2</sup> Means are reported with their standard errors; Group C, control group; Group E, experimental Biotropina-stimulated group (injection at the 35th and 40th days after delivery). \* Means in row/line for receptor show significant differences (*p* < 0.05); \*\* Means in row/line for receptor show highly statistically significant differences (*p* < 0.01).

The level of *TLR4* mRNA expression at 30 days of age was similar in both groups (Table 3). After Biotropina injection, *TLR4* mRNA expression increased by 41.33% (Group E), while the expression of *TLR4* mRNA in foals from Group C decreased by 20.22%. On the following days after immunostimulation, *TLR4* mRNA expression in the foals from Group E was higher, but we did

not find statistically significant differences between groups. Statistical analysis showed statistically significant differences in *TLR4* mRNA expression after immunostimulation.

The initial expression of *TLR7* mRNA before first suckling was higher in Group C. In Group E, the highest level of expression was observed during delivery at 7.57 ± 0.88 (Table 3). In Group C, the expression was higher during the experiment compared to Group E, and we found highly statistically significant differences between groups (*p* < 0.001). No statistically significant differences were found between the expression of *TLR7* mRNA and immunostimulation.

The Spearman's rank correlation test showed a strong negative correlation between the*TLR3* and *TLR4* genes, and lack of correlation between *TLR3* and *TLR7*as well as*TLR4* and *TLR7* (Table 4).


**Table 4.** Correlations (*p*-value) of the expression of *TLR3*, *TLR4,* and *TLR7* mRNA over the subsequent days of age of the foals from the control group.

<1 h, sample collected at delivery; 24 h, sample collected 24 h after delivery; 3 days, sample collected 3 days after delivery; 5 days, sample collected 5 days after delivery; 10 days, sample collected 10 days after delivery; 20 days, sample collected 20 days after delivery; 30 days, sample collected 30 days after delivery; 40 days, sample collected 40 days after delivery; 50 days, sample collected 50 days after delivery; 60 days, sample collected 60 days after delivery. Correlations (*p*-value) bolded show significant differences (*p* < 0.05) while underlined show highly significant differences (*p* < 0.01).

#### *3.3. Influence of Age and Stimulation with Biotropina on the Level of Blood Morphotic Elements*

Statistical analysis showed (Table 5):




<1, sample collected at delivery; 24 h, sample collected 24 h after delivery; 3 days, sample collected 3 days after delivery; 5 days, sample collected 5 days after delivery; 10 days, sample collected 10 days after delivery; 20 days, sample collected 20 days after delivery; 30 days, sample collected 30 days after delivery; 40 days, sample collected 40 days after delivery; 50 days, sample collected 50 days after delivery; 60 days, sample collected 60 days after delivery. 2 Means are reported with their standard errors. Group C, control group; Group E, experimental Biotropina-stimulated group (injection on the 35th and 40th days after delivery). \* Means in row/line for receptor show significant differences (*<sup>p</sup>* < 0.05); \*\* means in row/line for receptor show highly statistically significant differences (*<sup>p</sup>* < 0.01).

1

#### **4. Discussion**

Infectious diseases are common in foals between the first and fifth months of age. Analysis of the concentrations of immune system components during this period in healthy and infected foals may help understand the basics of the maturation of the immune system and also better understand infection mechanisms [22]. To the best of our knowledge, there are very limited reports where weekly collections and the expression of immune-related genes have been performed. Therefore, our results may be interesting for better understanding the changes during the first weeks of a foal's life and the changes it undergoes during this time. Most studies report data from the first 24 h of a foal's life, from the first 42 days, and very often from adult horses. Moreover, most data include thoroughbreds, while we performed our analysis on a primitive horse breed that is known for their adaptation to harsh conditions. While this study was performed on a primitive domestic breed, the results from this study may differ from potential results if performed on selectively breed domestic horses. Flaminio et al. [22] reported that healthy lymphocytes of healthy foals were the lowest at birth and that values increased until the sixth month of age. In our study, we obtained similar results; however, in Polish Pony, lymphocyte counts were higher (Table 5).

#### *4.1. Changes in TLR4 Gene Expression*

Because of the increased vulnerability of foals to some pathogens (e.g., *Rhodococcus equi*), it seems reasonable to analyze changes in the expression of the genes that are responsible for the recognition of the conserved constituents of pathogens [23]. In the present study, a highly significant influence of age and stimulation with Biotropina was observed for *TLR4*. Data available in the literature indicate that the influence of age on *TLR4* expression in horses is contradictory. Vendrig et al. [24] found no differences between the expression of *TLR4* in the blood mononuclear cells of foals at 12 h of life or in adult horses. Stimulation with lipopolysaccharides (LPS) resulted in higher *TLR4* mRNA expression in adult horses, while no response to LPS stimulation was found in foals in an in vivo study [24]. In contrast, Tessier et al. [25], having compared *TLR4* mRNA expression in umbilical cord blood and peripheral blood from adult horses, demonstrated higher*TLR4* mRNA expression in umbilical cord blood in response to LPS administration. The influence of age on the expression of *TLR4*was observed by Hansen et al. [26] in horses aged between 5 and 27 years. Higher *TLR4* mRNA levels were found in younger horses, but the decrease in mRNA levels with age were not statistically significant. *TLR4* mRNA levels were higher in blood mononuclear cells compared to the same cells from pulmonary vascular secretions. Osorio et al. [27] and Strong et al. [28] evaluated the expression of the *TLR4* gene in the first weeks of a calf's life, and their results were similar to the one produced in our study. The highest level of*TLR4* gene expression was observed after birth, and a statistically significant decrease in expression was reported during the following days. A similar trend was identified in our study. Yerkovich et al. [29] and Levy et al. [30] showed that the expression of the *TLR4* gene was significantly higher in peripheral blood in premature and full-term infants than in adults, both before and after LPS stimulation [31]. This trend, indicating a decreasing expression of *TLR4* with age, was also found in humans and mice [32,33]. On the other hand, some results illustrate a higher expression of *TLR4* in newborns compared to adults [34] or decreasing expression of *TLR4* after stimulation [35]. The differences in *TLR*4 expression revealed in these studies and in our experiment may be caused by different concentrations of LPS in the stimulation.

#### *4.2. Changes in TLR3 Gene Expression*

We found that *TLR3* mRNA levels increased with the age of the foals. Our results are in agreement with other reports about horses and mostly human newborns [25,36–38]. Interestingly, there are reports proving epigenetics control *TLR3* expression mechanisms [36]. As Porras et al. [36] reported in their results from healthy donors, it can be presumed that a low level of *TLR3* in newborns is a developmentally desirable trend. In a mouse model, Zhang et al. [37] reported higher abortion rates linked with higher *TLR3* levels. It was also mentioned that *TLR3* expression was age-dependent, which can be confirmed by our results. As *TLR3* binds double-stranded RNAs (dsRNAs), its decreased level may increase susceptibility to viral infection in young foals. In our study, the lowest level was recorded before the first suckling, which is in agreement with other reports in premature infants [38] and newborns [39,40]. The data presented above, and the results obtained in our study of the *TLR3* gene, may explain the higher incidence of equine herpesvirus-1 (EHV-1) and equine herpesvirus-4 (EHV-4), responsible for massive respiratory tract infections in foals and young horses. A severe course and high mortality due to contracting equine viral arteritis (EVA) in young horses may also be a result of the decreased expression of *TLR3*. Hussey et al. [41] suggested that *TLR3* plays an essential role in recognizing EHV-1 infections. In our study, a decrease in *TLR3* gene expression was observed after stimulation with Biotropina. We analyzed the level of expression in foals up to 60 days of age, but the literature indicates that the immune system of horses develops most intensively until about 90 days of life [42]. Foals have all of the components of an immune system characteristic of adult horses—but many mechanisms of the immune response have yet to mature. The results indicate that activation of horse monocytes by ligands for the *TLR2* and *TLR4* genes increases their expression, but not that of *TLR3*. Additionally, *TLR3* gene expression decreases with the increase of *TLR4* gene expression after the stimulation of monocytes [43].

#### *4.3. Changes in TLR7 Gene Expression*

*TLR7* is responsible for recognizing guanidine-rich, single-stranded viral RNA (ssRNA) and is an important mediator of the peripheral immune response. Asquith at al. [41] and Slavica at al. [39] observed no effect of age on*TLR7* gene expression in newborns, similar to Talmadge et al. [22] in horses. Belnoue et al. [44] reported significantly higher levels of *TLR7* gene expression in two-week-old foals compared to adult horses. Harrington et al. [9] found neither an age-dependent pattern in the expression of the *TLR7*/8 genes, nor did they detect the effect of imidazoquinol R848 stimulation on its expression, despite increasing the levels of IL-6 and IL-8. Our results also did not confirm any relationship between a foal's age and expression of *TLR7*. Until now, little was known about the signaling mechanisms of Toll-like receptors in foals. Identifying the receptors and describing ligands that react with them can provide new insights into immunological responses and can also point to new pathways in the field of therapy and prevention of diseases, particularly infectious ones.

#### **5. Conclusions**

In summary, on the basis of the results obtained, it was concluded that the expression of the *TLR3* and *TLR4* genes in peripheral blood cells is dependent on age. The expression of the *TLR3* and *TLR4* genes, as well as the levels of immunoglobulins and interleukins, can be modulated by stimulation with the pharmacological agent Biotropina. This experiment demonstrated a strong negative correlation between *TLR3* and *TLR4* gene expression. Detailed knowledge of the molecular mechanisms of immunoglobulin synthesis appears necessary for a better understanding of foal immunity maturity and its influencing factors. At the same time, this experiment encourages studies regarding the influence of the signaling cascade's proteins on the primary immunological response, providing an opportunity to develop extremely precise methods of regulating acquired immunity. There is still little information about the maturity of a horse's immune system in the pre- and postnatal period, which negatively affects the planning of health protection strategies for foals.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/11/1966/s1, Supplementary File 1: TEMPUS SPIN manufacturer protocol and Morphotic Blood Parameters, Figure S1: Birth system alarm, Table S1: Mastermix reaction, Table S2: Reaction efficiency.

**Author Contributions:** Conceptualization, A.M., A.O. and A.C.-S.; methodology A.M. and M.O.; collection and analysis of samples, A.M., Ł.M. and A.O.; Data interpretation A.M., Ł.M., M.O. and A.C.-S.; writing A.M. and Ł.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financed by the Ministry of Science and Higher Education of the Republic of Poland (funds for statutory activity, DS 3208 and SUB.2015-D201).

**Conflicts of Interest:** The authors declare that they have no competing interests. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **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* **Genome-Wide DNA Methylation Changes of Perirenal Adipose Tissue in Rabbits Fed a High-Fat Diet**

#### **Jiahao Shao 1,**†**, Xue Bai 1,**†**, Ting Pan 2, Yanhong Li 1, Xianbo Jia 1, Jie Wang <sup>1</sup> and Songjia Lai 1,\***


Received: 3 November 2020; Accepted: 24 November 2020; Published: 26 November 2020

**Simple Summary:** Obesity is spreading rapidly in most countries and regions, becoming a considerable public health concern because it is associated with type II diabetes mellitus, fatty liver disease, hypertension, and even certain cancers. The biological effects of caloric restriction are closely related to epigenetic mechanisms, including DNA methylation. Here, rabbits were used as a model to study the effect of a high-fat diet on the DNA methylation profile of perirenal adipose tissue. The results indicate that 2906 genes associated with differentially methylated regions were obtained and were involved in the PI3K-AKT signaling pathway (KO04151), linoleic acid metabolism (KO00591), DNA replication (KO03030), and MAPK signaling pathway (KO04010). In conclusion, high-fat diet may cause changes in the DNA methylation profile of adipose tissue and lead to obesity.

**Abstract:** DNA methylation is an epigenetic mechanism that plays an important role in gene regulation without an altered DNA sequence. Previous studies have demonstrated that diet affects obesity by partially mediating DNA methylation. Our study investigated the genome-wide DNA methylation of perirenal adipose tissue in rabbits to identify the epigenetic changes of high-fat diet-mediated obesity. Two libraries were constructed pooling DNA of rabbits fed a standard normal diet (SND) and DNA of rabbits fed a high-fat diet (HFD). Differentially methylated regions (DMRs) were identified using the option of the sliding window method, and online software DAVID Bioinformatics Resources 6.7 was used to perform Gene Ontology (GO) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of DMRs-associated genes. A total of 12,230 DMRs were obtained, of which 2305 (1207 up-regulated, 1098 down-regulated) and 601 (368 up-regulated, 233 down-regulated) of identified DMRs were observed in the gene body and promoter regions, respectively. GO analysis revealed that the DMRs-associated genes were involved in developmental process (GO:0032502), cell differentiation (GO:0030154), and lipid binding (GO:0008289), and KEGG pathway enrichment analysis revealed the DMRs-associated genes were enriched in linoleic acid metabolism (KO00591), DNA replication (KO03030), and MAPK signaling pathway (KO04010). Our study further elucidates the possible functions of DMRs-associated genes in rabbit adipogenesis, contributing to the understanding of HFD-mediated obesity.

**Keywords:** DNA methylation; high-fat diet; rabbits

#### **1. Introduction**

From the last 5 decades, the incidence of obesity has sharply increased, becoming one of the most considerable threats to human health because it is associated with the risk of type II diabetes mellitus, fatty liver disease, hypertension, and even certain cancers [1]. Obesity is a multifactorial pathological process, and genetic, environmental, and behavioral factors influence the development of obesity [2]. Nowadays, an imbalance between energy intake and expenditure is a major contributor to fat deposition in individuals predisposed to obesity [3]. Fat deposition is characterized by an increase in the number and size of adipocytes, and its process is closely related to physiological homeostasis, far beyond simple fat storage [4]. HFD has been shown to induce obesity in animal models and humans, and further induce a variety of obesity-related clinical diseases, such as osteoporosis, inflammation, and even neurodegeneration [5–7]. Perirenal fat, as part of abdominal visceral fat, is often used to elucidate the molecular and pathophysiological mechanisms of metabolic disorders associated with obesity or adipose development, because it is closely related to kidney injury, metabolism of triacylglycerol, and other metabolic regulation [8]. For example, detailed studies have shown that the perirenal fat thickness in obese patients could be a valuable marker to define the risk of developing hypertension and kidney dysfunction [9,10]. The expression profile of perirenal fat microRNA was changed during different growth stages of rabbits, and the differential microRNA expression was enriched for the MAPK signaling pathway, Wnt signaling pathway, aldosterone synthesis, and secretion pathways [11].

First proposed by Waddington in 1942, epigenetics refers to heritable changes in gene expression without an altered DNA sequence [12]. Epigenetics is caused by the interaction of environmental factors and intracellular genetic material, such as dietary factors, microRNA, and genomic imprinting, etc. Noteworthily, the biological effects of caloric intake are closely related to epigenetic mechanisms, including chromatin remodeling and DNA methylation [13]. DNA methylation of leptin and adiponectin promoters in obese children is associated with BMI, dyslipidemia, and insulin resistance [14]. These observations support the hypothesis that epigenetic modifications might underpin the development of obesity and related metabolic disorders. Hypermethylation of the pro-opiomelanocortin and serotonin transporter genes has been positively associated with childhood or adult obesity [15]. HFD changes the methylation status of *Casp1* and *Ndufb9* genes in obese mice, which are related to liver lipid metabolism and liver steatosis [16]. In addition, the leptin promoter was hypermethylated and *Ppar-*α promoter was hypomethylated in oocytes of mice fed with HFD, and the same changes were also observed in the liver of female offspring [17]. However, few studies have reported the changes in perirenal adipose tissue methylation profile in HFD-induced obese rabbits.

To further understand the epigenetic mechanisms influencing fat metabolism in obese rabbits, we investigated the role of DNA methylation in perirenal adipose tissue by sequencing and analyzing DNA methylation libraries from rabbits fed a standard normal diet (SND) and a high-fat diet (HFD).

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

#### *2.1. Animals*

A total of 24 female Tianfu black rabbits from a strain breed at the Sichuan Agricultural University in China were randomly divided into two groups and fed either a standard normal diet (SND) or a high-fat diet (HFD; 10% lard was added to the standard normal diet) for four weeks. The composition and nutrient content of the standard normal diet (SND) and the high-fat diet (HFD) were described in our previous report [18]. At the beginning of the trial, all rabbits were 35 days of age and housed individually in a clean iron cage (600 × 600 × 500 mm) and kept in an environmentally controlled room. Rabbits were free to access water and fed twice a day. At the end of the trial, rabbits were screened for obesity using the body mass index (BMI; BMI = bodyweight (kg)/height2 (m)), and three rabbits from each group meeting the experimental requirements were selected for sampling. All experimental protocols were performed under the direction of the Institutional Animal Care and

Use Committee from the College of Animal Science and Technology, Sichuan Agricultural University, China (DKY-B2019202015, 5 December 2019).

#### *2.2. DNA Extraction*

Perirenal adipose tissue samples were collected immediately after rabbits were euthanized (shock and bleed treatment). Tissue blocks were placed in 4 mL EP tubes and stored in a −80 ◦C freezer. Total DNA from perirenal adipose tissue was extracted using a commercial TIANamp Genomic DNA extraction kit (Tiangen, Beijing, China), following the manufacturer's instructions. Subsequently, the purity and concentration of DNA were assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, Carlsbad, CA, USA), and only DNA meeting quality criteria (thresholds: A260/A280 ≈ 1.8; concentration ≥ 200 ng/μL) was used for the trial.

#### *2.3. DNA Methylation Library Construction and Sequencing*

To identify genome-wide DNA methylation changes in perirenal adipose tissue induced by HFD, two libraries were constructed by pooling the DNA samples from three SND rabbits and three HFD rabbits. Briefly, DNA was fragmented by sonication to 100 to 500 bp fragments. The fragments were end-repaired using T4 DNA polymerase and Klenow enzyme and adaptors were ligated after generating 3'dA overhangs. Bisulfite treatment was conducted using the ZYMO EZ DNA Methylation-Gold kit (Zymo Research, Orange, CA, USA), following the manufacturer's protocol. After desalting, fragments of sizes ranging from 220 to 320 bp were isolated using a 15% PAGE gel and amplified by adaptor-mediated PCR. Lastly, the libraries were sequenced using the Illumina HiSeq 4000 platform (Illumina, San Diego, CA, USA) by Chengdu Life Baseline Technology Corporation, China.

#### *2.4. Processing and Comparison of Sequencing Data*

By removing adapter sequences and low-quality reads containing more than 50% low-quality bases (quality score < 5), clean reads were retained. Clean reads were aligned to the rabbit reference genome (GCF\_000003625.3) with software BSMAP 2.90 (http://code.google.com/p/bsmap). Two forward strands, i.e., BSW (++) and BSC (−+) were used as references. The accuracy of DNA methylation detection depends on the conversion efficiency of cytosine, and the incomplete transformation of cytosine in sequences may lead to false-positive results. Here, lambda phage DNA was used as a control group to calculate the bisulfite conversion rate.

#### *2.5. Methylation Site Detection*

The methylation C sites were determined using the method described in a previous study [19]. Briefly, a binomial distribution test was performed for methylated reads number and non-methylated reads number at C sites. C sites were identified as the methylation C sites when the number of reads was greater than or equal to the binomial distribution expected value and the total effective coverage was greater than or equal to four.

#### *2.6. Methylation Level Analysis*

The average genome-wide methylation level reflects the overall characteristics of the methylation pattern of the genome. DNA methylation occurs in three sequence contexts: CG, CHG, and CHH (H = A, C, or T). The average methylation levels of CG, CHG, and CHH were calculated based on the percentage of methylated cytosine in the entire genome, chromosome, and genomic functional elements. For each type of sequence (CG, CHG, and CHH), the average methylation level was calculated according to the following formula: the average methylation level = methylated reads / (methylated reads + non-methylated reads) × 100%. To assess the association between sequence characteristics and methylation bias, we calculated the methylation percentage of nine bases upstream and downstream of the methylation site.

#### *2.7. Searching for Di*ff*erentially Methylated Regions (DMRs)*

DMRs were identified using the option of a sliding window method. Briefly, the sliding windows, which were used for further analysis, had to meet the following criteria: (a) the depth in each cytosine should be more than four in each sample, and each C site should cover at least four methylation reads; (b) the number of selected cytosine should be larger than five; (c) after calculating mean methylation level of each sample, the fold change of mean methylation level between the two samples should be larger than two. After repeating extension steps, the merged regions with *p* < 0.05 were defined as DMRs.

#### *2.8. Functional Enrichment Analysis of Di*ff*erentially Methylated Genes*

To explore the role of epigenetic variation in biological processes and pathways, online software DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/home.jsp) was used to perform Gene Ontology (GO) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of DMRs-associated genes. GO analysis can be used to identify the performance of the gene product and contains three types of information: cellular component, molecular function, and biological processes. KEGG is the main public database that integrates the genome, chemistry, and system function information, particularly the set of genes associated with the systemic functions of cells, organisms, and ecosystems. Differences were considered to be statistically significant at *p* < 0.05.

#### **3. Results**

#### *3.1. Quality Assessment of Sequencing Data*

After raw reads were processed, a total of 1,221,455,488 clean reads were obtained from methylation sequencing libraries (Table S1). The clean reads were mapped to the rabbit reference genome, and the mapping rate was 84.910% in the SND group and 84.730% in the HFD group, respectively. The bisulfite conversion rate was 99.550% for SND, and 99.520% for HFD. In addition, the effective coverage rate of C base in each chromosome ranged from 89.892% to 97.577% and ranged from 91.822% to 97.804% in different functional genomic elements (Tables S2 and S3).

#### *3.2. Methylation Level Analysis*

Genome-wide methylation level analysis showed that the methylation level of C, CHG, and CHH in the HFD group was higher than in the SND group but the CG methylation level in the HFD group was lower than in the SND group (Table S4). Results of the methylation level C, CG, CHG, and CHH on different chromosomes are shown in Table S5. The greatest differences in C, CG, CHG, and CHH between the two groups were found on chromosome 20, chromosome X, chromosome 1, and chromosome 11, reaching 0.569%, 2.736%, 0.056%, and 0.047%, respectively. In addition, we classified the various functional genomic elements into promoter, CDS, intron, mRNA, downstream, CpGIsland, ncRNA, and transposons. Compared with the SND group, the methylation level of C, CHG, and CHH in each functional genomic element was increased in the HFD group (Table S6). However, based on comparison with the SND group, promoter, intron, mRNA, downstream, and ncRNA methylation levels were decreased in CG, and only CDS, CpGIsland, and transposons methylation levels were increased in CG.

#### *3.3. Genome-Wide Characteristics of Methylated C Bases*

The percentage of methylated C bases in CG were highest, reaching 94.795% (SND) and 94.843% (HFD) but rarely cytosine methylation was found in CHH and CHG. In addition, we calculated the methylation percentage of nine bases (methylated C at the fourth base) upstream and downstream of the methylated site. As shown in Figure 1, CG, CAG, and CAC were the most likely sites to be methylated in both SND and HFD groups.

**Figure 1.** Genome-wide characteristics of methylated C bases. Sequence characteristics of bases near mCG, mCHG, and mCHH in the SND and HFD group.

#### *3.4. Analysis of Di*ff*erentially Methylated Regions (DMRs)*

A total of 12,230 DMRs were identified in the genome of the HFD group compared to the SND group. Chromosome 21 was the chromosome with the least amount of DMRs and chromosome 13 was the chromosome with the most amount of DMRs (Figure 2a,b). The total length of DMRs in each chromosome is shown in Table S7. In addition, the DMRs were mapped to the gene body and promoter regions, and 2305 (1207 up-regulated, 1098 down-regulated) and 601 (368 up-regulated, 233 down-regulated) methylated genes were obtained, respectively. Some genes involved in adipocyte growth and development have also been identified, including *ACE2*, *AGTR1*, *IGF1R*, and *ACSL4*.

#### *3.5. GO and KEGG Enrichment Analysis*

To better study the biological functions of the DMRs-associated genes, we used online software DAVID Bioinformatics Resources 6.7 (http://david.abcc.ncifcrf.gov/home.jsp) to carry out gene ontology (GO) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis. GO analysis of the overlapping DMRs-associated genes in the gene body regions found a total of 6310 enriched GO terms (4796 biological processes (BP), 579 cellular components (CC), and 935 molecular functions (MF)), of which 12.570% were significantly enriched (*p* < 0.05) (Table S8). The main GO terms involved in overlapping DMRs-associated genes in the gene body regions included the developmental process (GO:0032502), cell differentiation (GO:0030154), and lipid binding (GO:0008289). The top 10 significantly enriched terms in the BP, CC, and MF categories are shown in Figure 3a. The KEGG pathway analysis showed that overlapping DMRs-associated genes in the gene body regions were enriched in 314 pathways including the PI3K-AKT signaling pathway (KO04151), linoleic acid metabolism (KO00591), and pathways for DNA replication (KO03030). Thirty-nine of these pathways (12.420%) were significantly enriched (*p* < 0.05, Table S9). In addition, a scatter analysis was carried out for the 20 most significant pathways to intuitively show the significance of these pathways (Figure 3b).

**Figure 2.** Analysis of differentially methylated regions (DMRs). Two chromosomes with the least (**a**) and the most (**b**) amount of DMRs. The outer ring represents the position of the genomic chromosome; the second circle is the DMRs region: the red area represents the higher methylation level of HFD compared to the SND group and the green area represents the lower methylation level of HFD compared to the SND group; the third circle represents the methylation rate of each site of sample HFD; the fourth circle represents the methylation rate of each site of sample SND; the fifth circle represents the difference of methylation rate.

**Figure 3.** GO and KEGG enrichment analysis. (**a**) GO analysis of the overlapping DMRs-associated genes in the gene body regions and promoter regions. (**b**) KEGG pathway analysis of the overlapping DMRs-associated genes in the gene body regions and promoter regions. Rich factor = (DMRs-associated genes annotation in term/genes annotation in term)/(DMRs-associated genes with KEGG annotation / all genes with KEGG annotation).

GO analysis of the overlapping DMRs-associated genes in the promoter regions showed enrichment of 2223 biological processes (BP), 311 cellular components (CC), and 405 molecular functions (MF), of which 173 BP (7.780%), 20 CC (6.430%), and 37 MF (9.140%) were significantly enriched (Table S10). The significantly enriched GO terms mainly include positive regulation of lipid biosynthetic process (GO:0046889), regulation of cholesterol metabolic process (GO:0090181), and regulation of lipid biosynthetic process (GO:0046890). The top 10 significantly enriched GO terms in the BP, CC, and MF categories of GO analysis are shown in Figure 3a. KEGG pathway analysis found 266 enriched pathways including the MAPK signaling pathway (KO04010) (Table S11). The top 20 significantly enriched pathways are presented in Figure 3b.

#### **4. Discussion**

DNA methylation represents an important epigenetic marker because it is associated with chromosomal structural changes, embryonic development, expression of imprinted genes, and causing corresponding diseases, including X chromosome inactivation and DNA unwinding [20–22]. Nowadays, obesity prevention and treatment strategies have been unsuccessful, and DNA methylation is one of the epigenetic modifications associated with obesity [23]. The rabbit is an ideal material to study obesity due to its lipid metabolism and obesity-related clinical manifestations similar to those of humans [24,25]. Thus far, some DNA methylation related studies have been investigated in rabbit models but studies of the changes in perirenal adipose tissue methylation profile in HFD-induced obese rabbits have not been carried out. In this study, DNA methylation patterns were investigated in rabbit models to understand obesity-related DNA methylation changes.

Here, the mapping rates were 84.910% and 84.730% in the SND group and HDF group, respectively. The bisulfite conversion rates were 99.550% (SND) and 99.520% (HFD) in the two groups, which was consistent with previous research, indicating that the libraries were high quality and reliable [26,27]. Methylation is a dynamic process in cells, which can be regulated by methylation and demethylation. The average methylation level of the whole genome reflects the overall characteristics of the genome methylation profile. The results of the genome-wide methylation level analysis in this study were similar to those in mice [28]. The methylation level of CG was higher than the methylation level in C, CHG, and CHH. However, this result is different from that of plant Arabidopsis Thaliana. The plant genome has extensive methylation at the CHG site [29]. CG methylation was maintained by Dnmt1. CHH methylation and some CHG methylation is usually maintained by the activity of the conserved Dnmt3. The high level of CHG methylation seen in Arabidopsis thaliana is maintained by plant-specific methyltransferase [30]. In addition, research on chickens suggests that promoter DNA methylation generally affects chromatin structure and is a signal to inhibit gene transcription, and promoter regions are lowly methylated [31]. Our study also found that promoter regions showed a lower methylation level than other regions. However, a study in mice fed with HFD showed that promoter regions are hypermethylated [32]. Therefore, we hypothesized that differences in methylation level may be species-specific.

The results of genome-wide characteristics of methylated C bases showed that the proportion of mCG was the highest, while the cytosine methylation was low in CHH and CHG. Some studies have shown that no enzyme can maintain mCHG during DNA replication in animals, so the sites of CHG type in animal cells generally show a very low level of methylation. CHH can only rely on the methylation mechanism, so CHH methylation is easily lost in the process of DNA replication and is generally in the state of hypomethylation [33]. The results of this study showed that the characteristics of methylation in the rabbit genome were similar to those in other animals.

DMRs refer to the regions of DNA molecules with different methylation status in two samples. The identification of DMRs is the first step towards the study of DMRs-associated genes [34]. In our study, a total of 2906 DMRs were identified, and 2305 (1207 up-regulated, 1098 down-regulated) and 601 (368 up-regulated, 233 down-regulated) methylated genes were associated with differentially methylated regions. Many genes are related to adipocyte growth and development. For example, as the members of the renin-angiotensin system (RAS), *ACE2* and *AGTR1* were reported to participate in the development and progression of obesity [35,36]. *PPAR*γ and *aP2* are important transcription factors in the development and function of the adipose tissue and marker of lipogenesis [37]. Previous studies

showed that inhibition of *IGF1R* decreased the expression of PPARγ, thereby inhibiting lipogenesis [38]. Moreover, *ACSL4* plays a role in the regulation of lipid metabolism. *ACSL4* was expressed throughout the entire differentiation process in pig preadipocytes and showed a similar expression trend with lipogenesis-associated genes *PPAR*γ and *aP2* [39].

Gene ontology (GO) analysis is a reliable bioinformatics tool for understanding the characteristics of genes and gene products. The significantly enriched terms in the BP, CC, and MF categories indicated the possible roles of the DMRs-associated genes in regulating obesity. The significantly enriched GO terms showed correlation with adipocyte lipid metabolism and metabolisms, such as lipid binding (GO:0008289), positive regulation of lipid biosynthetic process (GO:0046889), regulation of cholesterol metabolic process (GO:0090181), developmental process (GO:0032502), and cell differentiation (GO:0030154). Some terms were related to adipocyte development, including cytoskeletal protein binding (GO:0008092), tubulin binding (GO:0015631), calcium ion binding (GO:0005509). Cytoskeletal remodeling and cell–cell interaction are a necessary step in the transformation of preadipocytes into mature adipocytes, and adipocyte development is dependent on α-tubulin [40,41]. Calcium is a complex mediator in adipogenesis because it regulates numerous cellular processes [42]. Furthermore, other GO items related to hormones and enzymes were also significantly enriched, such as regulation of glucocorticoid secretion (GO:2000849), N-acetyltransferase activity (GO:0008080), and phosphoric diester hydrolase activity (GO:0008081). Increasing evidence suggests that excess glucocorticoids leads to increased fat mass and obesity through the accumulation of adipocytes [43]. Acetyltransferase is a regulator of adipogenesis and lipid metabolism, and its regulatory mechanism is mainly transcription and post-translation modifications [44]. Phosphoric diester hydrolase is a regulator of systemic glucose and insulin homeostasis [45]. Interference of phosphoric diester hydrolase expression in 3T3-L1 adipocytes caused a dramatic decrease in adipocyte differentiation key gene (*PPAR*γ, *aP2*) and lipid accumulation [46].

Adipogenesis is a complex process involving an elaborate network of transcription factors and signaling pathways. Results of KEGG analysis showed that DMRs-associated genes were mainly involved in the PI3K-AKT signaling pathway (ko04151), linoleic acid metabolism (KO00591), DNA replication (KO03030), and MAPK signaling pathway (KO04010). The PI3K-AKT signaling pathway is a key regulator in cell proliferation, differentiation, and apoptosis [47]. Activation of the PI3K-AKT signaling pathway promotes the expression of marker genes involved in adipogenesis and glucose uptake [48]. In our study, 70 DMRs-associated genes were enriched in the PI3K-AKT signaling pathway, thereby revealing that these DMRs-associated genes may be essential for adipogenesis. Linoleic acid metabolism (KO00591) is also associated with adipogenesis. Linoleic acid can be converted to the metabolically active arachidonic acid, which has roles in inducing inflammation and adipogenesis. Excessive intake of linoleic acid results in increasing magnitudes of adiposity, inflammatory cytokines, and insulin resistance [49]. In addition, it is becoming clear that DNA replication (KO03030) and the MAPK signaling pathway (KO04010) play an important role in adipocyte growth and development [50,51]. Thus, the results of our study indicate that these DMRs-associated genes might be an important regulator in adipogenesis. However, due to the limitation of experimental conditions, such as pooled samples, only one library per group, sequencing methods, etc., functional verification of these DMRs-associated genes will be important to consider in the future.

#### **5. Conclusions**

In conclusion, our study indicates that a high-fat diet may affect genes associated with adipogenesis by altering DNA methylation patterns. We identified 2906 methylated genes, of which, *ACE2*, *AGTR1*, *IGF1R*, and *ACSL4* may have a key role in adipogenesis. These genes may be involved in the regulation of adipogenesis through the PI3K-AKT signaling pathway (KO04151), linoleic acid metabolism (KO00591), DNA replication (KO03030), and MAPK signaling pathway (KO04010).

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/12/2213/s1, Table S1: Data output and comparative statistics, Table S2: The effective coverage rate of C base in each chromosome, Table S3: The effective coverage rate of C base in different functional genomic elements, Table S4: Results of genome-wide methylation level analysis, Table S5: Results of the methylation level C, CG, CHG, and CHH on different chromosomes, Table S6: Methylation level in various functional genomic elements between standard normal diet (SND) and high-fat diet (HFD) group, Table S7: DMRs number and length statistics of standard normal diet (SND) and high-fat diet (HFD) group, Table S8: GO analysis of the overlapping DMRs-associated genes in the gene body regions, Table S9: KEGG pathway analysis of the overlapping DMRs-associated genes in the gene body regions, Table S10: GO analysis of the overlapping DMRs-associated genes in the promoter regions, Table S11: KEGG pathway analysis of the overlapping DMRs-associated genes in the promoter regions.

**Author Contributions:** X.J., J.W., and S.L. conceived and designed the study; T.P. and Y.L. collected data and conducted the research; J.S. and X.B. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the China Agricultural Research System (Grant No. CARS-44-A-2).

**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* **Embryonic Thermal Manipulation A**ff**ects the Antioxidant Response to Post-Hatch Thermal Exposure in Broiler Chickens**

#### **Khaled M. M. Saleh 1, Amneh H. Tarkhan <sup>1</sup> and Mohammad Borhan Al-Zghoul 2,\***


Received: 23 November 2019; Accepted: 10 January 2020; Published: 13 January 2020

**Simple Summary:** The broiler chicken is one of the most important livestock species in the world, as it occupies a major role in the modern human diet. Due to uneven artificial selection pressures, the broiler has increased in size over the past few decades at the expense of its ability to withstand oxidative damage, the latter of which is often a byproduct of thermal stress. In order to attenuate the effects of heat stress, thermal manipulation (TM), which involves changes in incubation temperature at certain points of embryonic development, is increasingly being presented as a way in which to improve broiler thermotolerance. Therefore, the objective of this study was to investigate how TM might affect broiler response to post-hatch thermal stress in the context of the genes that help combat oxidative damage, namely the catalase, NADPH oxidase 4 (*NOX4*), and superoxide dismutase 2 (*SOD2*) genes. Expression of all three aforementioned genes differed significantly between TM and control chickens after exposure to cold and heat stress. Conclusively, TM may act as a viable mode of preventative treatment for broilers at risk of thermally induced oxidative stress.

**Abstract:** Thermal stress is a major source of oxidative damage in the broiler chicken (*Gallus gallus domesticus*) due to the latter's impaired metabolic function. While heat stress has been extensively studied in broilers, the effects of cold stress on broiler physiologic and oxidative function are still relatively unknown. The present study aimed to understand how thermal manipulation (TM) might affect a broiler's oxidative response to post-hatch thermal stress in terms of the mRNA expression of the catalase, NADPH oxidase 4 (*NOX4*), and superoxide dismutase 2 (*SOD2*) genes. During embryonic days 10 to 18, TM was carried out by raising the temperature to 39 ◦C at 65% relative humidity for 18 h/day. To induce heat stress, room temperature was raised from 21 to 35 ◦C during post-hatch days (PD) 28 to 35, while cold stress was induced during PD 32 to 37 by lowering the room temperature from 21 to 16 ◦C. At the end of the thermal stress periods, a number of chickens were euthanized to extract hepatic and splenic tissue from the heat-stressed group and cardiac, hepatic, muscular, and splenic tissue from the cold-stressed group. Catalase, *NOX4*, and *SOD2* expression in the heart, liver, and spleen were decreased in TM chickens compared to controls after both cold and heat stress. In contrast, the expression levels of these genes in the breast muscles of the TM group were increased or not affected. Moreover, TM chicks possessed an increased body weight (BW) and decreased cloacal temperature (TC) compared to controls on PD 37. In addition, TM led to increased BW and lower TC after both cold and heat stress. Conclusively, our findings suggest that TM has a significant effect on the oxidative function of thermally stressed broilers.

**Keywords:** broiler; thermal manipulation; antioxidant; heat stress; cold stress

#### **1. Introduction**

The term broiler refers to any member of the red junglefowl subspecies, *Gallus gallus domesticus*, that has been reared for the purposes of meat production and consumption [1]. Constituting the largest standing avian population, the broiler has become a vital component of modern human nutrition as a result of the rapid industrialization of the poultry production process [2]. Since the mid-twentieth century, broilers have undergone intensive breeding in order to enhance their growth rates, meat yield, and feed conversion ratios [3]. However, these artificial selection pressures have been largely consumer-driven, focusing solely on improving certain commercially attractive parameters at the expense of immune, metabolic, and skeletal function [4]. As a result, the modern broiler has become increasingly susceptible to the effects of thermal and, in turn, oxidative stresses [5].

Due to impaired metabolic function and expensive energetics, broilers are especially vulnerable to heat stress, which occurs when the broiler is unable to adequately dissipate body heat to the environment [6]. Rising global temperatures have consolidated the threats of heat stress to the development and wellbeing of broiler chickens, and such increases in temperature are exacerbated during the hotter seasons [7]. In addition, broiler heat stress can be caused by certain stages of the poultry production process, especially during their transport from rearing to processing facilities [8]. Broilers subject to heat stress will have a lower body weight due to decreased feed intake, and their innate immune function will be impaired as a result of decreased immune organ weight [9]. In fact, it has been illustrated that heat stress results in oxidative stress in broilers, resulting in a number of adverse metabolic changes [10].

Heat stress is a major cause of oxidative stress in broilers, and oxidative damage deteriorates the appearance, flavor, and nutritional value of broiler meat [11]. Oxidative stress can be defined as the imbalance that occurs when the amount of reactive oxygen species (ROS) in an animal cell exceeds the latter's antioxidant capacity [12]. To prevent oxidative stress, several genes are involved in the maintenance of cellular homeostasis, including NADPH oxidase 4 (*NOX4*), superoxide dismutase (*SOD2*), and catalase [13,14]. Primarily expressed in renal and vascular cells, the *NOX4* gene is constitutively active and codes for an oxygen-sensing enzyme that can also play a role in antimicrobial defense [15–17]. If overexpressed, *NOX4* leads to oxidative stress due to its production of superoxide (O2 <sup>−</sup>) radicals and hydrogen peroxide (H2O2) molecules [18,19]. To prevent NOX4-associated oxidative stress from occurring, SOD2 and catalase act to dismutate O2 <sup>−</sup> and break down H2O2, respectively [20].

Unlike heat stress, cold stress in broilers has not been the subject of much research in the context of its relation to oxidative stress. Nonetheless, cold stress has been found to induce oxidative stress and modulate immune function in broilers while also increasing their susceptibility to necrotic enteritis and ascites development [21–25]. Moreover, cold stress was found to affect the thigh muscle of broilers more severely than the breast muscle, and it significantly reduced the feed intake and body weights of broilers but increased their feed conversion ratios [26,27]. As outdoor rearing systems gain more popularity, preventing cold stress will become increasingly costly to the poultry industry, and such costs often fluctuate depending on fuel prices, season, and existing heating systems [28].

To mitigate the damage caused by heat and cold stress, thermal manipulation (TM), which involves embryonic exposure to high or low temperatures, has been found to improve thermotolerance and enhance physiological parameters of broilers [13,29–33]. However, further research needs to be carried out in order to understand the effects of heat- and cold-induced oxidative stress in thermally manipulated broilers. Therefore, the main purpose of the present study was to investigate the effects of both cold and heat stress on the antioxidant defense mechanisms of thermally manipulated broiler chickens.

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

Ethical approval for all experimental procedures was obtained from the Animal Care and Use Committee at Jordan University of Science and Technology (approval # 16/3/3/418).

#### *2.1. Egg Procurement and Incubation*

Fertile Cobb eggs (n = 600) were obtained from local distributors based in Madaba, Jordan. Before incubation, eggs were thoroughly examined, and eggs were excluded if they displayed abnormality or damage (n = 69). The remaining eggs (n = 531) were then randomly divided into two groups, control (n = 266) and thermal manipulation (TM) (n = 265), and incubated in semi-commercial incubators (Masalles S.L., Barcelona, Spain). In the control group, eggs were incubated under standard conditions (37.8 ◦C and 56% relative humidity (RH)) throughout embryogenesis. In contrast, TM eggs were only incubated under standard conditions from embryonic days (ED) 1 to 9 and 19 to 21, as TM was applied from ED 10 to 18 by incubating the eggs at 39 ◦C and 65% RH for 18 h/day. On ED 7, candling was performed on each egg in order to exclude infertile and/or nonviable eggs.

#### *2.2. Hatchery Management*

On hatch day, the hatchability, which is the percentage of fertile eggs that hatch, was calculated according to the following equation: hatchability = (number of hatched chicks/total number of incubated eggs) × 100. Chicks were left in the incubator to dry for the first 24 h of their post-hatch life, after which they were transported to a special area designated for the field experiments. On post-hatch days (PD) 1 and 37, the cloacal temperatures (TC) and body weights (BW) were recorded, and the number of chicks that died within the whole field experimental period was noted. Dead chickens were histopathologically examined, but no significant obvious findings were reported. Before exposure to thermal stress, chicks were randomly distributed into their coops in groups of ten. In the first week, the temperature of the enclosures was kept at 33 ± 1 ◦C and was steadily reduced to 24 ◦C by the end of the third week. The RH during the rearing period was maintained within a range of 45%–52%. Water and appropriate feed were supplied to the chicks ad libitum during the whole field experiment period. On PD 8 and 20, chicks were vaccinated against Newcastle disease, and, on PD 15, the chicks were vaccinated against infectious bursal disease. The overall experimental design is illustrated in Figure 1.

**Figure 1.** Experimental design showing the three main phases: thermal manipulation, heat stress, and cold stress.

#### *2.3. Experiment 1: Post-Hatch Heat Exposure*

On PD 26, male chicks (n = 60) were randomly selected from each of the control and TM groups to be transported to the experimental room for the induction of heat stress. On PD 28, heat stress was induced by raising the temperature of the experimental room to 35 ◦C and 45%–52% RH until PD 35. During this period, male chicks (n = 60) from each of the control and TM groups were subject to normal conditions (not shown in Figure 1). The experimental and rearing rooms were located on the same floor in order to mitigate transport stress. After 0, 1, 3, 5, and 7 days of heat exposure, chicks (n = 8) were randomly chosen from the control and TM groups. The T<sup>C</sup> and BW of the chicks were recorded, after which they were euthanized in order to collect hepatic and splenic organs. Samples were snap-frozen on-site using liquid nitrogen, transferred to the laboratory, and stored at −80 ◦C.

#### *2.4. Experiment 2: Post-Hatch Cold Exposure*

On PD 32, chicks (n = 40) from each of the two incubation groups (control group and TM group) were randomly chosen and subdivided into four subgroups: control exposed to cold stress (CS), TM exposed to cold stress (TS), control exposed to normal conditions (CN) and TM exposed to normal conditions (TN). Cold stress was achieved by lowering the room temperature to 16 ◦C and 45%–52% RH from PD 32 to 37. At PD 37, BW and TC were recorded for chicks (n = 5) from each subgroup, and a number of chicks (n = 5) were humanely euthanized in order to collect the liver, spleen, heart, and breast muscle organs. Samples were snap-frozen on-site using liquid nitrogen, transferred to the laboratory, and stored at −80 ◦C.

#### *2.5. cDNA Synthesis*

The Direct-Zol™ RNA MiniPrep (Zymo Research, Irvine, CA, USA) was utilized alongside TRI Reagent® (Zymo Research, Irvine, CA, USA) in order to isolate total RNA from all the collected samples. The Biotek PowerWave XS2 Spectrophotometer (BioTek Instruments, Inc., Winooski, VT, USA) was employed to determine the quantity and quality of the samples, after which 2 μg of total RNA from each sample were inputted into the Superscript III cDNA Synthesis Kit (Invitrogen, Carlsbad, CA, USA) to synthesize cDNA.

#### *2.6. Primer Design and Relative mRNA Quantitation Analysis by Real-Time RT-PCR*

The primer sequences that were used for real-time RT-PCR analysis are listed in Table 1. Primers were taken from previous reports [34] and were designed using the PrimerQuest tool on the Integrated DNA Technologies website (Coralville, IA, USA) (https://eu.idtdna.com/pages) and the Nucleotide database on the NCBI (Bethesda, MA, USA) website (https://www.ncbi.nlm.nih.gov/nucleotide/). The QuantiFast SYBR® Green PCR Kit (Qiagen, Hilden, Germany) was utilized on a Rotor-Gene Q MDx 5 plex instrument (Qiagen, USA) according to the manufacturer's protocol. For the internal control, fold changes in gene expression were normalized against the 28S ribosomal RNA. Single target amplification specificity was ensured by the melting curve, and relative quantitation was calculated automatically by the software on the Rotor-Gene Q MDx 5 plex instrument.

**Table 1.** Primer sequences used for real-time RT-PCR analysis.


#### *2.7. Statistical Analysis*

IBM SPSS Statistics v23.0 (IBM, USA) was used for all statistical analyses performed in the current study. The chi-squared test was used to analyze hatchability and mortality rates. TC, BW, and the fold changes in mRNA levels of the catalase, *NOX4*, and *SOD2* genes are portrayed as means ± SD. An independent t-test compared between the control and TM groups with respect to several parameters at each time interval (PD 1, 3, 5, 7, 9, 11, 13, 15, 19, 22, 25, 28, 30, 33 and 35). Within the treatment group itself, two-way ANOVA was also used to compare between different parameters at specific time intervals after thermal stress. Statistical significance for parametric differences was set at 0.05.

#### **3. Results**

#### *3.1. E*ff*ect of Thermal Manipulation (TM) on Hatchability and Physiological Parameters of Broiler Chicks*

No significant effect was observed in either the mortality (control = 1.7; TM = 1.7) or the hatchability (control = 85.71; TM = 83.02) rates between the control and TM groups. However, TM led to significantly lower cloacal temperatures (TC) on PD 1 and 37 and to higher body weights (BW) on PD 37. However, no significant change was observed in hatchling BW (Table 2).

**Table 2.** Effects of embryonic thermal manipulation (TM) on post-hatch body weight (BW) and cloacal temperature (TC) of broiler chickens.


a,b within the same row, means ± SD with non-identical superscripts are significantly different.

#### *3.2. E*ff*ect of Post-Hatch Heat Stress on Physiological Parameters of Thermally Manipulated Broilers*

Table 3 illustrates the effects of heat stress for 7 days (PD 28 to 35) on TC, BW, and BW gain in the controls and TM broiler chickens. TM significantly decreased the mortality rate during post-hatch heat exposure (control = 12%; TM = 8%). Moreover, heat stress significantly increased the TC in both groups, but the T<sup>C</sup> of controls was significantly higher compared to TM chicks. On day 0 (PD 28) of heat stress, the BW of the TM group was significantly higher than that in controls. Similarly, the BW of controls was significantly lower compared to TM chicks on day 7 (PD 35) of heat stress, but the subgroups exposed to heat stress possessed significantly lower BW and BW gain compared to those exposed to normal conditions.

**Table 3.** Effects of post-hatch heat stress for 7 days (post-hatch day (PD) 28 to 35) on cloacal temperature (TC), body weight (BW), and BW gain in broiler chickens subjected to embryonic thermal manipulation and controls.


a–d within the same row, means ± SD with non-identical superscripts are significantly different.

#### *3.3. E*ff*ect of Post-Hatch Heat Stress on Antioxidant Enzyme mRNA Levels in Thermally Manipulated Broilers*

Figure 2 represents the effects of heat stress on the hepatic and splenic mRNA levels of certain antioxidant enzymes in broiler chicks subjected to embryonic TM. Tables S1 and S2 includes the mRNA levels of the same antioxidant genes for broiler chicks (TM and controls) not exposed to heat stress during the same timeframe.

**Figure 2.** Effect of post-hatch heat stress for 7 days (PD 28 to 35) on the mRNA levels of catalase, *NOX4*, and *SOD2* in the liver (**A**–**C**) and spleen (**D**–**F**) of TM broiler chicks (n = 5). \* within the same day, means <sup>±</sup> SD of TM and control chicks are significantly different. # within the control group, means <sup>±</sup> SD of non-identical days differ significantly. \$ within the TM group, means <sup>±</sup> SD of non-identical days differ significantly.

*Catalase*. On day 0 (PD 28) of heat stress, TM led to significantly lower catalase mRNA levels in the liver. In the control group, the hepatic mRNA levels of catalase were significantly higher after 5 (PD 33) and 7 (PD 35) days of heat stress compared to day 0, while, in the TM group, the level was significantly higher only after 1 day (PD 29) of heat exposure. The hepatic catalase mRNA level was significantly lower in TM chicks compared to controls after 5 (PD 33) and 7 (PD 35) days of heat stress.

The splenic mRNA level of catalase was not significantly different between TM and control chicks on day 0 (PD 28) of heat stress. However, the level was significantly lower in TM chicks compared to controls after 1 (PD 29) and 7 (PD 35) days of heat stress. Within the control group, the splenic mRNA level of catalase was significantly higher after 1 (PD 29), 5 (PD 33), and 7 (PD 35) days of heat exposure compared to day 0 (PD 28), whereas, in the TM group, the splenic mRNA level of catalase was significantly higher after 3 (PD 31) and 5 (PD 33) days (vs. day 0 (PD 28)).

*NOX4*. In the liver, the mRNA level of *NOX4* was not significantly different between the TM and control groups on day 0 (PD 28) of heat stress. In contrast, the *NOX4* mRNA level was significantly lower in TM chicks compared to controls after 3 (PD 31), 5 (PD 33), and 7 (PD 35) days of heat stress. Within the control group, the mRNA level was significantly increased after 3 (PD 31), 5 (PD 33), and 7 (PD 35) days of heat stress (vs. day 0 (PD 28)), but, in the TM group, the level did not significantly change during heat stress compared to day 0 (PD 28).

Similarly, the splenic mRNA level of *NOX4* was not significantly different between TM and control chicks on day 0 (PD 28) of heat stress. However, the level was significantly lower in the TM group compared to controls after 5 days (PD 33) of heat stress. Within the control group, the mRNA level of *NOX4* was significantly increased after 5 days (PD 33) of heat exposure in comparison with

day 0 (PD 28), while in the TM group, the level did not significantly change during heat exposure in comparison with day 0 (PD 28).

*SOD2*. The liver mRNA level of *SOD2* was not significantly different between the TM and control groups on day 0 (PD 28) of heat stress. However, *SOD2* levels were significantly lower in TM chicks compared to controls after 3 (PD 31), 5 (PD 33), and 7 (PD 35) days of heat stress. Within the control group, the hepatic mRNA level of *SOD2* was significantly higher after 7 days (PD 35) of heat stress compared to day 0 (PD 28), while, in the TM group, the level significantly increased only after 1 day (PD 29) of heat exposure (vs. day 0 (PD 28)).

The splenic mRNA level of *SOD2* did not significantly differ between the TM and control groups on day 0 (PD 28) of heat stress. Contrastingly, the splenic level was significantly lower in TM chicks compared to controls after 7 days (PD 35) of heat exposure. Within the control group, the splenic mRNA level of *SOD2* was significantly higher after 7 days (PD 35) of heat stress compared to day 0 (PD 28), whereas, in the TM group, the level was significantly higher only after 5 days (PD 33) of heat exposure (vs. day 0 (PD 35)).

#### *3.4. E*ff*ect of Post-Hatch Cold Stress on Physiological Parameters of Thermally Manipulated Broilers*

Table 4 represents the effects of cold stress for 5 days (PD 32 to 37) on TC, BW, and BW gain in thermally manipulated broiler chickens and controls. Application of TM significantly decreased the mortality rate during post-hatch exposure to cold stress (control = 5%; TM = 0). In contrast, cold stress did not significantly affect TC, but, in both the TC and TN subgroups, controls exhibited significantly higher T<sup>C</sup> compared to TM chicks. On day 0 (PD 32) of cold exposure, there was no significant change was observed in BW between the control and TM groups. After 5 days (PD 37) of cold stress, the BW of controls was significantly lower in comparison with control chicks exposed to cold stress. Furthermore, cold stress significantly decreased the BW gain in both the control and TM chicks, although the weight gain was significantly lower in controls compared to TM chicks.


**Table 4.** Effect of post-hatch cold stress for 5 days (post-hatch day (PD) 32 to 37) on cloacal temperature (TC), body weight (BW), and body weight gain in broiler chickens subjected to embryonic thermal manipulation and controls.

a–d within the same row, means ± SD with non-identical superscripts are significantly different.

#### *3.5. E*ff*ect of Post-Hatch Cold stress on mRNA Levels of Antioxidant Enzymes in Thermally Manipulated Broilers*

Figure 3 represents the effects of cold stress on the mRNA levels of antioxidant enzymes in the liver, spleen heart and breast muscle of broiler chickens subjected to embryonic thermal manipulation.

**Figure 3.** Effect of post-hatch cold stress for 5 days (PD 32 to 37) on the mRNA levels of catalase, *NOX4*, and *SOD2* in the liver (**A**–**C**), spleen (**D**–**F**), heart (**G**–**I**), and breast muscle (**J**–**L**) in TM broiler chickens (n <sup>=</sup> 5). a–d means <sup>±</sup> SD with non-identical superscripts are significantly different.

*Catalase*. TM did not significantly change the cardiac, hepatic, and muscular mRNA levels of catalase in chicks kept under normal environmental temperatures. However, TM significantly decreased the catalase mRNA level in the spleen. Regarding those chicks exposed to cold stress, the TM group possessed a significantly lower mRNA level of catalase in the liver, spleen, and heart compared to controls.

*NOX4*. In the chicks of the TN subgroup, TM chicks possessed significantly lower splenic, hepatic, and cardiac mRNA levels of *NOX4* compared to controls. Despite this, muscular NOX4 mRNA levels were significantly higher in the TM chicks. Similar results were observed in the chicks exposed to cold stress.

*SOD2*. No significant changes were seen in the hepatic and muscular *SOD2* mRNA levels between the TM and control groups exposed to normal conditions. However, the splenic and cardiac levels of *SOD2* mRNA were significantly higher in controls compared to TM chicks. After cold exposure, the cardiac, hepatic, and splenic mRNA levels of *SOD2* were significantly higher in controls compared to TM chicks, but the level in breast muscle was significantly higher in the TM group.

#### **4. Discussion**

Oxidative damage is caused by excess reactive oxygen species (ROS), such as superoxide (O2 −) and hydrogen peroxide (H2O2), which are a necessary product of aerobic metabolism [35]. Heat stress is a major cause of oxidative damage in poultry, and it is associated with a modulation in the expression of antioxidant genes, including catalase, *NOX4*, and *SOD2* [10,13]. Thermal manipulation (TM) has often been suggested as a viable method of improving the acquisition of thermotolerance in heat-stressed broilers [31,36–38]. However, the effects of TM and subsequent heat challenge on broiler antioxidant capacity has not been extensively explored. Similarly, a dearth of information exists with regard to the effects of cold stress on broiler gene expression, especially within the context antioxidant gene expression. The objective of the present study was two-fold: it aimed to ascertain the effects of embryonic TM on broilers under conditions of post-hatch heat stress as well as cold stress.

During heat stress, the behavior of broilers is altered as they attempt to decrease their body temperature (TC), resulting in myriad negative effects on performance [39]. In the current study, TM was found to result in significantly lower T<sup>C</sup> on post-hatch days (PD) 1 and 37 and higher body weights (BW) on PD 37 compared to controls. Correspondingly, it has often been reported that TM treatments significantly increased broiler BW [37,40,41] and improved their abilities to regulate their T<sup>C</sup> in periods of heat challenge [31,32,42]. Lower TC during heat stress also improved feed conversion ratios in TM broilers, and it has been suggested that the lower TC in TM broilers is due to slower metabolic rates as a result of the TM treatment [43].

With regard to post-hatch heat stress, our findings show that TM chickens had significantly decreased mortality rates and BT as well as increased BW compared to controls. Heat stress has been extensively reported to affect broiler physiological parameters [8,44,45]. On a similar note, cold-stressed TM chickens had significantly lower mortality rates than cold-stressed controls. Previously, TM has been found to reduce the mortality rates of broilers during heat challenge [46,47]. Contrastingly, one study reported that heat-stressed TM broilers experienced higher mortality rates than their control counterparts [48]. These differences in findings may be attributed to the fact that there is no one single type of TM treatment, and different studies employ different periods and conditions of TM.

The catalase enzyme is found in the majority of aerobic organisms as well as in some obligate anaerobes [49]. Catalase is responsible for the breakdown of hydrogen peroxide (H2O2) into oxygen and water, thereby preventing oxidative damage from occurring in a cell [50]. In the present study, catalase expression was significantly modulated in heat- and cold-stressed TM and control chickens. In fact, heat stress resulted in decreased hepatic and splenic catalase expression in TM chickens compared to controls, while cold stress led to significantly lower cardiac, hepatic, and splenic catalase expression levels in the TM group. Compared to controls, heat-stressed TM chickens were previously reported to exhibit decreased catalase mRNA levels [13]. Moreover, broilers were found to exhibit higher levels of catalase activity during acute heat stress, but this antioxidant capacity decreased with age [51]. In female broilers, cardiac catalase activity was reduced after cold stimulation [52].

To maintain homeostasis, the NOX4 enzyme is heavily involved in the oxygen-sensing process, the latter of which causes it to generate significant amounts of ROS [17]. Additionally, *NOX4* over-expression is often associated with oxidative stress in a number of different organs [18,53,54]. In avian muscle cells, ROS production during heat stress and subsequent oxidative damage has been tentatively attributed to *NOX4* up-regulation [55]. In the present study, hepatic and splenic *NOX4* expression levels were significantly lower in heat-stressed TM chickens compared to controls. Similarly, after cold stress, TM chickens exhibited decreased cardiac, hepatic, and splenic but increased muscular *NOX4* expression than that in controls. Hepatic *NOX4* mRNA expression was previously reported to be lower in TM chickens exposed to heat stress compared to controls [13]. In cultured avian cells, heat stress was found to upregulated *NOX4* mRNA expression [55].

The SOD2 enzyme functions to transform the superoxide (O2 −) radical into hydrogen peroxide and water, and it plays an important cytoprotective role against oxidative stress [56]. Our findings indicate that both heat and cold exposure led to generally decreased *SOD2* expression in several organs. Compared to controls, heat-stressed TM chickens displayed lower hepatic and splenic *SOD2* expression levels, while cold-stressed TM chickens showed decreased cardiac, hepatic, and splenic *SOD2* expression. Like *NOX4*, however, muscular *SOD2* expression levels were higher in cold-stressed TM chickens compared to their control counterparts. A previous study found that hepatic *SOD2* expression and enzymatic activity were decreased in TM chickens exposed to heat stress [13]. In contrast, another study found that *SOD* mRNA levels in two broiler strains (Cobb and Hubbard) were unaffected by heat stress [57]. Additionally, *SOD2* levels remained unchanged in avian cell cultures exposed to heat stress [55]. The present findings may suggest that lower levels of oxidative *NOX4* expression may lead to lower expression of the anti-oxidative catalase and *SOD2* genes.

Interestingly, mRNA expression levels of the catalase, *NOX4*, and *SOD2* genes in the breast muscle differed from those in the heart, liver, and spleen in cold-stressed TM chickens. Such inter-organ variation in gene expression is to be expected, as expression varies to a larger degree between organs of a single species than between different species [58]. However, in broilers, the breast muscle in particular has been subject to rapid changes in size and conformation over the past few decades due to the artificial selection pressures applied by the commercial poultry industry [59]. This has resulted in a number of abnormalities and myopathies of the breast muscle that is estimated to affect up to 90% of broilers worldwide [60–62]. In fact, broiler breast muscle cells were suggested to constantly undergo hypoxic stress, as the transcriptional profiles of non-stressed broiler breast muscle and heat-stressed layer breast muscle were similar [63].

A number of strengths can be found in the current study. All samples were taken from male Cobb chicks in order to reduce inter-strain and inter-sex genetic variation. Moreover, any non-experimental stress was minimized by ensuring that the rearing and experimental rooms were in close proximity to one another. However, there are some limitations of the present study. Firstly, the effect of TM on the developmental parameters of broiler embryos was not investigated, requiring future research. Secondly, the oxidation levels of lipids, proteins, and DNA in different tissues must still be measured in order to ascertain the final balance of catalase, *NOX4*, and *SOD2* expression. Lastly, the exact impact of TM on embryonic mortality was not considered, which mandates future lines of research in this context.

#### **5. Conclusions**

Our findings indicate that TM at 39 ◦C and 65% RH for 18 h/day from days 10 to 18 of embryonic development might result in positive long-lasting effects on broiler antioxidant capacity. Future research should focus on the effects of TM and subsequent thermal challenge on various types of broiler muscle, as the expression dynamics of the breast muscle was found to differ from those of other organs.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/1/126/s1, Table S1: effects of heat stress on the hepatic mRNA levels of certain antioxidant enzymes in broiler chickens subjected to embryonic TM; Table S2: effects of heat stress on the splenic mRNA levels of certain antioxidant enzymes in broiler chickens subjected to embryonic TM.

**Author Contributions:** Conceptualization, K.M.M.S., A.H.T. and M.B.A.-Z.; methodology, K.M.M.S.; software, K.M.M.S. and M.B.A.-Z.; validation K.M.M.S. and M.B.A.-Z.; formal analysis, K.M.M.S. and M.B.A.-Z.; investigation, K.M.M.S.; resources, M.B.A.-Z.; data curation, K.M.M.S. and M.B.A.-Z.; writing—original draft preparation, K.M.M.S. and A.H.T.; writing—review and editing, K.M.M.S., A.H.T. and M.B.A.-Z.; visualization, K.M.M.S., A.H.T. and M.B.A.-Z.; supervision, M.B.A.-Z.; project administration, M.B.A.-Z. and K.M.M.S.; funding acquisition, M.B.A.-Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deanship of Research/Jordan University of Science & Technology, grant number 44/2019.

**Acknowledgments:** The authors would like to thank Eng. Ibrahim Alsukhni for his excellent technical assistance and valuable comments.

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

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Animals* Editorial Office E-mail: animals@mdpi.com www.mdpi.com/journal/animals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18