Search tips
Search criteria 


Logo of oncotargetLink to Publisher's site
Oncotarget. 2017 August 15; 8(33): 54416–54433.
Published online 2017 April 13. doi:  10.18632/oncotarget.17081
PMCID: PMC5589591

Analysis of Bos taurus and Sus scrofa X and Y chromosome transcriptome highlights reproductive driver genes


The biology of sperm, its capability of fertilizing an egg and its role in sex ratio are the major biological questions in reproductive biology. To answer these question we integrated X and Y chromosome transcriptome across different species: Bos taurus and Sus scrofa and identified reproductive driver genes based on Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm. Our strategy resulted in 11007 and 10445 unique genes consisting of 9 and 11 reproductive modules in Bos taurus and Sus scrofa, respectively. The consensus module calculation yields an overall 167 overlapped genes which were mapped to 846 DEGs in Bos taurus to finally get a list of 67 dual feature genes. We develop gene co-expression network of selected 67 genes that consists of 58 nodes (27 down-regulated and 31 up-regulated genes) enriched to 66 GO biological process (BP) including 6 GO annotations related to reproduction and two KEGG pathways. Moreover, we searched significantly related TF (ISRE, AP1FJ, RP58, CREL) and miRNAs (bta-miR-181a, bta-miR-17-5p, bta-miR-146b, bta-miR-146a) which targeted the genes in co-expression network. In addition we performed genetic analysis including phylogenetic, functional domain identification, epigenetic modifications, mutation analysis of the most important reproductive driver genes PRM1, PPP2R2B and PAFAH1B1 and finally performed a protein docking analysis to visualize their therapeutic and gene expression regulation ability.

Keywords: spermatogenesis, spermiogenesis, fertilization, WGCNA, reproduction


The sex chromosomes XY or ZW in case of mammals or birds varies both phylogenetically and structurally among different taxonomical classes with many similar features and is under selective pressure that shapes their evolution [1, 2, 3]. The karyotypes shows that the hemizygous chromosome Y or W as short, heterochromatic containing few genes while on the other hand the homozygous chromosomes X or Z is of same size to that of autosomes and contains similar amount of genes although it tends to cluster functionally similar genes together most of which are enriched for reproduction and differentiation [4].

Reproduction is the fundamental biological process by which animals produce their offspring. In mammals it is done by sexual reproduction where a haploid sperm from male fertilizes a haploid oocyte to restore the chromosome number. To explore the reproductive driver genes that can improve the reproductive health and fertility is at the core of reproductive biology. The recent advances in sequencing technologies and its analysis have a strong potential to find out genes that can improve the overall reproductive performance of farm animals including cow and pig. To this end, we first integrated X and Y chromosome transcriptome across different species: Bos taurus and Sus scrofa to identify reproductive driver genes based on Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm. The schematic diagram for a multi-step strategy to identify reproductive driver genes is shown in Figure Figure11.

Figure 1
A multi-step Strategy

The paternal genetic material across different species is transferred to oocyte by sperm after several dramatic morphological and physical changes [5]. A global transcriptional silencing occurs after meiosis as sperm nuclear basic proteins replaces somatic histones and cause tight DNA compaction [6, 7], and also most of the cytoplasmic material including ribosomes are discarded [8]. At this stage post translational mediators during the process of morphogenesis and motility of sperm become most important as there is no global transcription and translation [5]. The sperm motility is crucial for its ability to fertilize oocyte where PP1 a sperm specific phosphatase in mouse [9, 10] and GSP 3/4 in C. elegans have been shown important phosphatases for sperm motility [5]. The process of male and female gametes i.e. sperm and oocyte to unite and form a zygote is termed as fertilization. The sperm passes several developmental stages before it can reach oocyte in the female body to fertilize it. One such process is when the spermatids start to develop tails which will help them to swim in female reproductive tract to reach oocyte. This process is under strict control of gene expression. It has been observed that a change in gene expression level preferentially benefits X or Y chromosome bearing spermatozoa to develop tails and causes a skewed female or male sex ratio. In this scenario, genes that are responsible for sperm maturity takes a central position as candidate genes for sex ratio control. The WGCNA analysis also determines a cohort of genes that determine the process of sperm maturity and their preferential selection for tail development of specific X or Y chromosome bearing spermatozoa hence providing grounds for their selection as candidate genes in sex ratio control [5].

The genes involved in sexual development, spermatogenesis, spermatid development and differentiation and pregnancy are of vital importance in reproductive biology. Spermatogenesis is a complex process that is involved in the development of highly differentiated haploid sperms from diploid spermatogonial stem cells that is further divided into three phases of almost equal length, the spermatogoniogenesis [11], meiosis that generates four haploid spermatids [12] and spermiogenesis a process of dramatic morphological changes from spermatid to mature spermatozoa [13]. Spermatogenesis sequence of happenings is well defined in mammals and takes place in seminiferous tubules of the testes that is initiated at puberty, a continuous process in reproductively active adults to maintain the continuous sperm production with a fully controlled timing, in mouse it is 35 days [14].

