Comparative RNA-Seq Analysis of Differentially Expressed Genes in the Epididymides of Yak and Cattleyak

Background: Cattleyak are the Fl hybrids between (♀) yak (Bos grunniens) and (♂) cattle (Bos taurus). Cattleyak exhibit higher capability in adaptability to a harsh environment and display much higher performances in production than the yak and cattle. The cattleyak, however, are females fertile but males sterile. All previous studies greatly focused on testes tissues to study the mechanism of male infer-tility in cattleyak. However, so far, no transcriptomic study has been conducted on the epididymides of yak and cattleyak. Objective: Our objective was to perform comparative transcriptome analysis between the epididymides of yak and cattleyak and predict the etiology of male infertility in cattleyak. Methods: We performed comparative transcriptome profiles analysis by mRNA sequencing in the epidi-dymides of yak and cattleyak. Results: In total 3008 differentially expressed genes (DEGs) were identified in cattleyak, out of which 1645 DEGs were up-regulated and 1363 DEGs were down-regulated. Thirteen DEGs were validated by quantitative real-time PCR. DEGs included certain genes that were associated with spermatozoal matura-tion, motility, male fertility, water and ion channels, and beta-defensins. LCN9, SPINT4, CES5A, CD52, CST11, SERPINA1, CTSK, FABP4, CCR5, GRIA2, ENTPD3, LOC523530 and DEFB129, DEFB128, DEFB127, DEFB126, DEFB124, DEFB122A, DEFB122, DEFB119 were all downregu-lated, whereas NRIP1 and TMEM212 among top 30 DEGs were upregulated. Furthermore, protein processing in endoplasmic reticulum pathway was ranked at top-listed three significantly enriched KEGG pathways that as a consequence of abnormal expression of ER-associated genes in the entire ER protein processing pathway might have been disrupted in male cattleyak which resulted in the down-regulation of several important genes. All the DEGs enriched in this pathway were downregulated ex-cept NEF. Conclusion: Taken together, our findings revealed that there were marked differences in the epididymal transcriptomic profiles of yak and cattleyak. The DEGs were involved in spermatozoal maturation, mo-tility, male fertility, water and ion channels, and beta-defensins. Abnormal expression of ER-associated genes in the entire ER protein processing pathway may have disrupted protein processing pathway in male cattleyak resulting in the downregulation of several important genes involved in sperm maturation, motility and defense.


INTRODUCTION
Cattleyak (CY) are the Fl hybrids between yak (YK) (♀) and cattle (CL) (♂). The cattleyak, also known as "khainag", are females fertile but males sterile. The backcrossing of F1 female cattleyak with either male yak or cattle respectively produces B1 backcross animal called "ortoom", that is also male sterile but female fertile. Moreover, backcrossing of B1 animal females either male yak or cattle produces B2 backcross animal called the "usanguzees", which are both male and female fertile [1]. Cattleyak display higher potential in *Address correspondence to this author at the School of Life Science and Engineering, Southwest University of Science and Technology, Mianyang 621010, Sichuan, China; Tel: +86 08162419792; E-mail: caixin2323@126.com adapting to harsh environment and exhibit much higher performances in production such as economic traits, milk and meat production as compared to both yak and cattle. Cattleyak, therefore, has a great impact on the socio-economic development in the yak breeders in plateau regions, although how to make F1 male cattleyak fertile to utilize its hybrids is the key task in yak breeding [2].
Following spermatogenesis, spermatozoa formed in the testis are carried to the epididymis, where post-testicular maturation of the sperm cells occurs. Post-testicular maturation of spermatozoa mainly involves motility acquisition and the competence to undergo capacitation, leading to the ability to fertilize an egg [3,4]. The epididymis is functionally and structurally divided into three distinctive anatomical segments: the caput, the corpus and the cauda. Each segment possesses a distinct pattern of gene expression associated with physiological functions needed for the sperm maturation and motility [5][6][7]. The epididymis performs three major functions. The first main function is the transportation of sperm, which is accomplished by the smooth muscle cells contractions around the epididymal tube, and by constant movement and fluid secretions from the testis. The second main function is the maturation of spermatozoa. In the course of their passage through the epididymis, spermatozoa become fully mature and motile and get the ability to fertilize the oocyte. Maturation goes through a series of changes such as nuclear condensation and modifications in the proteins and lipids composition of the sperm plasma membrane. The third important function is the storage of spermatozoa, especially in the caudal epididymal segment [8,9]. During the maturation of spermatozoa, they undergo a series of molecular and biochemical changes as they pass through the epididymis. These changes are believed to be initiated by the transport of ions and the addition of proteins secreted in the epididymal duct lumen [8][9][10]. Genes expressed in the epididymis have a very important role in sperm maturation and motility; alterations in the expression levels of these epididymal specific genes may lead to infertility.
In recent studies, the certain number of genes have been identified to have putative roles in sperm maturation and motility. Examples of such genes are ROS proto-oncogene 1, receptor tyrosine kinase (Ros1) [11][12][13], glutathione peroxidase 5 (Gpx5) [14,15], murine epididymal retinoic acidbinding protein (Erabp), Acidic epididymal glycoprotein (Aeg1) [16,17], human epididymal proteins types 1-6 [18,19] and Bin1b [20], while many genes yet remain unknown [21,22]. Male Cattleyak infertility is one of the major hindrances in the yak breeding. Over the last few decades, researches have been conducted on the testes of YK and CY in order to investigate the mechanisms of infertility in male CY. No transcriptomic study, however, has been conducted on the epididymides of yak and cattleyak so far. In order to predict the etiology of infertility and to determine the differential gene expression patterns in the epididymal transcriptome of the YK and CY, we performed comparative mRNA sequencing (RNA-seq) analysis between the epididymal tissues of YK and CY. This study will enhance our understanding of the molecular mechanisms and differences in genes expression profiles between yak and cattleyak epididymal transcriptomes.

