Identifications of immune-responsive genes for adaptative traits by comparative transcriptome analysis of spleen tissue from Kazakh and Suffolk sheep

Aridity and heat are significant environmental stressors that affect sheep adaptation and adaptability, thus influencing immunity, growth, reproduction, production performance, and profitability. The aim of this study was to profile mRNA expression levels in the spleen of indigenous Kazakh sheep breed for comparative analysis with the exotic Suffolk breed. Spleen histomorphology was observed in indigenous Kazakh sheep and exotic Suffolk sheep raised in Xinjiang China. Transcriptome sequencing of spleen tissue from the two breeds were performed via Illumina high-throughput sequencing technology and validated by RT-qPCR. Blood cytokine and IgG levels differed between the two breeds and IgG and IL-1β were significantly higher in Kazakh sheep than in Suffolk sheep (p < 0.05), though spleen tissue morphology was the same. A total of 52.04 Gb clean reads were obtained and the clean reads were assembled into 67,271 unigenes using bioinformatics analysis. Profiling analysis of differential gene expression showed that 1158 differentially expressed genes were found when comparing Suffolk with Kazakh sheep, including 246 up-regulated genes and 912 down-regulated genes. Utilizing gene ontology annotation and pathway analysis, 21 immune- responsive genes were identified as spleen-specific genes associated with adaptive traits and were significantly enriched in hematopoietic cell lineage, natural killer cell-mediated cytotoxicity, complement and coagulation cascades, and in the intestinal immune network for IgA production. Four pathways and up-regulated genes associated with immune responses in indigenous sheep played indispensable and promoting roles in arid and hot environments. Overall, this study provides valuable transcriptome data on the immunological mechanisms related to adaptive traits in indigenous and exotic sheep and offers a foundation for research into adaptive evolution.


Identification of DEGs in spleen.
A Venn diagram was constructed showing expressed genes in the two libraries. The TPM method was used to identify expression of 19,724 and 17,997 reference genes in the Kazakh and Suffolk sheep libraries, respectively, and 15,465 shared genes (Fig. 3). A total of 1,158 significantly DEGs were identified between the two libraries, of which 246 genes were up-regulated and 912 genes were down-regulated in Suffolk sheep (Fig. 4). The significant DEGs between two libraries were listed in Supplementary Table S1.
Functional enrichment of DEGs. To gain insights into biological implications, the enrichment of DEGs in terms of GO were tested in this study (Fig. 5). The summary of all significantly enriched GO assignments of DEGs are shown in the results. In the cellular component category, a significant percentage of clusters were Table 1. Contents of IgG and 8 Cytokines in serum of Suffolk and Kazakh sheep*. *Means ± SEM marked with different lower-case letter either a or b in the same line were significant at P < 0.05.