The present study explores the integrated X and Y chromosome trancriptome across Bos taurus and Sus scrofa and identified important reproductive driver genes including PRM1, PPP2R2B and PAFAH1B1. Interestingly various immune related genes are also expressed together with genes involved in reproduction suggesting its probable role in reproduction specifically in preferential selection of X and Y chromosome bearing spermatozoa before fertilization. The preferential selection of X or Y chromosome bearing spermatozoa by female can have a profound importance for farm animals where farmers can select animals gender without the risk of losing subsequent fertilization ability of sperm in female oviduct that are exposed to sorted sperm [15]. Furthermore, the present study conducts complete genetic analysis of the genes to figure out its comprehensive role in reproduction. Moreover, protein docking analysis is performed to assess its therapeutic potential in case of infertility or its potential use as sex ratio control as well as to regulate the gene expression profile. Hence the present study is designed to explore an important economic aspect in animal reproduction with wide applicability.


Data pre-processing

In present study, we got 11007 and 10445 unique genes in Bos taurus and Sus scrofa, respectively, where we normalized expression data with Limma package. Expression data, before and after normalization can be found in Supplementary Table 1 (GSE47139-before.txt”, “GSE47139-norm.txt”, “xy-before.txt” and “xy-norm.txt) and Supplementary Figure 1. Bos taurus and Sus scrofa expression data shared 555 homologous genes which were listed in Supplementary Table 2.

Consensus modules identification of Bos taurus and Sus scrofa

Definition of adjacency function

In order to satisfy the preconditions of scale-free network distribution, we need to explore the value of the adjacency matrix weight parameter: power. The power plot was shown in Figure Figure2.2. The higher the square of the correlation coefficient is, the closer the network is to the network-free distribution. We choose power = 20 as the power cutoff, when the square of the correlation coefficient reached 0.9 for the first time.

Figure 2
Adjacency matrix weight parameter power plot

Identification consensus modules of Bos taurus

We regarded Bos taurus gene expression data as training dataset, constructed network and mined gene modules based on 555 shared genes in Bos taurus. Firstly, the coefficients of dissimilarity between genes were calculated, and the system clustering tree was obtained. Then cut the clustering tree according to standard for hybrid dynamic shear trees. We set minimum number of genes per module to 30 genes and the pruning height to cutHeight = 0.9. We finally got 9 different modules in Bos taurus gene expression data, shown in Figure Figure3A3A.

Figure 3
Module clustering tree of Bos taurus (A) and Sus scrofa (B). Different colours in module bar mean different modules. D1 and D2 refer to Bos taurus and Sus scrofa datasets, respectively. M1 to M11 mean different module numbers.

Screening consensus modules in Sus scrofa

We also screened gene modules in Sus scrofa by using same method and parameters as in Bos taurus. We obtained 11 gene modules in Sus scrofa, 2 more modules than that in Bos taurus, as shown in Figure Figure3B.3B. Then we calculate consensus of every two modules in Bos taurus and Sus scrofa. The consistency of modules was displayed in Figure Figure4.4. Information for each module and overlapped genes were listed in Table Table1.1. Genes and their corresponding colors of Bos taurus and Sus scrofa can be found in Supplementary Table 3.

Figure 4
Correspondence of Bos taurus-specific modules and Sus scrofa-specific modules
Table 1
Characterization of Bos taurus and Sus scrofa modules

We collected significantly overlapped genes in consensus modules for further analysis, that is to say, a total of 167 genes in column 6 of Table Table1,1, which were listed in Supplementary Table 4.

Gene significance analysis

A gene set with 846 gene which met cutoff (p<0.05 and |log2FC|>0.585) were screened as DEGs in Bos taurus by using Limma package in R language, statistic test volcano plot and heatmap were shown as Figure Figure5A5A and and5B,5B, DEGs were listed in Supplementary Table 5. Then we mapped DEGs into 167 module genes selected previously, we finally found 67 genes which had dual features: significantly expressed between X and Y chromosome transcriptome, as well as significantly overlapped in consensus modules. The 67 dual featured genes were included into further analysis (gene list can be found in Supplementary Table 6).

Figure 5
(A) Volcano plot of statistic test for gene expression in Bos taurus. (B) Heatmap of DEGs in Bos taurus.

Construction of co-expression network of shared genes