Sequencing and Assembly of Yak and Cattlyak Transcriptome
The cDNA libraries were generated using the mRNA from six epididymal tissue samples of yak and cattleyak (i.e., M1, M2, M3 P1, P2 and P3). Six biological samples of epididymis, i.e., a total of five samples were completed, three from the yak (control) and two replicates from cattleyak (treatment), were separately subjected to RNA-sequencing analysis. A total of 99.28 Gb of clean data were obtained from the five epididymal samples altogether, with Clean Data of 0.07Gb for each sample and 0.01% for Q30 ≥ bases and above. The efficiency of aligning reads with specified reference genome UMD3.1 of each sample ranged from 88.77% to 91.24% and matched to either unique or multiple genomic locations ( Table 1). In total, 109,496 unigenes including 21,681 new genes were mined in the bovine transcriptome. According to the databases sequence analysis of the newly discovered genes with Non-redundant protein sequences, NR [23], Swiss-Prot [24], Gene Ontology GO [25], Orthologous Groups COG [26] and Kyoto Encyclopedia of Genes and Genomes KEGG [27] databases using BLAST [28] total 18576 new genes were annotated in these databases. The replicates of yak and cattleyak repeated well with each other (r2=0.98~1.00) [29].

Comparative Analysis of Gene Expression
Comparative analysis of gene expression was performed for the yak and cattleyak epididymal tissues. In total 3008 unigenes, differential expressed genes were screened (Fig.  1). Three thousand and eight unigenes/DEGs, including 1363 down-regulated and 1645 up-regulated ones, were differentially expressed with false discovery rate (FDR) <0.05 and fold change (FC) ≥2 (Fig. 2). For the validation of RNA-seq data, we randomly selected 13 DEGs (APOL3, SYT4,  TMEM212, UBE2U, CES5A, COL11A2, CST11, Ooep,  CD52, DEFB124, LCN9, PRODH2 and MYT1L). The expression patterns of all thirteen differentially expressed genes (DEGs) were validated by quantitative real-time PCR (qPCR) and the results were consistent with those of the mRNA-Seq analysis ( Table 2). These results demonstrated that there is a significant difference in the number of differentially expressed genes in yak epididymal tissue compared to cattleyak epididymal tissue. Results of the functional   annotation indicated that all the DEGs had significant matches in the Nr database, and some of them also had annotation information in the Swiss-Prot, GO, KEGG and COG databases ( Table 3).
Based on the functional annotation obtained from GO enrichment analysis, the yak and cattleyak DEGs were separated into 54 functional groups, which belong to three main categories: biological processes (60 DEGs), cellular components (54 DEGs), and molecular functions (68 DEGs) (Fig.  3). In the biological process category, greater percentages of DEGs were involved in cell aggregation, cell killing, biological phase, growth, locomotion, localization, reproductive process, developmental process, multicellular organismal process, response to stimulus, multi-organism process and biogenesis compared to all unigenes of yak and cattleyak epididymal tissue. More proportions of DEGs in the cellular component category were localized to the collagen trimer, extracellular matrix part, synapse part, organelle part, synapse, membrane, membrane part and organelle. The DEGs in the molecular function category were more enriched in GO class of electron carrier activity, chemoattractant activity, antioxidant activity, transporter activity, enzyme regulator activity, receptor regulatory activity, channel regulator activity and nucleic acid binding transcription factor activity than compare to all unigenes.
Analysis of KEGG pathways was performed to determine the biological functions of the DEGs. Two hundred and sixty-nine (269) DEGs were allocated to two hundred and seventy-one (271) KEGG pathways (Fig. 4). The taste transduction pathway accounted for the highest number of DEGs (58 DEGs) followed by Huntington's disease, Protein processing in endoplasmic reticulum, Oxidative phosphorylation and Alzheimer's disease. In total, forty DEGs were allocated in protein processing in the endoplasmic reticulum pathway. ko04141 (Protein processing in the endoplasmic reticulum). All the DEGs enriched in ER protein processing pathway were downregulated except NEF which was upregulated in CY compare to YK. COG annotation grouped the DEGs into twenty-three COG functional classes. The top five COG classes with the highest number of DEGs included those involved in general function prediction only, posttranslational modification, protein turnover, chaperones, translation, ribosomal structure and biogenesis, signal transduction mechanisms, energy pro-duction and conversion and replication recombination and repair (Fig. 5).

