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 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-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-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-13], glutathione peroxidase 5 (Gpx5) [14, 15], murine epididymal retinoic acid-binding 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.

RESULT

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

The top 30 DEGs with at least a -8.22 log2-fold change difference are listed in Table , of the top 30 most DEGs, twenty-eight 28 were significantly downregulated in cattleyak with respect to yak. Whereas, the difference in the expression level of LY6G5C (lymphocyte antigen 6 family

member G5C) was the largest (-13.87 log2-fold change). Among the top 30 DEGs, majority belonged to ß - defensins members of gene family i.e DEFB129 (defensin beta 129), DEFB124 (Defensin. beta 124), DEFB122A (beta-defensin 122a), DEFB122 (beta-defensin 122) and CST11(cystatin-11 precursor), CLEC4D (C-type lectin domain family 4 member D), LY6G5C (lymphocyte antigen 6 family member G5C), LY6G5B (lymphocyte antigen 6 family member G5B), CLEC4D (C-type lectin domain family 4 member D) and UBQLN2 (ubiquilin 2) were all associated with antimicrobial activity, anti-mycobacterial immunity, sperm protection, motility and auto-immune responses [30-34].

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. ). 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. ). 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 production and conversion and replication recombination and repair (Fig. ).

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 up-regulated 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 down-regulated 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-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-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 Huntington's, Alzheimer's and Parkinson's disease [60-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. and ) and the oxidative phosphorylation was the top fourth enriched pathway in KEGG database analysis (Fig. ).

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-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 ER-associated proteins, the entire ER protein processing pathway might have been disrupted in male cattleyak (Fig. ).

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 co-chaperones 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 (XTP3-transactivated 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.

MATERIALS AND METHODS

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 paired-end 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:

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 t-test 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.

Table 1

Summary of read numbers aligned onto the bovine reference genome in the study. The number in brackets indicates the percentage of total reads aligned onto the bovine reference genome and/or matched at either unique or multiple genomic locations.

Samples	Total Reads	Mapped Reads	Uniquely Mapped Reads	Multiple Mapped Reads
YAK-M1	130422090	118831604 (91.11%)	114139684 (87.52%)	4691920 (3.60%)
YAK-M2	145932742	132900739 (91.07%)	128501066 (88.05%)	4399673 (3.01%)
YAK-M3	113814354	101028672 (88.77%)	97706628 (85.85%)	3322044 (2.92%)
CATTLEYAK-P2	126917080	115587937 (91.07%)	112254260 (88.45%)	3333677 (2.63%)
CATTLEYAK-P3	145833598	133053036 (91.24%)	129007683 (88.46%)	4045353 (2.77%)

Table 2

Validation of mRNA-Seq data of a selected set of yak and cattleyak genes with quantitative real-time PCR (qPCR) to confirm gene expression changes at the contact stage (treatment vs. control). *P<0.05. FC, fold change (treatment vs. control).

Gene ID	Symbol	Annotation	Log2FC
			qPCR	mRNA-seq
DEGs
ENSBTAG00000040244	APOL3	apolipoprotein L3 [Bos taurus]	6.21		5.83
ENSBTAG00000001801	SYT4	synaptotagmin-4 [Bos taurus]	14.98		5.71
ENSBTAG00000006174	TMEM212	PREDICTED: transmembrane protein 212 isoform X1 [Bos taurus]	5.04		7.46
ENSBTAG00000040268	UBE2U	PREDICTED: ubiquitin-conjugating enzyme E2 U [Bison bison bison]	12.91		4.56
ENSBTAG00000011164	CES5A	carboxylesterase 5A precursor [Bos taurus]	0.24		-10.96
ENSBTAG00000000601	COL11A2	TPA: collagen alpha-2(XI) chain precursor [Bos taurus]	6.09		4.95
ENSBTAG00000013415	CST11	cystatin-11 precursor [Bos taurus]	0.15		-10.86
ENSBTAG00000019600	Ooep	oocyte-expressed protein homolog [Bos taurus]	7.57		6.88
ENSBTAG00000005628	CD52	CAMPATH-1 antigen precursor [Bos taurus]	0.03		-9.15
ENSBTAG00000031254	DEFB124	beta-defensin 124 precursor [Bos taurus]	0.13		-12.48
ENSBTAG00000038972	LCN9	Epididymal-specific lipocalin-9, partial [Bos mutus]	0.76		-13.33
ENSBTAG00000031772	PRODH2	TPA: probable proline dehydrogenase 2 [Bos taurus]	0.27		6.83
ENSBTAG00000003540	MYT1L	myelin transcription factor 1-like protein [Bos taurus]	16.17		4.44

Table 3

Number of functional annotations of the differentially expressed genes (DEGs) of yak and cattleyak respectively, in the study.

Annotated Databases	Annotated	COG	GO	KEGG	KOG	NR	Swiss-Prot
L01_L02_L03_vs_L04_L05	2905	570	2327	1267	1565	2904	2103

Table 4

List of the 30 most differentially expressed genes between epididymal samples of cattleyak (C) and yak (Y).

Gene Symbol	Ensemble ID	Annotation	Log2-FC	FDR
LY6G5C	ENSBTAG00000025449	lymphocyte antigen 6 complex locus protein G5c [Bubalus bubalis]	-13.87	0.006590395
LCN9	ENSBTAG00000038972	Epididymal-specific lipocalin-9, partial [Bos mutus]	-13.33	0.000100162
B124	ENSBTAG00000031254	beta-defensin 124 precursor [Bos taurus]	-12.48	0.000368901
DEFB127	Bos_taurus_newGene_201576	beta-defensin 127 [Bos taurus]	-11.99	7.1802E-10
SPINT4	ENSBTAG00000009626	Kunitz-type protease inhibitor 4, partial [Bos mutus]	-11.92	0.00065468
DEFB126	Bos_taurus_newGene_201580	beta-defensin 126 [Bos mutus]	-11.00	0.011069676
CES5A	ENSBTAG00000011164	carboxylesterase 5A precursor [Bos taurus]	-10.964	0.001093997
LOC103003268	Bos_taurus_newGene_203067	colostrum trypsin inhibitor-like isoform X1 [Bos taurus]	-10.94	0.003923811
Tctp	ENSBTAG00000048229	translationally-controlled tumor protein [Bos taurus]	-10.88	0.018211012
CST11	ENSBTAG00000013415	cystatin-11 precursor [Bos taurus]	-10.86	0.018686585
LOC102282407	Bos_taurus_newGene_185995	uncharacterized protein [Bos mutus]	-10.342	0.0084299
TEDDM1	ENSBTAG00000021432	transmembrane epididymal protein 1 [Bos taurus]	-10.242	0.003437677
M91_02466	ENSBTAG00000023301	hypothetical protein [Bos mutus]	-10.15	0.000925333
#N/A	Bos_taurus_newGene_85005	TPA: ribonuclease A K2 [Bos taurus]	-9.40	0.010653548
DEFB129	ENSBTAG00000048288	beta-defensin 129 [Bos taurus]	-9.40	0.003773411
CD52	ENSBTAG00000005628	CAMPATH-1 antigen precursor [Bos taurus]	-9.15	0.014821799
CRYBA4	ENSBTAG00000019282	beta-crystallin A4 [Ovis aries]	-9.12	0.000509967
ATF4	ENSBTAG00000017462	TPA: cyclic AMP-dependent transcription factor ATF-4 [Bos taurus]	-9.07	6.04979E-23
LY6G5B	ENSBTAG00000039740	lymphocyte antigen 6 complex locus protein G5b precursor [Bos taurus]	-9.05	0.000955208
DEFB128	Bos_taurus_newGene_201574	Beta-defensin 128 [Bos mutus]	-9.02	1.58001E-33
TKTL1	ENSBTAG00000020127	transketolase-like protein 1 [Bos taurus]	-8.85	0.003542531
#N/A	Bos_taurus_newGene_20024	reverse transcriptase-like [Bos taurus]	8.69	1.06E-10
PHALS_10256	Bos_taurus_newGene_203836	hypothetical protein M91_10256, partial [Bos mutus]	8.64	0.031058549
LOC111701397	ENSBTAG00000023300	sperm acrosomal protein FSA-ACR.1-like [Bos taurus]	-8.48	0.002642282
LOC107133408	Bos_taurus_newGene_42683	endogenous retrovirus group K member 18 Env polyprotein-like [Bos taurus]	-8.47	0.0212064
LOC111095245	ENSBTAG00000019584	epididymis-specific alpha-mannosidase-like isoform X3 [Bos taurus]	-8.36	0.005175811
TIMP	ENSBTAG00000012131	Metalloproteinase inhibitor 1, partial [Bos mutus]	-8.30	0.009554343
GDF3	ENSBTAG00000012171	Growth differentiation factor 3 [Bos taurus]	-8.30	0.025215041
CLEC4D	ENSBTAG00000039347	C-type lectin domain family 4 member D [Bos taurus]	-8.26	2.43E-11
C2H2orf80	ENSBTAG00000012888	chromosome 2 C2orf80 homolog	-8.22	0.017896104