We calculated expression correlation coefficient of 67 interested genes and ticked out gene pairs whose correlation coefficient below 0.8 (correlation coefficient matrix was in Supplementary Table 7), finally constructed co-expression network among 67 genes, as shown in Figure Figure6.6. The co-expression network was consist of 58 nodes (27 down-regulated and 31 up-regulated genes) and 497 edges.

Figure 6
Co-expression network of 67 selected genes

Functional and KEGG pathway analysis of the network genes

We enriched genes in the co-expression network to GO and KEGG pathway, and got 66 significantly related GO biology process (BP) annotation and 2 KEGG pathways with cutoff of p<0.05. We extracted 6 GO annotations which were related to reproduction to display, as shown in Figure Figure7,7, listed in Table 2-1. The KEGG pathway analysis list can be found in Table 2-2 The complete table for GO annotations can be found in Supplementary Table 8.

Figure 7
Barplot for reproduction related GO annotations
Table 2-1
Reproduction related GO annotations list
Table 2-2
KEGG pathway annotations list

Related miRNA and TF search

We searched significantly related TF using The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 and got 4 related TFs: ISRE, AP1FJ, RP58, CREL, see Table Table3.3. In addition, we extracted 4 related miRNAs which targeted genes in co-expression network from miRTarBase (2016 version): bta-miR-181a, bta-miR-17-5p, bta-miR-146b, bta-miR-146a, listed in Supplementary Table 9. We integrated miRNA and their target information into co-expression network to construct miRNA-targeted co-expression network, as shown in Figure Figure88.

Figure 8
miRNA-targeted co-expression network
Table 3
Related TF list

It can be observed that RP58, ISRE and APIFJ are significantly related to 30, 27 and 19 genes respectively including PAFAH1B1 and PPP2R2B. PLAU, a gene responsible for creating suitable environment for embryo implantation and SCD, a gene responsible for fatty acids have a miRNA binding sites for bta-miR-181a. LDLR that has a demonstrated role in sperm storage tubules in hen [16] and its lack of presence in female mouse is associated with low fertility [17]. We have found 3 miRNA binding sites including bta-miR-17-5p, bta-miR-146b, bta-miR-146a associated with LDLR that can prove beneficial target for improving fertility. The expression of LDLR, PLAU and SCD in sperm transcriptome suggests it's crucial role in fertilization and activating energy related pathways.

Association of reproductive genes with immunity and synaptic transmission genes

The reproduction have several predominant biological processes like spermatogenesis, sexual development and pregnancy. The genes associated with these processes are co-expressed with genes associated with immunity in respective modules indicating interesting cross talk among reproductive and immune related genes. The reproductive driver genes PPP2R2B and PRM1 lies in brown module along with important immunity related genes including TREM-1, TNF and TRAF-6, whereas PAFAH1B1 lies in blue module along with immune related genes IL4, IL6 etc. in Bos taurus suggesting a cross talk between reproductive and immune related genes.

It is worth noting that the reproductive driver gene PAFAH1B1 is also associated with synaptic transmission along with other genes including MAOA, CACNB4, CHRNA1, FGF2 suggesting its pivotal role as neuro-transmitter that activates the receptors of another neuron indicating functional role in reproductive endocrinology.

Genetic analysis of vital reproduction genes across Bos taurus and Sus scrofa

We found 6 reproduction related GO annotations (listed in Table 2-1) which involved 10 genes: PPP2R2B, PAFAH1B1, PRM1, HMGB2, OVGP1, GPX4, EPOR, PLAU, RPL29, INHBB. Among the 10 genes, PPP2R2B, PRM1 and PAFAH1B1 participated in 5 reproduction related GO annotations at the same time. Therefore, PPP2R2B, PRM1 and PAFAH1B1 are considered as vital reproduction related genes. We did a series of genetic analysis for the selected three genes: PPP2R2B, PRM1 and PAFAH1B1.

Genetic phylogenetic tree of PPP2R2B, PRM1 and PAFAH1B1 in Bos taurus and Sus scrofa

We downloaded protein sequences of these three genes in Bos taurus and Sus scrofa from NCBI database (Supplementary Table 10), and aligned them with ClustalW in Mega6 (18), and made a phylogenetic tree as shown in Figure Figure9.9. It can seen from the phylogenetic tree that genes with the same names in Bos taurus and Sus scrofa have a closer relationship in evolution.

Figure 9
Phylogenetic tree of PPP2R2B, PRM1 and PAFAH1B1 in Bos taurus and Sus scrofa

Consensus sequences of PPP2R2B, PRM1 and PAFAH1B1 in Bos taurus and Sus scrofa