DISCUSSION
The comparative transcriptome profiles analysis provided an outline of the transcriptomic diversity in yak and   cattleyak epididymal DEGs. Cattleyak exhibit higher capability in adaptability to a harsh environment and display much higher performances in production than the yak and cattle. Previous studies conducted on male cattleyak infertility have focused on testes whereas, so far, none of such studies have been conducted on the epididymal tissues in the male cattleyak to investigate the mechanism of infertility [35]. In this study, we performed comparative transcriptome analysis in the epididymides of yak and cattleyak to predict the etiology of male cattleyak infertility.
Given the fact that a number of genes could be differentially expressed in the epididymis of yak and cattleyak, transcriptomic comparison between the epididymes of yak and cattleyak is an effective approach to study or discover the key genes involved in different metabolic pathways that are specifically involved in sperm maturation and fertility.
Comparative analysis of gene expression between yak and cattleyak epididymal tissues was performed resulting in 3008 DEGs; 1363 down-regulated DEGs and 1645 upregulated DEGs. Our results demonstrated a significant difference in the number of DEGs between cattleyak and yak epididymal tissue. Functional annotation analysis performed for DEGs indicated that all the DEGs had significant matches in the Nr database while some DEGs additionally had annotation information in the Swiss-Prot, GO, KEGG and COG databases.
Interestingly, our results show that group of genes downregulated in the epididymis of infertile male cattleyak compared with yak encodes for proteins that are well known to play pivotal role in the host defense and spermatozoal maturation such as, DEFB129, DEFB128, DEFB127, DEFB126, DEFB124, DEFB122A, DEFB122, DEFB119 [36], CRISP1 [37] and SPINLW1 [38]. Since the ß -defensins were first identified in the epithelial cells of cattle airway as antimicrobials, more than 30 β-defensins genes have been described in bovines, human and other species [39][40][41]. Several are constantly expressed in the reproductive tract of male adult tissues and are particularly found to be abundant in the testis and epididymis [39, 42,43]. In some specific cases, their expression pattern is almost very abundant or restricted in the adult epididymal epithelium, where the majority of them depend on the β-defensin evaluated exhibit unique cell and region-specific expression patterns. It is also believed that Beta-defensins serve a very important role in both fertility and host defense [44][45][46]. Such as in innate immunity function of sperm coating proteins, inhibition of proteases may directly attack the sperm plasma membrane and antimicrobial activity. CRISP1 has been shown to play an important function in egg-sperm interaction and inhibition of the uptake of the ions, for instance, Ca2þ, is needed for spermatozoal capacitation [47]. Additionally, SPINLW1 binds with semeno-gelin in seminal plasma and on spermatozoa after ejaculation. It is believed that this protein complex is a part of a larger protein network on the sperm surface that acts as a protective shield in the female reproductive tract before capacitation [48].
Furthermore, another important group of genes, Lipocalins such as LCN9 (lipocalin 9), LCN6 (lipocalin 6), LCN12 (lipocalin 12) [49], which are all downregulated in the epididymis of infertile male cattleyak as compared to yak, consist of a family of proteins that play a role in the transportation of hydrophobic ligands, for instance, fatty acids retinoids and steroids to both the epithelium lining of the epididymis and the spermatozoa. LCN9 (-12.80 log2-fold change) which is reported to function in the transportation of luminal hydrophobic protein ligands to specific cells of epithelium lining of the epididymis [31] has an important but non-redundant function in male fertility. The functions of all the other epididymal Lipocalins, however, have not been so far well understood [50]. Interestingly, other genes having proteolytic activity, such as CTSK and ADAMDEC1, are also downregulated in epididymis of infertile male cattleyak compared with yak. CTSK is proteases that function in the degradation and processing of proteins and sperm maturation [51], while ADAMDEC1 is a proteolytically active metzincin metalloprotease [51]. Conversely, SPINT4 is an epididymis-specific protein having anti-proteolytic activity [52].
The remaining top 30 most DEGs were associated with CES5A (carboxylesterase 5A), an epididymis-specific gene which is crucial in sperm maturation and male fertility [53]. CD52 (CD52 molecule) is mainly secreted in the lumen of the epididymis and via the glycosyl-phosphatidylinositol (GPI) anchor. It is introduced into the sperm membrane during the passage of the spermatozoa through the epididymis. This is necessary for the maturation and mobility of sperm. Downregulation of CD52 is strongly associated with strong immobility and also lack of ability of sperm to penetrate the zona pellucida [54,55]. SERPINA1 (serpin family A member 1), lack of protein C inhibitor Serpina5 were found to cause infertility [56] while LOC523530 (nuclear RNA export factor 3), an epididymis-specific expression in principle cells regulates the expression of TGF-β3 in an mRNA export activity [57]. FABP4 (fatty acid-binding protein 4) associated with the stabilization of the sperm flagellum [58], CCR5 (C-C chemokine receptor type 5) essential for sperm motility, hyper-activation and acrosome reaction [18], GRIA2 (glutamate ionotropic receptor AMPA type subunit 2) functions in ion channels [31] and ENTPD3 (Ectonucleoside triphosphate) participates in smooth muscle contraction, steroidogenesis and spermatogenesis, and also regulates trans epithelial transport [59].
In the KEGG enrichment analysis of the DEGs, we selected the top four pathways obtained from the KEGG database analysis that were the taste transduction, Huntington's disease, protein processing in the endoplasmic reticulum, and the oxidative phosphorylation. Huntington's disease was identified as the top second in the KEGG database analysis. In the previous studies, it has been reported that ER stress is associated with many conformational diseases, such as Hun-tington's, Alzheimer's and Parkinson's disease [60][61][62], however, in our findings endoplasmic reticulum protein processing pathway was ranked at top three position in the KEGG database analysis, in which most of the enriched DEGs were downregulated (Figs. 4 and 5) and the oxidative phosphorylation was the top fourth enriched pathway in KEGG database analysis (Fig. 4).
The endoplasmic reticulum (ER) is membrane-bounded cytoplasmic organelle in the eukaryotic cells which serves many general functions such as synthesis and folding of proteins, lipids, and transportation of these synthesized proteins and lipids in vesicles to the Golgi apparatus. ER stress can be briefly defined as a disproportion between the protein folding ability of the ER and the protein load, consequential leads to misfolded proteins accumulation [63]. In these studies, it has been speculated that ER stress is involved in certain conformational diseases, such as Huntington's, Alzheimer's and Parkinson's disease [60][61][62]. Our findings showed that all the DEGs enriched in the "protein processing in endoplasmic reticulum pathway" process were downregulated except NEF in male cattleyak as compared to yak. We speculate that as a consequence of abnormal expression of ERassociated proteins, the entire ER protein processing pathway might have been disrupted in male cattleyak (Fig. 5).
ER chaperones comprise of molecular chaperones and folding enzymes localized in the ER, the function of which is to assist in the unfolding or folding and assembly of nascent proteins [64]. BIP (Binding protein) is a member of the HSP70 family, which is an important molecular chaperone by binding to the hydrophobic region of unfolded proteins and helps in folding [65]. Moreover, SEC63 and DNAJC3 are also molecular chaperones that have been reported as cochaperones which can regulate the function of BIP, downregulated in male cattleyak epididymal tissues. Furthermore, in our findings, numerous protein folding enzymes, such as OST (dolichyl-diphosphooligosaccharide protein glycosyltransferase) and PDI (protein disulfide isomerases) were reduced in the male cattleyak. These protein folding enzymes help in the oxidization of cysteine residues of nascent proteins and coordinate to make disulfide bonds. However, ERO1 (ER oxidoreduction) is involved in the reoxidation of reduced folding enzymes which was also down-regulated [66]. Additionally, the downregulation of XTP3B (XTP3transactivated gene B) and OS9 (osteosarcoma 9) which are very important ERAD (ER-associated degradation) complex components responsible for the identification of unfolded proteins, was downregulated in male cattleyak compare to yak [67-69]. Together, our findings suggest that downregulation of most of the DEGs involved in protein processing in the ER pathway may cooperate to cause the accumulation of some misfolded proteins in the ER, that leads to deficiency of vital proteins such as beta-defensins Lipocalins, etc. which are important in the sperm maturation, motility and host immunity, that might be the cause of infertility in male cattleyak. Additionally, the subsets of DEGS from this and other altered pathways could be further investigated first for technical validation and then for clinical applications.

