The effect of crossbreeding by the Montbeliard cattle on the transcriptome of the Sistani cattle

Objective: Transcriptome study of the desired crossbreed and comparing it with the pure parental population is a method for evaluating and selecting cross with appropriate products, yet disease resistance, high tolerance to environmental stresses, and compatibility with the climate of the breeding tract. RNA-Sequencing (RNA-seq) technology was employed to investigate the effect of crossbreeding by the Montbeliard cattle on the “Sistani cattle” transcriptome. Blood samples were collected from 90 cows (45 pure and 45 crossbreeds). Total RNA was extracted in the Media Teb Gene Laboratory of Iran. The mRNA selection, library preparation, and sequencing were performed at BGI Shenzhen, China on an Illumina HiSeq 2500 sequencer. After assembling the transcripts with the bovine reference genome, Cuffdiff was used to detect differential gene expression in purebred and its crossbred. Results: In total, 45 signi ﬁ cantly differentially expressed genes were identi�ed (false discovery rate q< 0.05) between purebred Sistani and “Sistani×Montbeliard” crossbreed. These genes are reported to have a role in immune response, calving, fertility, and resistance to mastitis. Besides, they may have a role in different levels of heat tolerance and disease resistance. Therefore, Sistani×Montbeliard crossbred might be suitable for the climatic conditions of the Sistan region in Iran.


Introduction
There has been increased attention to crossbreeding between the indigenous cattle populations (Bos indicus) and the imported breeds (Bos Taurus) in rural areas in Iran, mainly to increase milk production in recent years. Crossbreeding is suggested as a shortcut to improving performance characteristics through adapting to harsh tropical conditions [19,20].
Sistani cattle (Bos indicus), as the only known meat-type breed in Iran, has special genetic features adapted to areas not suitable to crop agriculture [24]. In recent years, most of the area's breeders and villagers are reluctant to continue the pure breeding of these cattle such that there has been a tendency to crossbreed with foreign breeds such as Montbeliard, Simmental, and Holstein. Due to shrinking the free grazing land since producing concentrate feed is not economically justi ed, rearing local cattle breed is not pro table to farmers [17].
Evaluating crossbreeding schemes needs to collect the performance records and some adaptive indices that are hard to measure. The transcriptome pro ling and the comparison between pure parental animals and the crossbred may be an alternative method [24]. The best genetic combination from crossbreeding can be decided based on differential gene expression analysis and the desired crossbreeding goals (high productivity or higher adaptation and resistance). The application of RNA-seq was studied to investigate the genetic contribution of two Italian breeds (Maremmana and Chianina) with a different history of selection/divergence in meat production and quality by Bongiorni et al. [20] and Bongiorni et al. [7].They uncovered several differentially expressed genes that encode for proteins belonging to a family of tripartite motif proteins involved in growth, cell differentiation, and apoptosis. Therefore, this study aimed to investigate the effects of crossbreeding by the Montbeliard cattle on the transcriptome of the Sistani cattle.

Sample collection and RNA extraction
This experiment was approved and conducted under the supervision of the Scienti c-Technical Committee of Animal Science Research Institute (Karaj, Iran) with the ethic code of 2-13-13-010-960191. Ninety healthy pure and crossbred cows (pure Sistani breed (n = 45) and Sistani×Montbeliarde crossbreed (n = 45)) were selected for a sample collection from the cow farm in Zabol, Sistan province, Iran. All experimental animals, which belonged to a national crossbreeding program, were reared in the same location (Zabol city, Sistan province, Iran), farm, food, physiological, and management conditions. RNA was extracted from the blood samples by RNeasy (Qiagen) method and according to the manufacturer's protocols in the Iranian Media Gene Medicine Laboratory.

Library Preparation and Sequencing
According to the manufacturer, the mRNA selection, library preparation, and sequencing were performed at BGI Shenzhen, China on an Illumina HiSeqTM 2500 sequencer and 150bp paired-end reads speci cations. The next step, Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System were used in the quanti cation and quali cation of the sample library during the QC steps.

Quality Control and Alignment
Quality control was done using FastQC (version 1.6; http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc). Adapter sequences were removed from reads