The calculated order of the most frequent residues that is of nucleotides or amino acids found at each position in sequence alignment is termed as consensus sequence. The related sequences in the multiple alignment are compared to each other and similar sequence motifs are calculated. We aligned downloaded protein sequences with ClustalW in Mega6 and found consensus sequences of PPP2R2B, PRM1 and PAFAH1B1 across Bos taurus and Sus scrofa, results were shown in Figure Figure10.10. The genes with the same names in Bos taurus and Sus scrofa have a higher consistency in the sequence.

Figure 10
Consensus sequences of PAFAH1B1 (A) PRM1 (B) and PPP2R2B (C) in Bos taurus and Sus scrofa. Sequences with black background are consensus sequences.

Domains in PPP2R2B, PRM1 and PAFAH1B1

In our present integrated transcriptome study of Bos taurus and Sus scrofa we identified PPP2R2B, PAFAH1B1 and PRM1 as candidate genes involved in spermatid differentiation and development. We then searched domains in PPP2R2B, PRM1 and PAFAH1B1 with InterProScan [19] in EBI, and finally got 3, 1 and 2 domains, respectively. Information of each domain was listed in Table Table4.4. The visualization result are shown in Figure Figure1111.

Figure 11
The visualization result of Domain in PPP2R2B, PRM1 and PAFAH1B1
Table 4
Domain information list of PPP2R2B, PRM1 and PAFAH1B1

“PR55 a domain in PPP2R2B” has diverse functions including substrate recognition, target an enzyme to correct subcellular localization, have a role in spindle cell assembly and centrosome attachment to nuclei. The expression of PPP2R2B hence indicates its critical role in reproduction especially spermatogenesis. Apart from PR55, WD-REPEATS domain is also present in PPP2R2B. PRM1 contains PROTAMINE domain which plays important roles sperm chromatin is substituted for histones during haploid spermatogenic phase. The sperm DNA is packed into a highly condensed, stable and inactive complex by protamines hence maintaining sperm resistant to lethal mutations. LISH and WD-REPEATS domains correspond to PAFAH1B1. In present study we found PAFAH1B1 as an important phosphatase for sperm motility and its fertilizing capability hence regarded as critical factor in male fertility. The 33-residue LIS1 homology (LisH) motif is found in eukaryotic intracellular proteins that is involved in microtubule dynamics, cell migration, nucleokinesis and most importantly in chromosome segregation when we look for its function in reproduction.

Epigenetic predictions for PPP2R2B, PRM1 and PAFAH1B1 in Bos taurus and Sus scrofa

Cytosines in CpG dinucleotides can be methylated to form 5-methylcytosine. In mammals, a methylating cytosine within a gene can change its expression pattern, hence regulating gene epigenetically. We searched CpG islands of PPP2R2B, PRM1 and PAFAH1B1 in Bos taurus and Sus scrofa. In our prediction, we used Methprimer [20] to identify CpG islands. CpG islands are defined as sequence ranges where the Obs/Exp value is greater than 0.6 and the GC content is greater than 60%. We found no CpG islands in PRM1. The results of PPP2R2B and PAFAH1B1 are shown in Figure Figure1212.

Figure 12
CpG islands in PPP2R2B and PAFAH1B1 in Bos taurus and Sus scrofa

Mutations in PPP2R2B, PRM1 and PAFAH1B1

We searched SNPs in PPP2R2B, PRM1 and PAFAH1B1 from dbSNP, shown in Table Table5.5. In addition, we also searched CNVs of each gene from Database of Genomic Variants (DGV), listed in Table Table6.6. A total of five SNPs were found in PPP2R2B including one missense rs11547494, while PAFAH1B1 was found to have three SNPs rs12143448254, rs12143448454, rs12143448654, all of which are missense from dbSNP database. There was no SNP found in PRM1 in dbSNP database.

Table 5
Sequence variations from dbSNP
Table 6
Structural variations from Database of Genomic Variants (DGV)

Protein docking analysis of vital reproduction genes

Protein structures of PAFAH1B1 (PDB ID: 3DT6), PPP2R2B (PDB ID: 4MEW) and PRM1 (PDB ID: 4Y0Q) were downloaded from PDB and docked against MPD3 database ready to dock library using MOE docking tool. Ten conformations were provided by MOE for each ligand. These conformations were arranged according to S score. From 2500 MPD3 database ligands, only top seven individual conformations for every ligand which showed minimum S score and RMSD value were chosen. Along with minimum S score, these selected ligands also showed strong binding with targeted proteins. Among the strongly bonded selected ligands Tannic Acid, 5,6,7,4'-Tetrahydroxyflavanone 6,7-diglucoside, Oolonghomobisflavan-A showed most favorable results with PAFAH1B, PPP2R2B, and PRM1 respectively. The detailed results of docking including information about interacting residues are shown in Table Table77 along with interacting residue figures in Figure Figure1313 where 1A, 1B, 1C are showing protein-ligand complexes (ligand is in magenta color), 2A, 2B, 2C showing 2D docked interactions of ligand-protein and 3A, 3B, 3C are showing binding pocket mode of proteins.