Animals and Epididymal Sample Collection
Sample collection was carried out under license in accordance with the Guidelines for Care and Use of Laboratory Animals of China and all protocols were approved by the Institutional Review Board of Southwest University of Science and Technology [70]. Male cattleyak (Holstein Maiwa yak; n ¼ 3, named P1, P2, and P3) and male yak (Maiwa yak; n ¼ 3, named M1, M2, and M3) aged 12 months were sampled from a Maiwa yak population fed on a pasture in Hongyuan county, Sichuan province of China, in which cattleyak were F1 generations of Holstein and Maiwa yak. The epididymis of each animal was obtained by veterinary surgical operation and fat and connective tissues were removed. Epididymes were separated apart from testis by fine-scale dissection and were snap-frozen in liquid nitrogen (-196°C), transported to the laboratory and stored at -80°C until RNA isolation.

RNA Extraction, Library Preparation and Sequencing
Three biological replicates were used for yak and two biological replicates for cattleyak as in the previous studies [71,72]. Total RNA extraction was performed from each epididymis frozen samples of yak and cattleyak. RNA was extracted using Trizol (Life Technologies, LT) as per the manufacturer's protocol [73,74]. RNA degradation and DNA contamination was monitored on 1.5% agarose gel. The RNA concentration and purity were confirmed by using the NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) according to the manufacturer's instructions. The RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, USA).
A total amount of 1.5 µg RNA per sample was used as an input material for rRNA removal using the Ribo-Zero rRNA Removal Kit (Epicenter, Madison, WI, USA). Libraries sequences were generated using NEBNext R Ultra ™ Directional RNA Library Prep Kit for Illumina R (NEB, USA) following the manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primer and Reverse Transcriptase. Second-strand cDNA synthesis was succeeded using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities after adenylating of 3' ends of DNA fragments. NEBNext Adapter with hairpin loop structures was ligated to prepare for hybridization. In order to insert fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP beads (Beckman Coulter, Beverly, USA). Then 3 µL of USER Enzyme (NEB, USA) was used with size-selected, adapter-ligated cDNA at 37°C for 15 min before PCR. PCR was then performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 and qPCR.
An AcBot Cluster Generation System was used for the clustering of the index-coded samples by using TruSeq PE Cluster Kitv3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the libraries were sequenced on an Illumina Hiseq Xten platform (Biomarker Technologies Co. LTD, Beijing, China) and pairedend reads were generated.

