|Home | About | Journals | Submit | Contact Us | Français|
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.
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 .
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.
The paternal genetic material across different species is transferred to oocyte by sperm after several dramatic morphological and physical changes . 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 . 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 . 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 . 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 .
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 , meiosis that generates four haploid spermatids  and spermiogenesis a process of dramatic morphological changes from spermatid to mature spermatozoa . 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 .
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 . 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.
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.
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.
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.
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.
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).
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.
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.
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.
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  and its lack of presence in female mouse is associated with low fertility . 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.
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.
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.
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.
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.
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  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.
“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.
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  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.
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.
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.
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 . 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 . 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 . 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 . There are three possible ways for distorted sex ratio in genetic terms 1) non random sex chromosome segregation at meiotic division  2) post meiotic differences in haploid gametes causes sperm specific function differences in response to sex chromosomes  or 3) sex specific embryonic lethality caused by post zygotic mechanisms . 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 . 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 . 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.
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  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 , including background correction and normalized simultaneously. The data sets are all normalized by using Limma package in R language.
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 . 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:
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.
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.
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.
We constructed hierarchical clustering tree based on dmn defined before, tree branches represented different gene modules.
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.
Filtered expression data of Bos taurus were analyzed with the test method in Limma package , 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.
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 .
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  with a cutoff of p lower than 0.05.
We searched significantly related TF by using The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 , and extracted related miRNAs of genes in co-expression network from miRTarBase (2016 version) . Then constructed TF/miRNA -targeted co-expression network by integrating miRNA/TF and co-expression information.
Phylogenetic tree was made using Mega 6 after aligning the downloaded sequences in ClustalW.We searchedwith InterProScan  in EBI, we identify CpG islands using Methprimer , and finally searched SNPs and CNVs from dbSNP and DGV respectively.
Molecular docking analysis were performed using MOE (Molecular Operating Environment) software . Protein structures of vitally important genes were downloaded from PDB (Protein Data Bank)  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  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.
Author contributionsF.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.
CONFLICTS OF INTEREST
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.