Figure 13
Interaction analys is of docked proteins; (A) (1, 2, 3) PAFAH1B1, (B) (1, 2, 3) PPP2R2B, (C) (1, 2, 3) PRM1.
Table 7
Docking results of vital reproductive genes


The advances in a recent day molecular biology provides ample information suggesting that phenotypes are actually a result of several genes that works in a network. Reproduction is a complex biological process that involves a number of genes and pathways to work together in a coherent way to bring out the desired phenotypes or germ cells capable of fertilization. WGCNA has been widely used to find out candidate or signature genes for disease condition [21]. We have used this state of the art analysis tool to find out reproductive driver genes across Bos taurus and Sus scrofa. Such an analysis has a wide application in animal reproductive biology where one can target specific set of genes as per own desire to generate desired phenotypes for higher reproductive performance, fertility or sex ratio control.

The gene enrichment in co-expression network of genes in present study demonstrates 66 different biological processes including genes related to immunity indicating its connected role in reproduction. It is worth noting that the immunity genes co-express with reproductive genes provide a safe microenvironment for the sperm to develop as well as to fertilize oocyte in germ free environment. Another probable role of those genes could be assumed as generating a suitable environment in female reproductive tract to select specific X or Y bearing spermatozoa [22]. Such genes can be exploited in favor of female or male sex ratio control.

There are mechanisms that favor the X chromosome bearing spermatozoa to fertilize oocyte. One such mechanism is that at time of spermiogenesis, the tail formation of sperm where certain gene expression benefits X bearing spermatozoa can make the female ratio higher than expected 1:1 [5, 23]. In present study PAFAH1B1, PPP2R2B and PRM1 is involved with spermatid development, where PAFAH1B1 have a role with spermiogenesis and sperm maturation. At this stage PAFAH1B1 is responsible for tail formation where it can differentially benefit X bearing spermatozoa to develop tail that can increase the ratio of female from normal 1:1.

It is proposed in Mendel's first law that an XY cross will produce a 1:1 male to female ratio in sexually reproducing organisms as evolutionary forces favor parents that spend equal amount of resources on both sexes and an imbalance in this ratio will immediately counter balance the effect in favor of under-presented sex [24]. It has been observed in a population of metazoans that it transmits its sex chromosomes unequally that cause broods to produce highly skewed sex ratio and was termed as a result of asymmetric spermatocyte division [23]. There are three possible ways for distorted sex ratio in genetic terms 1) non random sex chromosome segregation at meiotic division [25] 2) post meiotic differences in haploid gametes causes sperm specific function differences in response to sex chromosomes [26] or 3) sex specific embryonic lethality caused by post zygotic mechanisms [27]. Recently two spermatogenesis based modifications to the cellular program are reported in nematode Rhabditis sp. sB347 that sire >5% male progeny. Firstly the meiosis is modified where sister chromatids of unpaired X chromosome separate prematurely at meiosis 1 stage, and secondly during anaphase II of meiosis II, cellular components that are essential for sperm motility are partitioned exclusively to X-bearing sperm [23]. In present study, we found that the PRM1, PPP2R2B and PAFAH1B1 are associated with spermatogenesis specially the LiSH domain containing protein PAFAH1B1 have role in chromosome segregation and can be used as target for sex ratio increase in favor of female. More interestingly several genes that are associated with motility are co-expressed with reproductive driver genes suggesting its combined role in producing healthy sperm capable of fertilization; and those genes can be made exclusively partitioned proteins responsible for motility to X bearing spermatozoa in controlled breeding programs. Furthermore, specific alterations in gene expression can partition the components of essential sperm motility exclusively to X or Y chromosome hence change normal sex ratio.

Sperm needs to undergo many post meiotic changes of the spermatogenic process before it starts its journey into the reproductive tract of female. The haploid spermatid undergo striking morphological changes including restructuring and compaction of its chromatin material. This makes the almost complete substitution of histone based nucleosomal and histone based structure with protamine based structure. Protamine protects the paternal genetic material and remains key factor in male fertility [28]. In our present study PRM1 is differentially expressed and is present in reproductive module and hence is an important reproductive driver gene. The consequence of mutations in protamine based chromatin can result in infertility in mammals, in our results, the dbSNP search shows no SNP records in PRM1, although we found two structural variations from Database of Genomic Variants (DGV) showing mutations in PRM1 gene is highly protected to safeguard paternal genetic material.