Data Processing, and Assembly of Yak and Cattleyak Transcriptome
Raw data (raw reads) of fast format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyzes were based on clean data with high quality. The specified reference genome used was cow UMD3.1, for the sequence alignment of the clean reads from yak and cattleyak and subsequent analysis. Cow (Bos taurus) UMD3.1 (GCA_000003055.3) reference genome was retrieved from ensemble database: http://www.ensembl.org/ Bos_taurus/Info/Index. The alignment of yak and cattleyak clean reads to reference genome UMD3.1 was performed using HISAT2 in an efficient alignment system. The transcriptome was assembled using the StringTie based on the reads mapped to the reference genome. StringTie (1.3.1) https://ccb.jhu.edu/software/stringtie/index.shtml was used to calculate FPKMs (fragments per kilobase of transcript per million fragments mapped), of transcripts abundance in each sample. The number of fragments extracted from one transcript is related to the amount of sequenced data, the length of the transcript, and the expression level of the transcript. In order for the number of fragments to truly reflect the expression level of the transcript, the number of fragments in the Mapped Reads Number and transcript length were normalized. StringTie uses Fragments per Kilobase of transcript per Million fragments mapped as a measure of transcript or gene expression level, and FPKM were calculated as follows: !"#$ = !"#$ !"#$%&'() !"##$% !"#$%&'() !"##"$%& ×!"#$%&"'(! !"#n!"ℎ(!") The gff compare program was used to compare the gene and transcripts with the annotation. Differential expression analysis of two groups of yak and cattleyak samples was performed using the DESeq R package (1.10.1). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value <0.01 and an absolute value of log2 (fold change) > 2 found by DESeq were assigned as differentially expressed. We finally chose Fold Change ≥ 2 and FDR < 0.05 as the screening criteria.