Detection of differentially expressed genes
Following the assembling of the transcripts, differential expression analysis of two conditions (pure Sistani and crossbreed Sistani×Montbeliard) was performed using Cuffdiff (version 2.2.1.0). The Cuffdiff 2 method [26], implemented in Cu inks software [26], was used with default settings. Differentially expressed genes were considered at a FDR (false discovery rate corrected P-values) P < 0.05. We used CummeRbund (https://bioconductor.org/packages/devel/bioc/html/cummeRbund.html) for visualizing the DEGs.
Gene ontology and KEGG pathway enrichment analysis of differentially expressed genes As an international standardized gene functional classi cation system, GO offers a dynamic-updated controlled vocabulary and a strictly de ned concept to describe the properties of genes and their products comprehensively. To identify the signi cantly enriched GO terms, we mapped all DEGs to GO terms in the GO database (http://amp.pharm.mssm.edu/Enrichr/). In this process, the FDR (P < 0.05) was used as the threshold. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database is a resource for understanding high-level functions and utilities of biological systems, such as cells, organisms, and ecosystems from a molecular level. Also, it can be used to analyse large-scale Molecular datasets generated by genome sequencing and other high throughput experimental Technologies. KEGG pathway analyses were performed with the KEGG (https://www.genome.jp/kegg/kaas/).

RNA sequencing and assemble transcripts
The quality control performed using FastQC showed that the samples were of good quality based on the Phred scores (probability of incorrect base call). The analysis were carried out without trimming the reads because of the acceptable quality of the cDNA libraries and the fact that mapping algorithms use only reads that perfectly match a position on the reference genome. In total, we acquired 60575588-75529030 paired-end reads per sample. Based on Tophate2 results, the observed percentages of mapped reads to the bovine reference genome per sample were around 70.5-80.6%, and the percentages of unique aligned reads per sample were around 72.9-78.1%. The summary of reads mapping is provided in Table S1.

Differentially expressed genes between purebred Sistani and Sistani×Montbeliard crossbred
The differentially expressed genes identi ed by Cuffdiff 2 are shown in Table S2. Based on Cuffdiff 2 output, 45 gene were detected as signi cantly differentially expressed genes (false discovery rate q < 0.05) between pure Sistani breed and Sistani×Montbeliard crossbreed (Fig. 1).

Gene Ontology enrichment and pathway analysis
The GO enrichment analysis of the DEGs revealed 15 enriched pathways between pure Sistani breed and Sistani×Montbeliard crossbreed (Table 3). Based on the result of GO enrichment, revealed only the in ammation mediated by chemokine and cytokine signalling pathway signi cant GO terms (FDR q < 0.05). Cytokines and chemokines play a key role in the immune response to a myriad of pathogens in species (Fig. 2).

Discussion
Crossing the native breeds with the high-producing exotic cattle improves the production level of crossbreeds. The number of crossbred cattle has been increasing rapidly during the last ten years in Iran.
The population of native cattle in Iran includes Sarabi, Golpayegani, Sistani, Dashtiari, Najdi and Taleshi belonging to B. indicus [17]. Although the indigenous populations show good disease resistance and are adapted to harsh environmental conditions, their milk yields are low. On the other hand, breed selection will take many generations to improve production to a desirable level [24]. Therefore, one way to achieve a signi cant increase in production is the implementation of crossbreeding with the high-performing breed to combine the desirable characteristics of both breeds.In this regard, we compared transcripts of pure Sistani breed and Sistani×Montbeliard crossbreed populations to understand the effects of crossbreeding on gene expression.
In this study, two genes (i.e., PPPIR1B and ARID4A) were highly expressed in both populations (Heatmap correlation graph of expressed transcripts (Fig. 3). These genes are involved in the immune pathway and cell proliferation. The PPPIR1B (protein phosphatase 1 regulatory inhibitor subunit 1B) is involved in the cAMP signalling pathway. DARPP-32 and t-Darpp are proteins encoded by PPP1R1B [9] that they are overexpressed in several cancers, including those derived from gastric, colon, breast, and lung tissues [2,6]. The protein encoded by the ARID4A (AT-rich interaction domain 4A) is a ubiquitously expressed nuclear protein. This protein binds directly with several other proteins to retinoblastoma protein (pRB), regulating cell proliferation.
In this study, GO enrichment analysis revealed that the differentially expressed genes are enriched in many GO terms (Table S3). This study focuses on in ammation mediated by chemokine and cytokine signalling pathway, which is observed as the highest enriched pathway compared to the other pathways in this study. Chemokines and cytokines are two immune-modulating agent and play an essential role in controlling and regulating an immune response [13,21]. The chemokine system has shown a key role in both homeostasis, such as lymphoid organogenesis and leukocyte maturation and disease mechanisms [1,4,28].
We observed that PF4, vWF, CXCR4, CCL5, and RGS1 genes are deferentially expressed in the crossbreed Sistani and Montbeliard populations. These genes are involved in the in ammation mediated by chemokine and cytokine signalling pathways and are involved in resistance to disease and inappropriate environmental conditions. The Platelet factor 4(PF4), a small cytokine belonging to the CXC chemokine family is an abundant protein stored in α-granules of megakaryocytes (MK) and platelets. PF4 plays an important role in hemostasis and this factor acts as a lineage-speci c marker of megakaryocytic differentiation [12,18]. PF4 as a signi cant player in the initiation and development of in ammatory diseases [15,27].
Von Willebrand factor (vWF), is a glycoprotein that plays an important role in platelet adhesion and thrombosis formation. This factor speci cally synthesizes endothelial cells and platelet precursors. The species-speci c properties by functional characterization of bovine vWF gene promoter in the bovine endothelial cells have been demonstrated in a study by Janel et al. [16]. They reported that the sequences, which regulates the cell-speci c transcription, are not conserved in the bovine counterpart. It was also determined that the cis-acting elements do not behave identically in the bovine vWF promoter.
CXCR4 (c-x-c chemokine receptor type 4) is a molecule endowed with potent chemotactic activity for lymphocytes that this gene is involved in some developmental processes, including haematopoiesis, angiogenesis, and organogenesis [8,31]. This gene is highly expressed in the thymus, particularly in immature CD4 + and CD8 + cells [31]. Ashley et al. [5] reported that maternal CXCR4 was up regulated in sheep during early pregnancy due to implantation and placentation. Expression of C-X-C chemokine receptor type 4 (CXCR4) was reported in association with bovine viral diarrhea virus (BVDV) infection [25,29].
CC motif chemokine ligand 5 (CCL5) is an immune-associated gene belonging to the CC chemokine family [3]. This gene chemotactic for monocytes, T cells, and eosinophils. It also activates eosinophils and the release of histamine from basophils. Splice variants, transcript expression, and splicingassociated SNPs of the CCL5 gene in healthy and mastitis mammary glands in Holstein cows were studied by Yang et al. [30]. Based on their results, it was determined the chemokine CCL5 was down regulated in the mammary glands in the late stages of natural mastitis infection.
Regulator of G-protein signalling 1 (RGS1) is a member of the regulator of G-protein signaling family that accelerates Gai GTPase activity and acts to down-regulate the response to sustained chemokine activation [11,22]. Rgs1 is involved in chemokine receptor signalling in both T and B cells [10,23]. In particular, deletion of Rgs1 modi es the migratory behaviour of B cells within follicular areas and to increase germinal centre formation [10,14,22,23]. Also, the role of RGS1 in the control of lymphocyte homeostasis have been identi ed in other studies [10,22].
The consequences of cattle crossbreeding programmed on gene expression, signalling, and metabolic pathways, molecular networks, and biological functions are still indistinctive.

Conclusion
In conclusion, the presented study provided important information on the effect of crossbreeding by the Montbeliard cattle on the transcriptome of the Sistani cattle. Comparing crossbred vs. purebred, integrated analysis of differential gene expression and GO enrichment analysis allowed identifying disease-associated pathways such as infection of Mammalia. These genes are involved in pathways leading to immune and transcriptional responses, and they are associated with calving, fertility, and resistance to mastitis. Moreover, these pathways may have a role in different heat tolerance levels and disease resistance in the Sistani×Montbeliard crossbreeds.  Figure 1 Differential expression of genes between purebred Sistani and Sistani×Montbeliard crossbred: Volcano plots (the red point denotes the signi cantly differential expression and the points in the same region show the same expression)

Figure 2
In ammation mediated by chemokine and cytokine signalling pathway: The differentially expressed genes identi ed in the analysed cow blood and their involvement in the in ammation mediated by chemokine and cytokine signalling pathway