To further elucidate the important roles of PRM1, PPP2R2B, PAFAH1B1 we made genetic analysis of these genes. We made a phylogenetic tree, predicted important domains, carried out epigenetic analysis and mutation analysis. Furthermore, to regulate the expression of these gene we performed protein docking so that we can recommend chemicals or food items that are beneficial or harmful for reproduction.


The Holstein sperm was separated by flow cytometry and after achieving purity of 95% and stored at -80 C. Then we extracted RNA according to manual performed RNA quality testing and then make the microarray geneChip.

Data and method

Data collection and pre-processing

The hybridized microarrays results in images that were used to generate datasets which were pre processed to get meaningful data sets for further analysis. A particular type of preprocessing is termed as normalization that accounts for systematic difference across different datasets. To identify reproductive driver genes, two datasets were included in present study: 1) X and Y chromosome Affimetrix transcriptome data of Bos taurus generated in our lab (only one sample for each chromosome). 2) NCBI raw transcriptome data of Sus scrofa with accession number of GSE47139 [22] was downloaded from GEO website for analysis (a total of eight samples in GSE47139, 4 X and 4 Y chromosome Affimetrix transcriptomes of Sus scrofa.). Both datasets were pre-processed by using R-language Affy package [29], including background correction and normalized simultaneously. The data sets are all normalized by using Limma package in R language.

Consensus modules identification of Bos taurus and Sus scrofa

We constructed co-expression network and mined reproductive related modules in Bos taurus and Sus scrofa respectively, based on weighed gene co-expression network analysis (WGCNA) algorithm, which is a typical system biology algorithm for constructing gene co-expression network [30]. We extracted overlapped genes in Bos taurus and Sus scrofa for WGCNA, the network building and module identification process steps in Bos taurus transcriptome were as follows:

Definition of gene expression correlation matrix

Correlation coefficient of gene m and gene n was defined as:

Smn=|cor(m,n)|, We calculated expression correlation coefficient of every two genes to make a whole correlation coefficient matrix.

Definition of adjacency function

WGCNA used power index adjacency function as measurement of correlated indicators, the function was defined as:

amn, where, amn=power(Smn,β), then determined weighted factor β, we choose β above 0.9 as a power cutoff.

Node dissimilarity measurement

The adjacency matrix amn was transformed into a topological matrix Ω=wmn, where


lmn referred to adjacency coefficients sum of all connected nodes, km represented connection sum of gene m, dissimilarity of genes was defined as dmn, dmn = 1 - wmn.

Module identification

We constructed hierarchical clustering tree based on dmn defined before, tree branches represented different gene modules.

Screening consensus modules in Sus scrofa

We also identified reproductive related modules in Sus scrofa based on WGCNA with same processes and cutoffs as in Bos taurus described above. Then we calculated overlapped genes of each consensus modules in Bos taurus and Sus scrofa, the selected overlapped genes in consensus modules were used for further analysis.

Gene significance analysis

Filtered expression data of Bos taurus were analyzed with the test method in Limma package [31], implemented in R language, Bioconductor project. Meanwhile, fold change between groups were also calculated, because each sample may show intrinsic individual variability, the threshold for determining the fold change (FC) was set at 1.5. Significantly differentially expressed genes (DEG) were defined with the cutoff of p value <0.05 and |log2FC|>0.585. Then we selected shared genes by comparing DEGs and selected overlapped genes in consensus modules. The shared genes were used for further analyses.

Construction of co-expression network of shared genes

We calculated gene expression correlation coefficient of every two selected shared genes, and kept gene pairs with correlation coefficient higher than 0.8 to construct gene co-expression network. Network was visualized by Cytoscape 2.8.0 software [32].

Functional and KEGG pathway analysis of the network genes

To further analyze the functions of network genes, we enriched genes in the co-expression network to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways by using clusterProfiler [33] with a cutoff of p lower than 0.05.

Related miRNA and TF search

We searched significantly related TF by using The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 [34], and extracted related miRNAs of genes in co-expression network from miRTarBase (2016 version) [35]. Then constructed TF/miRNA -targeted co-expression network by integrating miRNA/TF and co-expression information.

Genetic analysis of vital reproduction genes across Bos taurus and Sus scrofa

Phylogenetic tree was made using Mega 6 after aligning the downloaded sequences in ClustalW.We searchedwith InterProScan [19] in EBI, we identify CpG islands using Methprimer [20], and finally searched SNPs and CNVs from dbSNP and DGV respectively.

Protein docking analysis of vital reproduction genes