Functional Annotation and Functional Enrichment Analysis of Differentially Expressed Genes
For the functional annotation of the transcripts, the unigenes were blasted against the databases of Nr, Swiss-Prot, KEGG, and COG using BLAST program with an E-value ≤ 1e-5. Te Blast2GO program was used to annotate the major GO categories of genes with an E-value ≤ 1e-5 [75].

Pathway Enrichment Analysis
GO enrichment analysis Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the topGO R packages. For KEGG pathway enrichment analysis, we used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways [76].

qPCR Validation of mRNA-Seq
To validate the repeatability and reproducibility of gene expression data from RNA-sequencing, quantitative reverse transcriptase PCR (qRT-PCR) analysis of 10 DEGs was performed on the same samples. Based on their potential functional importance, 12 and 10 DEGs were selected for validation by qPCR from yak and cattleyak, respectively. The cDNA was prepared from the remaining total RNA after transcriptome sequencing according to the instructions of the SuperScript™ III Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA). A SYBR Green assay was used to quantify the expression of each gene using SYBR Premix Ex Taq (TaKaRa, Dalian, China) in a CFX Connect™ Real-Time PCR Detection System (Bio-Rad, Munich, Germany) and Actin as the internal control, with the following cycle conditions: 95°C for 5 min, followed by 45 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 20 s. Three yak and cattleyak biological replicates were analyzed for each gene with three technical replicates. Data were processed using the 2−ΔΔCt method [77] and analyzed statistically using the Student's ttest in Sigmaplot v.8 software to compare the difference between the yak and the cattleyak samples at P <0.05. The efficiency of aligning Reads with reference genome of each sample ranged from 88.77% to 91.24%.

CONCLUSION
In conclusion, we have shown that several families of genes were downregulated in the epididymides of infertile male cattleyak compared to yak. Many of these genes appear to be specifically involved in the regulation of spermatozoal maturation, motility, host defense and the regulation of ions in the epididymal lumen. This transcriptome dataset will enrich the transcriptomic information for male cattleyak, and provide fundamental support for future research on the molecular mechanisms governing male cattleyak infertility. In addition, a large number of DEGs were obtained, which can be useful for further functional analysis studies in cattleyak breeding research. Furthermore, our finding also revealed very important pathways that are suggested to be associated with infertility of male cattleyak. Especially, protein processing endoplasmic reticulum and oxidative phosphorylation are suggested to be further investigated for a better understanding of the molecular mechanism and etiology of infertility in male cattleyak.

HUMAN AND ANIMAL RIGHTS
Sample collection was carried out under licence in accordance with the Guidelines for Care and Use of Laboratory Animals of China.

CONSENT FOR PUBLICATION
Not applicable.

AVAILABILITY OF DATA AND MATERIALS
The data supporting the findings of the article is available in the Biomarker at www.biomarker.com.cn, with reference number BMK170616-F748-0101.

FUNDING
This work was funded by the National Natural Science Foundation of China (31601946) and Longshan Academic Talent Research Supporting Program of SWUST (17LZ X671 and 18lzx661).