Kazakh (n = 94) Suffolk (n = 56) Rams (n = 65) Ewes (n = 94)
IgG (mg/ml) 3   www.nature.com/scientificreports/ assigned to membrane (42.4%), membrane part (20.3%), organelle (17.4%), and organelle part (16.8%). The main functional groups of DEGs in molecular function category were assigned to catalytic activity (9.7%) and binding (1.3%). In the biological process category, the terms were related to immune system process (7.1%),    www.nature.com/scientificreports/ metabolic process (2.4%), cellular process (2.4%), response to stimulus (2.6%), multi-organism process (1.6%), and cell killing (1.1%). KEGG pathway analysis was conducted to categorize and annotate 384 DEGs. Further comparisons in the KEGG pathway database identified four significantly enriched pathways that were associated with the immune system (Table 3), including hematopoietic cell lineage, natural killer cell mediated cytotoxicity, complement and coagulation cascades and intestinal immune network for IgA production. The most significant difference was in the hematopoietic cell lineage. Ten hematopoietic cell lineage genes (IgD, DQA2, DQB, DOB, CD3D, MS4A1, CSF1R, CD2, ITGA2B, and CD3E) were significantly enriched by comparing Suffolk with Kazakh sheep (P value = 0.0004). The genes DQA2 (HLA-DR), DQB (HLA-DR), CD3D (CD3), CD3E (CD3), and CD2 were significantly down-regulated (Fig. 6). Eight natural killer cell-mediated cytotoxicity genes (IgD, GRB10, GZMB, FAS, FCER1G, TYROBP, LCP2, and CD48) were significantly enriched by comparing Suffolk with Kazakh sheep (P value = 0.0075). The genes GZMB (Granzyme), FAS, FCER1G (FcεR1y), LCP2 (SLP-76), and CD48 were significantly down-regulated (Fig. 7). Four complement and coagulation cascade genes (C1s, HP, C3, and C4BPA) were significantly enriched by comparing Suffolk with Kazakh (P value = 0.0112). The genes C1s (C1qrs) and C3 were significantly up-regulated, and gene C4BPA (C4BP) was significantly down-regulated (Fig. 8). Four intestinal immune networks for IgA production genes (IgD, DQA2, DQB, and DOB) were significantly enriched by comparing Suffolk with Kazakh sheep (P value = 0.0158). The genes DQA2 (MHC) and DQB (MHC) were significantly down-regulated (Fig. 9).     Table 4 show that a total of 21 immune-responsive genes were identified in the spleen according to their KEGG assigned pathways including five up-regulated genes and 16 down-regulated genes in the Suffolk sheep spleen. Furthermore, two up-regulated genes (DQA2 and C1s) and three down-regulated genes (DQB, DOB and GZMB) showed the most significant differential expression with an absolute value of log2 ratio . Differentially expressed genes involved in natural killer cell mediated cytotoxicity in the spleen of two sheep breeds. The blue-labelled genes were significantly down-regulated (p < 0.05). The green-labelled genes were expressed in both sheep spleens. KEGG pathway annotation was conducted using the DAVID website (http://david .abcc.ncifc rf.gov/) and searching the KEGG database. Protein-protein interaction network of DEGs. In this study, the STRING online tool was used to construct the PPI of 21 DEGs associated with the immune system. In total, 16 immune-responsive genes were built into two interaction networks using Cytoscape 2.8.2 software (Fig. 10). Of these 12 down-regulated genes (TYROBP, Figure 8. Differentially expressed genes involved in complement and coagulation cascades in the spleen of two sheep breeds. The red-labelled genes were significantly up-regulated (p < 0.05). The blue-labelled genes were significantly down-regulated (p < 0.05). The green-labelled genes were expressed in both sheep spleens. KEGG pathway annotation was conducted using the DAVID website (http://david .abcc.ncifc rf.gov/) and searching the KEGG database.

Discussion
Adaptation can be defined as the level of tolerance to survive and reproduce under extreme living conditions 21 .
Sheep are an excellent model for researching the genetic mechanisms underlying the adaptation of livestock to extreme environments 13 , because they have spread and adapted to a wide range of agroecological conditions  Figure 10. The protein-protein interaction networks of the significantly differentially expressed genes associated with the immune system. The sketch represents the interaction network constructed for proteincoding DEGs. Green nodes indicate down-regulated genes, and red nodes indicate up-regulated genes. The lines indicate the interactions between the gene products. Research shows that most adaptive traits are at least moderately heritable 23 , suggesting that animal breeding to improve adaptation is feasible. In general, immune responses can be regulated by the immune factors related to animal adaptation to extreme environments 4 . According to previous studies, comparisons between exotic and local breeds of European sheep showed that the local Lacaune and Merino breeds had more resistance to natural infections with T. circumcincta than the imported, highly prolific Romanov sheep from the south of France 24 . Red Maasai sheep were more resistant (low FEC) and more resilient (high PCV) to gastro-intestinal (GI) nematode parasites than Dorper sheep in African sub-humid coastal environment 25 . Pelibuey sheep have been found to regulate body temperature after heat stress more efficiently than Suffolk sheep via a mechanism related to HSP-70 5 . On the Qinghai-Tibetan Plateau, Tibetan sheep have introgressed genomic regions related to the oxygen transportation system, sensory perception and morphological phenotypes from Argali (Ovis ammon), including HBB and RXFP2 genes, which have strong signs of adaptive introgression 12 . Previous research on sheep has mainly focused on specific immunity against parasitic infections and genetic mechanisms of adaption to different environments; however, research on the relationship between adaptation and immunity is lacking. In this study, immunological differences and immune-responsive genes for adaptive traits from indigenous and exotic sheep breeds were identified and compared. The Kazakh and Suffolk sheep appeared to be healthy as evidenced by histomorphology observation of spleens. ELISA results showed that the levels of eight cytokines and IgG were higher in the blood of Kazakh compared with Suffolk sheep, and the levels of IgG and IL-1β were significantly higher in Kazakh than in Suffolk sheep. However, they grow slowly and have excessive fat-deposition. Suffolk sheep grow quickly and have excellent meat quality, but do not adapt to harsh environments and have lower levels of IgG and cytokines. Thus, Kazakh sheep may have better immunity and adapt better to arid and extremely hot environments compared with Suffolk sheep. These results indicate that indigenous adapted breeds tend to be small with poor production characteristics, whereas high-performing breeds often have poor disease-resistance characteristics. Native or indigenous livestock breeds are well adapted to their regions and are of historical and economic importance 26 , but intensive farming practices widely use improved or cross-bred animals for increased productivity and profit, which results in significant declines of native breeds 27 . Native breeds may be more adapted to a more stringent environment and climate, such as that of high mountains and northern coastal areas that are relatively marginal for agriculture 26 . Therefore, preservation of genetic diversity of native breeds can play key evolutionary roles in smart agriculture practices to dealing with global climate change and sustainable farming in the long term 26 .
To better understand the molecular mechanisms related to the immunity of adaptive traits in sheep, RNA-Seq techniques were employed and transcriptome sequencing for the DEG identifications was performed in spleen tissues from Suffolk and Kazakh sheep. A total of 1,158 significantly DEGs were identified between the two libraries, thus providing a valuable resource for immunogenetics and genomic research on sheep spleen tissue. Furthermore, 25 GO terms of DEGs were enriched in the cellular components, biological processes and molecular functions. Immune-related proteins were included in the biological processing category with higher percentages in immune system processes and response to stimuli. Within the molecular function category, the spleen transcripts were mainly associated with catalytic activity and binding, which represent some immunerelated processes. The GO analysis followed a similar pattern to that obtained using pyrosequencing for defense mechanisms by injection with viral stimuli to increase the expression level of immune-related genes in Turbot 28 . These percentages suggest that a few immune-related genes were identified in this study.
Previous studies have reported that signaling pathways, such as toll-like receptor signaling pathway, Imd pathway, JAK/STAT pathway, JNK pathway 29 , TNF signaling pathway, B cell receptor signaling pathway, T cell receptor signaling pathway 30 , and NF-κB signaling pathway 31 , play an important role in immunity. In this study, comparison with the KEGG pathway database, showed that the toll-like receptor signaling pathway, JAK/STAT pathway, B cell receptor signaling pathway and T cell receptor signaling pathway were not significantly enriched; however, four immune system pathways were significantly enriched. Hematopoietic cell lineage refers to the differentiation and developmental processes of hematopoietic stem cells (HSCs), which divide into common myeloid progenitor cells by stimulation with IL-1, IL-3, IL-6, GM-CSF, and SCF. Cellular cytotoxicity is an important effector mechanism of the immune system to combat viral infections and is majorly mediated by www.nature.com/scientificreports/ cytotoxic T-cells and natural killer (NK) cells 32 . The complement system is a mediator of innate immunity and a nonspecific defense mechanism against pathogens. Complement and coagulation cascades may be implicated in the protective strategy of the brain for resistance to cold and environmental adaptation in the Himalayan marmot 33 . The intestine is the largest lymphoid tissue in the body, and intestinal immunity can generate many noninflammatory IgA antibodies, which are the first line of defense against microorganisms. Multiple cytokines, including TGF-β, IL-10, IL-4, IL-5, and IL-6, are necessary factors for IgA class switching and the terminal differentiation process in B cells 34 . The protein-protein interaction network of DEGs, the hematopoietic cell lineage, natural killer cell mediated cytotoxicity, and the intestinal immune network for IgA production pathway include 12 genes that coregulate the immune system in the spleen. The complement and coagulation cascades pathway includes four genes that regulate the immune system alone. These results indicate that hematopoietic cell lineage, natural killer cell mediated cytotoxicity, and the intestinal immune network for IgA production participate together in the regulation of the immune response. Results of this study suggest that the complement and coagulation cascade pathway likely plays an important role in the biological process of immune response and adaptation to arid and hot environments. DQA2, C1s, DQB, DOB, and GZMB genes of the 21 DEGs associated with immunity were listed by the threshold of FDR ≤ 0.001 and fold change ≥ 2. DQA2, DQB and DOB showed significant enrichment in hematopoietic cell lineage and intestinal immune network for IgA production. DQA and DQB located in the MHC class IIa region, and DOB, located in MHC class IIb region, control self/non-self-recognition of the immune system 35 , and play important roles in immune response to adaptive trait. Granzyme B (GZMB) is significantly enriched in natural killer cell-mediated cytotoxicity and play an important role in the ability of NK cells and CD8+ T cells to kill their targets and modulate inflammation 36 . GZMB levels in plasma might reflect the degree of pruritus and dermatitis in patients with Alzheimer's disease 37 . C1s is significantly enriched in complement and coagulation cascades, and can initiate the classical signaling pathway of complement activation and prevent of the formation and precipitation of large immune complexes 38 . In the current study, some genes were significantly down-regulated in Suffolk sheep, but up-regulated in Kazakh sheep, which had higher levels of blood IgG and cytokines, indicating that these genes might be conducive to adapting to arid and hot environment. DQB, DOB, GZMB, CD3D, MS4A1, FCER1G, CSF1R, TYROBP, LCP2, CD48, CD2, and CD3E (the 12 genes in the PPI network of hematopoietic cell lineage, natural killer cell-mediated cytotoxicity and the intestinal immune network for IgA production pathway), play roles in immune regulation. C1s, HP, and C3 in PPI network of the complement and coagulation cascades were significantly up-regulated in Suffolk sheep, but down-regulated in Kazakh sheep, suggesting that living in harsh environments for a long time might be not conducive to growth and development due to a suppressed immune pathway. In summary, these findings, combined with previous results, clearly indicate that many spleen-specific DEGs play roles in immune function in indigenous and exotic sheep breeds with different adaptabilities. Future studies are needed to determine functional characterization of the immune-response genes associated with adaptation in sheep.

Conclusion
In this study, the first comprehensive transcriptome sequencing and gene expression profiling analysis in spleen tissues based on content differences of IgG and cytokines between indigenous Kazakh sheep and exotic Suffolk sheep were performed. The results revealed the immunological and genetic differences associated with adaptation between indigenous and exotic sheep, with higher levels of blood IgG and cytokines in indigenous Kazakh sheep and 1,158 significantly differentially expressed genes in the spleens between two sheep breeds. Twenty-one immune-responsive genes were identified as spleen-specific genes associated with adaptive traits. Four pathways (hematopoietic cell lineage, natural killer cell mediated cytotoxicity, complement and coagulation cascades, and intestinal immune network for IgA production pathway) and up-regulated genes (including DQB, DOB, and GZMB) associated with immune response in indigenous sheep were helpful for adapting to harsh environments. Transcriptome and gene expression profiling data are valuable resources for future comparative genomic analysis on immunological mechanisms of adaptive trait and offer a foundation for research into genetic adaptive evolution.

Materials and methods
Animals and climatic conditions. The study was carried out in compliance with the ARRIVE guidelines. Samples. Peripheral blood was collected by jugular venipuncture. Serum was collected after centrifugation (2500 rpm for 10 min), and transported to the laboratory at 4 °C for ELISA assays. Three male Suffolk sheep and three male Kazakh sheep were slaughtered, and samples of the parenchyma of the spleen were collected. One part of each sample was immediately frozen in liquid nitrogen and stored at − 80 °C for subsequent total RNA extraction and another part was used for paraffin sections. www.nature.com/scientificreports/ Hematoxylin-eosin (HE) staining of spleen tissues. The spleen tissues were treated for histological imaging using standard H&E staining procedures. Spleens were fixed in freshly prepared 4% paraformaldehyde solution in PBS for 18 h. Spleens were then washed with PBS, dehydrated in gradient ethanol, made transparent with xylene, and embedded in paraffin wax. The wax blocks were fixed on a slicer and continuously sectioned at a thickness of 4 μm using a Leica CM1900 cryostat (Leica Microsystems, Wetzlar, Germany). Sections were dried at 37 °C, stained with a neutral resin, and observed under a microscope.
RNA extraction and quality assessment. Total RNA was isolated from spleen tissues using TRIzol . The mRNA was enriched from total RNA using oligo (dT) magnetic beads (Invitrogen, Carlsbad, CA, USA). The purified mRNA was pooled, and fragmentation buffer was added to disrupt the mRNA into short fragments of about 200 bp. A random hexamer primer was added to synthesize the first strand cDNA using the short fragments as templates. Buffer, dNTPs, RNase H, and DNA polymerase I were sequentially added to synthesize the second strand cDNA. The double-stranded cDNA was subsequently purified using a QiaQuick PCR extraction kit and resolved with elution buffer for end repair and addition of a poly (A) tail. Illumina sequencing adaptors were then ligated to the short fragments. The suitable fragments were purified by agarose gel electrophoresis and enriched using PCR amplification as templates. The average insert size for the paired-end libraries was 200 bp. Finally, two paired-end cDNA libraries were constructed and used for sequencing analysis on the Illumina HiSeq 2000 system. The imaging data output was transformed by base calling into raw reads and stored in FASTQ format. All clean reads went through a stringent filtering process. For obtaining clean reads, the program removes raw reads, reads with adaptors, and reads in which unknown bases were more than 10%, and low-quality reads (i.e. the percentage of the low quality bases of quality value ≤ 5 was more than 50% in a read). Transcriptome de novo assembly was carried out with short reads using Trinity software packages 39 . After combined reads with overlaps formed contigs, Trinity was used to connect the contigs to obtain sequences that could not be extended on either end. These contigs are subjected to further processing of sequence clustering to form longer sequences without N. Such sequences were defined as unigenes. Basically, unigenes are a collection of expressed sequences that are aligned or located to same position on the genome. The unigenes were annotated using the NR, NT, SwissProt, KEGG, COG, and Gene ontology (GO) databases and the number of annotated unigenes counted.

Expression profiling and analysis of DEGs. Digital gene expression profiling (DGE) was performed at
Beijing Genomics Institute (BGI, Shenzhen, China) using the Illumina Gene Expression Sample Prep Kit and Illumina sequencing chip (Illumina, San Diego, CA, USA) according to the manufacturer's instructions.
The procedure was modified from a published method 40 . The Illumina Cluster Station and Illumina HiSeq 2000 System were used. Experimental processes were as follows: mRNA was purified from 6 µg total RNA by using Oligo(dT) magnetic beads and the first and second-strand cDNA were synthesized using Oligo(dT) as primer. The 5′ ends of tags generated the NlaIII endonuclease site and the bead-bound cDNA was subsequently digested with restriction enzyme NlaIII, which recognized and cut off the CATG sites. The fragments separated from the 3′ cDNA fragments and connected to Oligo(dT) beads were washed away and the Illumina adaptor 1 was ligated to the sticky 5′ end of the digested bead-bound cDNA fragments. The junction of the adaptor 1 and CATG site was the recognition site of MmeI, which was cut at 17 bp downstream of the CATG site, producing tags with adaptor 1. After removing the 3′ fragments via magnetic bead precipitation, the Illumina adaptor 2 was ligated to the 3′ ends of the tags to acquire tags with different adaptors of both ends to form a tag library. After 15 cycles of linear PCR amplification, 105 bp fragments were purified using 6% TBE PAGE gel electrophoresis. After denaturation, the single-chain molecules were fixed onto the Illumina sequencing chip. Each molecule grew into a single-molecule cluster sequencing template through situ amplification. Then, four types of nucleotides labeled using four colors, were added, and the sequencing by synthesis (SBS) method was performed. Each tunnel generated millions of raw reads with sequencing lengths of 49 bp.
Sequencing-received raw image data were transformed via base calling into raw data, which were stored in FASTQ format. Raw sequences were transformed into clean tags by a stringent filtering process: 3′ adaptor sequences, empty reads, low quality tags with unknown sequences (N), and tags with a copy number of 1 were Protein-protein interaction network. The STRING Ver.11.0 (http://strin g-db.org/) database was used to construct a potential protein-protein interaction (PPI) network between proteins encoded by DEGs related to immunity. The PPI data obtained from the STRING database was imported into the Cytoscape software, and the nodes were used as a network center node.
Validation of RNA-Seq data using quantitative real-time PCR (RT-qPCR). Eleven candidate genes were selected and detected by using RT-qPCR. The primers for the 11 genes were designed using Primer-BLAST 49 ( Table 5). The amplification product of the expected size for each gene was verified by electrophoresis on a 2.0% agarose gel and confirmed by sequencing. A quantity of 2 μg total RNA was used to synthesize first strand cDNA using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa Biotech, Dalian, China), according to the manufacturer's instructions. The RT-qPCR reactions were performed in a 20 μL volume containing 10 μL 2 × LightCycler 480 SYBR Green I Master (Roche Diagnostics, Mannheim, Germany), 2 μL cDNA, 0.5 μM of each primer, and 7 μL ddH2O. The amplification conditions were as follows: 95 °C for 5 min of initial stage, followed by 40 cycles at 95 °C for 10 s, annealing temperature for 10 s, and 72 °C for 20 s, performed on the LightCycler 480II system (Roche Diagnostics GmbH, Mannheim, Germany). The GAPDH gene was used as an internal control. All RT-qPCR reactions were carried out in triplicate, with both technical and biological replicates. RNA samples from six Suffolk sheep and six indigenous Kazakh sheep were used for this analysis. Values of RT-qPCR were calculated by 2 −ΔΔCt method and data were analyzed with IBM SPSS 24.0 50 .

Ethics declarations.
The Guidelines for the Care and Use of Laboratory Animals were carefully followed during this study, which received approval from the Experimental Animal Care and Use Committee of Xinjiang Academy of Agricultural and Reclamation Sciences (Shihezi, China, approval number: XJNKKXY-AEP-039, 22 January 2012), and was approved by the Institutional Animal Care and Use Committee of University of Hawaii  TTT GAG TCT GGG GA, R: GGT AGC CGA CGT GAC AGT AG  61   SAA  F: AGC AGG TCA CAG CCA AAG GAT, R: TCC ACA TGT CTT TAG CCC CTT GAC 61   C3  F: GCA ATC AGG CAG CGA CAT GG, R: CAC GTA CAC CAT GAG GTC GAA  61   CD48  F: ACC CTG TGT ACA AGC CTC CT, R: TCT GTC CTG GTT AAG AGC CG  62   TYROBP  F: GAC TGT GGG TGG TCT CAG C, R: AGT ACA CAG CCA GAG CGA TG  62   FCER1G  F: TAC TGC CGA CTC AAG CTG C, R: AGC TAC TGT GGC GGT TTC TC  50   IFITM1  F: TGT GGT CCC TGT TCA  www.nature.com/scientificreports/ under the protocol No. 10-1039-2 and the Animal Welfare Assurance No. A3423-01. All efforts were made to follow the animal care and sampling procedures to minimize animal suffering.