Molecular docking analysis were performed using MOE (Molecular Operating Environment) software [36]. Protein structures of vitally important genes were downloaded from PDB (Protein Data Bank) [37] and optimized by minimizing their energies with parameters (Force Field: AMBER99, Gradient: 0.05). The minimized structures were used as the receptor protein for docking. MOE site finder tool was used to find out the active site where ligands can bind. To select a suitable ligand, selected proteins were docked against ready-to-dock library of MPD3 database [38] contacting 2500 potential compounds in 3D ready to dock format. MOE docking program with parameters (Rescoring function: London dG, Placement: Triangle matcher, Retain: 10, Refinement: Force field, Rescoring 2: London dG) was used to bind the selected ligands with proteins.



The authors extend sincere thanks to laboratory colleagues for their beneficial suggestions.


Differentially expressed genes
Kyoto Encyclopedia of Genes and Genomes
biological process
Gene Ontology
Weighted Gene Co-Expression Network Analysis
Protein Data Bank
Molecular Operating Environment
Trancription Factor
Contributed by

Author contributions

F.A.K. and S.Z. conceived and designed the experiments; F.A.K performed experiments bioinformatics data analysis; H.L. develop the microarray data of Bos taurus. F.A.K and N.S.P. wrote the paper; M.T.Q. performed protein docking analysis and F.A.K., N.S.P., H.Z., K.W. and S.Z. revised the paper, and all the authors approved the final version.


The authors declare that there is no potential conflicts of interest.


The authors are grateful for financial assistance from the project (The research platform construction of the ministry of agriculture international exchange and cooperation in science and technology) and EUFP7 projects (PIIFR-GA-2012-912205 and FP7-KBBE-2013-7-613689). The authors are grateful to all the laboratory colleagues constructive discussions.


1. Graves JAM, Ferguson-Smith MA, McLaren A, Mittwoch U, Renfree MB, Burgoyne P. The evolution of mammalian sex chromosomes and the origin of sex determining genes [and Discussion] Philosophical Transactions of the Royal Society of London B: Biological Sciences. 1995;350:305–312. [PubMed]
2. Graves JAM. Evolution of the mammalian Y chromosome and sex-determining genes. The Journal of experimental zoology. 1998;281:472–481. [PubMed]
3. Graves JAM. Sex chromosome specialization and degeneration in mammals. Cell. 2006;124:901–914. [PubMed]
4. Vallender EJ, Lahn BT. How mammalian sex chromosomes acquired their peculiar gene content. Bioessays. 2004;26:159–169. [PubMed]
5. Wu JC, Go AC, Samson M, Cintra T, Mirsoian S, Wu TF, Jow MM, Routman EJ, Chu DS. Sperm development and motility are regulated by PP1 phosphatases in Caenorhabditis elegans. Genetics. 2012;190:143–157. [PubMed]
6. Sassone-Corsi P. Unique chromatin remodeling and transcriptional regulation in spermatogenesis. Science. 2002;296:2176–2178. [PubMed]
7. Tanaka H, Baba T. Gene expression in spermiogenesis. Cell. Mol. Life Sci. 2005;62:344–354. [PubMed]
8. Miller D, Ostermeier GC. Towards a better understanding of RNA carriage by ejaculate spermatozoa. Hum. Reprod. Update. 2006;12:757–767. [PubMed]
9. Varmuza S, Jurisicova A, Okano K, Hudson J, Boekelheide K, Shipp EB. Spermiogenesis is impaired in mice bearing a targeted mutation in the protein phosphatase 1cγ gene. Developmental biology. 1999;205:98–110. [PubMed]
10. Chu DS, Liu H, Nix P, Wu TF, Ralston EJ, Yates III, JR, Meyer BJ. Sperm chromatin proteomics identifies evolutionarily conserved fertility factors. Nature. 2006;443:101–105. [PMC free article] [PubMed]
11. Hilscher B. Spermatogoniogenesis, an interacting proliferation process between stem cell spermatogonia and differentiating spermatogonia. Fortschritte der Andrologie. 1981;7:46–57.
12. Grootegoed JA, Jansen R, Van Der Molen HJ. The role of glucose, pyruvate and lactate in ATP production by rat spermatocytes and spermatids. Biochimica et Biophysica Acta (BBA)-Bioenergetics. 1984;767:248–256. [PubMed]
13. Chocu S, Calvel P, Rolland AD, Pineau C. Spermatogenesis in mammals: proteomic insights. Syst. Biol. Reprod. Med. 2012;58:179–190. [PubMed]
14. Oakberg EF. Duration of spermatogenesis in the mouse and timing of stages of the cycle of the seminiferous epithelium. Am J Anat. 1956;99:507–516. [PubMed]
15. Schenk JL, Cran DG, Everett RW, Seidel GE. Pregnancy rates in heifers and cows with cryopreserved sexed sperm: Effects of sperm numbers per inseminate, sorting pressure and sperm storage before sorting. Theriogenology. 2009;71:717–728. [PubMed]
16. Huang A, Isobe N, Obitsu T, Yoshimura Y. Expression of lipases and lipid receptors in sperm storage tubules and possible role of fatty acids in sperm survival in the hen oviduct. Theriogenology. 2016;85:1334–1342. [PubMed]
17. Guo T, Zhang L, Cheng D, Liu T, An L, Li WP, Zhang C. Low-density lipoprotein receptor affects the fertility of female mice. Reproduction, Fertility and Development. 2015;27:1222–1232. [PubMed]
18. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol. 2013;12:2725–9. [PMC free article] [PubMed]
19. Quevillon E, Silventoinen V, Pillai S, Harte N, Mulder N, Apweiler R, Lopez R. InterProScan: protein domains identifier. Nucleic Acids Res. 2005;33:W116–20. [PMC free article] [PubMed]
20. Li LC, Dahiya R. MethPrimer: designing primers for methylation PCRs. Bioinformatics. 2002;18:1427–31. [PubMed]
21. Sundarrajan S, Arumugam M. Comorbidities of Psoriasis-Exploring the Links by Network Approach. PloS one. 2016;11:e0149175. [PMC free article] [PubMed]
22. Almiñana C, Caballero I, Heath PR, Maleki-Dizaji S, Parrilla I, Cuello C, Gil MA, Vazquez JL, Vazquez JM, Roca J, Martinez EA, Holt WV, Fazeli A. The battle of the sexes starts in the oviduct: modulation of oviductal transcriptome by X and Y-bearing spermatozoa. BMC Genomics. 2014;15:293. [PMC free article] [PubMed]
23. Shakes DC, Neva BJ, Huynh H, Chaudhuri J, Pires-daSilva A. Asymmetric spermatocyte division as a mechanism for controlling sex ratios. Nature communications. 2011;2:157. [PubMed]
24. Robert T. Parental investment and sexual selection. Sexual Selection & the Descent of Man, Aldine de Gruyter. 1972. pp. 136–179.
25. Kubai DF. Meiosis in Sciara coprophila: structure of the spindle and chromosome behavior during the first meiotic division. J. Cell Biol. 1982;93:655–669. [PMC free article] [PubMed]
26. Cazemajor M, Joly D, Montchamp-Moreau C. Sex-ratio meiotic drive in Drosophila simulans is related to equational nondisjunction of the Y chromosome. Genetics. 2000;154:229–236. [PubMed]
27. Orr HA. Haldane's rule. Annu. Rev. Ecol. Syst. 1997;28:195–218.
28. Gunes S, Al-Sadaan M, Spermatogenesis Agarwal A. DNA damage and DNA repair mechanisms in male infertility. Reproductive biomedicine online. 2015;31:309–319. [PubMed]
29. Liu WM, Mei R, TB Di X Ryder, Hubbell E, Dee S, Webster TA, Harrington CA, Ho MH, Baid J, Smeekens SP. Analysis of high density expression microarrays with signed-rank call algorithms. Bioinformatics. 2002;18:1593–1599. [PubMed]
30. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. doi: 10.1186/1471-2105-9-559. [PMC free article] [PubMed] [Cross Ref]
31. Smyth G. Gentleman VCR, Dudoit S, Irizarry R, Huber W, editors. Limma: linear models for microarray data In ‘Bioinformatics and Computational Biology Solutions using R and Bioconductor’ 2005. pp. 397–420.
32. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011;27:431–2. [PMC free article] [PubMed]
33. Yu G, Wang LG, Han Y, He QY. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7. [PMC free article] [PubMed]
34. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID Bioinformatics Resources. Nature Protoc. 2009;4:44–57. [PubMed]
35. Chou CH, Chang NW, Shrestha S, Hsu SD, Lin YL, Lee WH, Yang CD, Hong HC, Wei TY, Tu SJ, Tsai TR, Ho SY, Jian TY, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016;44:D239–47. [PMC free article] [PubMed]
36. Vilar S, Cozza G, Moro S. Medicinal chemistry and the molecular operating environment (MOE): Application of QSAR and molecular docking to drug discovery. Curr Top Med Chem. 2008;8:1555–1572. [PubMed]
37. Bernstein FC, Koetzle TF, Williams GJ, Meyer EF, Brice MD, Rodgers JR, Kennard O, Shimanouchi T, Tasumi M. The protein data bank. European Journal of Biochemistry. 1977;80:319–324. [PubMed]
38. Mumtaz A, Ashfaq UA, ul Qamar MT, Anwar F, Gulzar F, Ali MA, Saari N, Pervez MT. MPD3: a useful medicinal plants database for drug designing. Natural Product Research. 2016. pp. 1–9. [PubMed]

Articles from Oncotarget are provided here courtesy of Impact Journals, LLC