Transcriptomic analysis of reproductive damage in the epididymis of male Kunming mice induced by chronic infection of Toxoplasma gondii PRU strain

Background Some researchers have reported that Toxoplasma gondii can cause serious reproductive impairment in male animals. Specifically, T. gondii destroy the quality of sperm in the epididymis, which affects their sexual ability. However, among such studies, none have investigated the male reproductive transcriptome. Therefore, to investigate the relationship between T. gondii and sperm maturation, we infected mice with T. gondii prugniaud (PRU) strain and performed transcriptome sequencing of the epididymis. Results Compared with the control group, 431 upregulated and 229 downregulated differentially expressed genes (DEGs) were found (P-value < 0.05, false discovery rate (FDR) < 0.05 and |log2 (fold change)| ≥ 1). According to results of a bioinformatics analysis, Gene Ontology (GO) function is divided into three categories: cellular component, molecular function and biological process. Upon performing GO analysis, we found that some DEGs correlated with an integral part of membrane, protein complex, cell surface, ATP binding, immune system process, signal transduction and metabolic process which are responsible for the epididymal injury. DEGs were mapped to 101 unique KEGG pathways. Pathways such as cytokine-cytokine receptor interaction, glycolysis/gluconeogenesis and apoptosis are closely related to sperm quality. Moreover, Tnfsf10 and spata18 can damage the mitochondria in sperm, which decreases sperm motility and morphology. Conclusions We sequenced the reproductive system of male mice chronically infected with T. gondii, which provides a new direction for research into male sterility caused by Toxoplasma infection. This work provides valuable information and a comprehensive database for future studies of the interaction between T. gondii infection and the male reproductive system.

the epididymis is important for identification of factors affecting sperm maturation. Previous studies have shown that reproductive pipeline damage [4] and hypogonadism [5] are associated with T. gondii infection. Furthermore, a large number of tachyzoites in semen of infertile patients and anti-sperm antibodies (AsAb) were detected in T. gondii infection [6]. Other studies have reported more advanced pathological changes in T. gondii infection cases, such as granular degeneration of vas deferens epithelial cells, luminal stenosis and interstitium infiltration with inflammatory cells [7]. The density and motility rate of sperm in epididymis infected with T. gondii were significantly lower than the control group; however, the rate of sperm deformity was increased [8,9].
Recently, various -omics technologies, including transcriptomics, proteomics and metabolomics, have been developed [10]. Transcriptomic analysis is one of the most common high-throughput techniques, which identify the types and copy numbers of mRNAs, while the cells are living in a functional state [11]. Genetic studies have shown that mRNA act as a "bridge" for genetic information transmission between DNA and protein.
Hence, it is a valuable source to identify the expression of all genes with a specific time and space in the cell.
Previous studies have investigated the association between T. gondii and male infertility using epidemiology, pathology and serology methods and techniques [12], detection of Toxoplasma DNA in semen [13] and observation of the tissue and cell damage [14]. To our knowledge, studies investigating the differentially expressed genes of host and T. gondii by transcriptome sequencing of the male reproductive system are limited. The main objective of this study was to examine differential gene expression by RNA sequencing (RNA-seq) technology in order to identify key genes associated with T. gondii (PRU strain) chronic infection of the epididymis in male Kunming mice.

Study population and experiment set-up
Thirty specific-pathogen-free eight-week-old Kunming male mice were purchased from the Laboratory Animal Center of Guangdong Province. Fifteen mice were subjected as the experiment groups, and 15 mice were subjected as the control group (to reduce individual differences, we set the epididymal tissue of 5 mice as a biological replicate, the experiment group and the control group were repeated three times). Mice in the experimental group were inoculated with four cysts of T. gondii PRU strain diluted with normal saline to 0.5 ml through the intragastric administration route (we have previously studied the optimal number of PRU strains infections per mice, which tried to minimize the case damage caused by T. gondii infection in mice). Meanwhile, the control groups were given the same amount of normal saline only. The male mice were sacrificed at 35 days post-infection. Under sterile conditions, the epididymides were harvested. Under the microscope, the peripheral adipose tissue and blood vessels of the harvested epididymides were carefully removed. The processed epididymides were subjected to quick freezing by storing them in liquid nitrogen at -80 °C for subsequent analysis.

Transcriptome sequencing, data analysis and verification
At present, transcriptome sequencing is a mature highthroughput second-generation sequencing method. A brief description of the experimental process is provided in Table 1; detailed methods of transcriptome sequencing, data analysis and verification can be found in Additional file 1: Text S1.

Pre-experiment
Before the histological analysis, we conducted a preliminary study on the sperm quality of T. gondii infection, including the statistics of sperm survival rate The sperm count per milliliter semen was calculated according to the formula: Sperm count per milliliter semen = (5 medium square count/80) × 400 × dilution multiple × 10,000. We diluted sperm 10 times for microscopic examination. The experimental data are shown in Additional file 2: Figure S1; there were significant differences in total sperm count between the experimental group and the control group (P = 0.0058). The data showed that T. gondii infection had a negative effect on sperm quality.

Identification of infection in mice
The outward appearance changes of infected mice comprised mainly a stress reaction about a week after infection, which included depression and the fur is not smooth. Thirty days after infection, 20 mice were selected for tail vein blood collection and qPCR was used to verify the infection of mice (T. gondii ROP18 virulence gene was detected and using the β-actin as a reference); the result is shown in Additional file 3: Figure S2. All of the 20 mice were infected. Fifteen mice were randomly selected for dissection and epididymis tissue was quickly removed. Meanwhile, the mouse brain was grinded and diluted with 100 µl saline solution; 5 µl of brain homogenate was taken for microscopic examination and the number of cysts was counted. Cysts from each mouse were counted three times and averaged. Statistical results showed that 5 µl brain homogenate contained three or more cysts. Therefore, we determined that all of the experimental mice were infected with T. gondii.

Quality analysis of RNA
The purity of total RNA was high, and the ratios OD 260 / OD 280 and OD 260 /OD 230 were both higher than 2.0 ( Table 2). The RNA integrity number (RIN) was higher than 7, and the ratio of 28S/18S rRNA was higher than 0.7. This indicated that the RNA had good integrity and met the requirements of the subsequent experiments.

Sequencing data preprocessing
The sample library constructed by Illumina Hiseq 2500 sequencing platform. We obtained 321,183,388 raw reads (48.17G). After trimming, we obtained 314,632,244 clean reads (47.21G). Among them, the amounts of sequencing data of each sample were more than 6G. In high throughput sequencing, each measurement of a basic group gives a corresponding value, which is a measure of the accuracy ( Table 3).

Analysis of gene expression
RPKM/FPKM was used to measure the abundance of gene expression in the transcriptome. We used FPKM to estimate the abundance of known genes expresses in different samples. After the construction of new transcripts, 69,907 transcripts and 29,128 unigenes were obtained.  The five-demographic data (distribution) of experimental and control groups are shown in Fig. 1.

Analysis of DEGs
Significant DEGs from the differential genes with P < 0.05, |log2 (fold change)| ≥ 1 and the FDR corrected P-value < 0.05 were selected. The overall distribution of DEGs was observed by drawing a volcano plot (Fig. 2). Compared with the control group, the number of up-and downregulated genes was 431 and 229, respectively. Partial DEGs are listed in Tables 4 and 5.

Bioinformatics analysis for the DEGs
GO database analysis was conducted on the DEGs, including biological processes (1746), cellular components (1275) and molecular functions (785); the distribution diagram of GO terms with significant changes in classification is shown in Fig. 3. The main categories in the biological processes include immune system process, response to stress and signal transduction. As for the cellular components, the GO terms with the most differentially expressed genes were cell, intracellular and organelle. The molecular functions were involved mainly in the protein binding followed by ion binding, and transmembrane transporter activity. KEGG database analysis showed that the DEGs were enriched in 101 pathways. The most enriched pathways of DEGs are exhibited in Fig. 4. Markedly, cell adhesion (31), phagosome (27) and natural killer cell mediated cytotoxicity (17) were the major three enriched pathways (Fig. 4). However, the pathway directly related to male reproduction was not found.

Verification results by qPCR analyzing and western blotting
The eight genes related to reproduction were identified by qPCR (Table 6); these expression trends were coincident with RNA-seq results (Fig. 5). The result of qPCR did not coincide with the RNA-seq perfectly; however, it is still possible to increase the reliability of these genes to a certain extent. DEGs including Piwil2, Spata18 and Tnfsf10 were selected for Western blot analysis. The results of the Western blot showed the same changing trend as the results of RNA-seq, although the specific amount of expression was different (Fig. 6).

Discussion
We investigated the pathogenesis of T. gondii in the male reproductive tract of mice at the stage of sperm maturation using RNA-Seq technology. Additionally,  we performed enrichment methodologies analysis to screen the DEGs that might be related to the reproductive functions. Results were confirmed by using qPCR and Western blot methods. Our findings showed that the RNA-Seq is accurate.
The analysis of DEGs suggested that some genes were interrelated with the infection of T. gondii, such as Iigp1 (interferon-inducible GTPase 1), Tgtp1 (T-cell-specific guanine nucleotide triphosphate-binding protein 1) and Tgtp2. Iigp1 found in the parasitophorous vacuoles (PV) plays a role in the destruction of PV and death of the parasite, when the Iigp1 gene is overexpressed [15]. We found that Iigp1 was in regulation of autophagy in the GO classification. Hence, it is possible that Iigp1 is involved in autophagy, although the specific mechanism is not clear. During infection of the type II of T. gondii, Tgtp1 and Tgtp2 were recruited to the PV membrane, PV vesiculation, rupture and digestion of the parasite within the cytosol [15,16]. In our study, the Iigp1, Tgtp1 and Tgtp2 were significantly upregulated by 94-fold, 48-fold and 23-fold, respectively. This could be the result of host resistance to T. gondii. However, after T. gondii infection, Tgtp1 and Tgtp2 had low levels of gene expression abundance. The average abundance of Tgtp1 was 13.7 for the experimental group and 0.58 for the control group and the average abundances of Tgtp 2 were 15.39 and 0.31, respectively. Tgtp1 and Tgtp2 are mainly expressed in the thymus and lymph nodes, predominantly T-cells and IFNG-stimulated macrophages [17]. In other words, these genes are tissue-specific. After T. gondii infection, the immune system was stimulated, and the expression level of these genes was greatly increased.

Glycolytic pathway and sperm motility
There are several enriched pathways in our data analysis. Glycolysis was thought to be the oldest and most primitive way of biological energy harvesting. We focused on DEGs that affect glycolysis such as phosphoglycerate kinase 2 (Pgk2), phosphoglycerate mutase (Pgam) and lactate dehydrogenase C (LDHC). Phosphoglycerate kinase is the first regulation factor of ATP generated in the glycolysis pathway. In the reproductive system, the  Spink2 Serine protease inhibitor Kazal type 2 − 1.02159 Required for maintenance of normal spermatogenesis, involved in regulating serine protease-dependent germ cell apoptosis examination of sperm in elderly and young male patients with asthenospermia showed that the level of Pgk2 was reduced. Pgk2, which is related to the quality of sperm, is mainly expressed in able-bodied spermatids [18]. Pgk2 deficiency leads to a decrease in sperm motility and ATP levels [19]. Pgam is an important enzyme of glycolysis and gluconeogenesis pathway, catalyzes the 3-phosphoglycerate (3-PGA) into 2-phosphoglycerate (2-PGA), which dominates the balance between glycolysis and gluconeogenesis. Moreover, Pgam2 is related to the fibrous sheath of the main segment in the sperm flagellum [20]. The fibrous sheath on the sperm flagellum is the source of ATP and an essential structure of sperm motility [21]. Therefore, Pgam2 is the key determinant to maintaining sperm motility and morphology. LDHC is an enzyme that controls the conversion of pyruvate to L-lactic acid, which plays a key role in the process of spermatogenesis and capacitation by participating in cellular  Fig. 4 The most enriched pathways of DEGs energy metabolism [22]. LDHC was found in the human spermatozoa and spermatogenic cells and was increased gradually with meiosis. LDHC in mature sperm is mainly located in the fibrous sheath of the main segment of the tail of sperm, which is involved in the regulation of sperm motility. In addition, LDHC is also participates in the homeostasis of ATP. Deficiency/destruction of germlinespecific LDHC leads to infertility, including a reduction in sperm ATP content and reduced sperm motility [22,23]. Conclusively, the enzymes are closely linked with sperm activity and morphology.
In the present study, the DEGs of glycolysis indicated that T. gondii infection might also result in dysregulation of those enzymes in humans, thus promoting a negative influence on the quality of male sperm. Normally, the energy production efficiency of oxidative phosphorylation is much higher than that of glycolysis, but glycolytic pathways have their particular importance in sperm energy sources. Mukai & Okuno [24] noticed that when the oxidative phosphorylation was blocked by rotenone of mitochondrial respiratory chain inhibitors, the sperm keeps moving. Therefore, glycolysis is not only an important source of energy for sperm but also plays a role in sperm capacitation [25] and acrosome reaction [26,27].

Mitochondrial damage and apoptosis
Spermatogenesis associated 18 (Spata18) is a key regulator of mitochondrial quality that mediates repair or degradation of unhealthy mitochondria [28]. There is no pathway for Spata18, but previous research has shown that it is a monitor of cell differentiation process, which is involved in the maturation of spermatids into spermatozoa [29].
In our study, the expression of Tnfsf10 gene was upregulated to 2.89-fold in the pathway of cytokine-cytokine receptor interaction, apoptosis and natural killer cell Table 6 Design of primers used for qPCR analysis   Gene ID  Gene  Forward primer  Reverse primer  Product length (bp)   57746  piwil2  GGC AGA GGC CTT GTG TTT AG  CGT TTT GAA GGA GGC CGG AA  134   73472  Spata18  TTA TCA CGT GTG GCC TGC TC  TTA AGC TCC TGC TAC GGC AC  125   22035  Tnfsf10  ACG CTT CCA AGA TGG TCT CA  ACA GTC CGT ACT CGG CAT CT  147   21822  Tgtp1  CGC CAC GTC TTC TCA CTG TC  CAC CAA GTG GAA TGG TGG CT  130   60440  Iigp1  TCA CTA TGA CTT CCC CGT CCT  TGC TTC AGA AAT TGC CGC TT  132   18663  Pgk2  GCC TGT [30]. Apoptosis, or programmed cell death (PCD), is a biological process controlled by genes and plays a vital role in maintaining the morphology of normal tissues and organs. So far, the mitochondrial pathway is considered to be the classical pathway of apoptotic signal pathways. Mitochondria control the life and death of cells in aerobic environments, which is the main place of ATP production, and the major regulatory site of cell apoptosis. Mitochondrial permeability transition pores (MPTP) exist in the inner and outer membranes of mitochondria, and the permeability transition plays a vital role in cell apoptosis [31]. The pathway of cytokinecytokine receptor interaction shows Tnfsf10 as a cytokine and interacts with the cytokine receptor Tnfsf10b. In the apoptosis pathway, Tnfsf10b (2.3-fold upregulation) as a death receptor, interacts with Tnfsf10 causing production of the death inducing signaling complex (DISC), which initiates the caspaseprotease family (aspartate-specific cysteine proteases) and activates apoptosis through DNA fragmentation, low synthesis of ploy, and destruction of nuclear membrane integrity [32]. At the same time, the translocation and permeabilization of mitochondrial membrane could induce apoptosis through Bax/Bak [33]. MPTP opening resulting from exogenous injury factors (such as Ca 2+ overload) would lead to necrosis and apoptosis of cells. Some scholars have found that the synergism of melittin with Trail could increase the number of apoptosis in hepatoma cell through activating Ca 2+ / calmodulin-dependent protein kinase [34]. This is consistent with the finding of Chen et al. [35] who stated that the Ca 2+ dependent calpain and/or Ca 2+ independent cathepsins participate in apoptosis-like PCD and necrosis-like PCD. Trail gene-deficient (Trail2/2) mice were significantly reduced by 54% compared to the wild type mice, which explains the negative effects on sperm maturation in the reproductive system [30]. A high concentration of soluble Tnfsf10 was detected in the seminal plasma, which revealed a variable expression of Trail receptors in the sperm cellular fraction among different participants [36]. Those results suggest that Tnfsf10 may play an important role in sperm maturation. The cytotoxicity of natural killer (NK) cells is motivated by the activating cell receptors. It can switch on the cascade reaction of downstream signal transduction. In our study, the Src family kinase (Fcgr4) was related to the cytotoxic pathway which is mediated by NK cells. The Fcgr4 also could active the mitogen activated protein kinase (MAPK) pathway, which could produce the cytotoxicity and secret proinflammatory cytokines (such as IFN-γ, TNF-α and GMCS). The upregulation of piwi and Itgb2 genes mediates the target cell apoptosis by Tnfsf10, which is caused by perforin and granzyme. The Gzmb (Granzyme B) gene is responsible for aspartate-specific cysteine proteases (DNA fragmentation, low synthesis of ploy and destruction of nuclear membrane integrity) which cause apoptosis.
From these pathways, we consider that Tnfsf10 is directly or indirectly connected with the male reproductive system. Chronic infection of T. gondii gives rise to dysregulation genes in these pathways, which in turn have a negative consequence on the integrity and motility of sperm. At the stage where spermatids mature into spermatozoa, mitochondria formed the mitochondrial sheath by gathering around the tail and wrapping around the shaft spirally. Finally, it can provide the energy for rotational motion of the tail in the way of making full use of energy in plasma drops. The fertility was influenced by various channels to create the apoptosis of germ cells.

Testis-specific genes in the epididymis
Our data analysis revealed that many DEGs were upregulated in the epididymis; however, previous studies have shown that they are testicular-specific genes, such as Piwil2 and RGS22 (regulator of G-protein signaling 22) genes. Piwil2 is a member of the Argonaute family and PIWI (the P-element-induced wimpy testis) subfamily. PIWI family proteins are only expressed in animal gonad tissues and have tissue specificity, these are closely related to the maintenance of germ stem cells and spermatogenesis [37]. The Piwil2 gene, expressed in tumor stem cells and germ cells, plays an  [38]. Similarly, Rgs22 is a specific gene in the male testis, which can be expressed simultaneously in a variety of epithelial tumors. In testis, Rgs22 protein is associated with spermatogenesis, which can be expressed in spermatogenic cells involved in the spermatogonia differentiation into sperm cells and Sertoli cells. Hu et al. [39] found that the expression of Rgs22 was downregulated in patients with azoospermia. Proteins are located in spermatogenic cells and Leydig cells, which interact with the guanine-nucleotide-binding proteins [39]. Thus, we consider that Rgs22 may be related to the meiotic stage of spermatogenesis.
In summary, Piwil2 and Rgs22 genes are testis-specific genes, mainly expressed in spermatogenic cells and play a key role in the spermatogenesis process. Other genes that were identified (see Table 7) are testis-specific and highly expressed in the epididymis. One possible reason for the latter observation is the invasion of T. gondii, which contributes to the destruction of the blood-testis barrier of the testis, and disruption of the spermatogenesis homeostasis [40]. Furthermore, spermatogenic cells in the seminiferous tubules were disordered, resulting in the increase of testis-specific gene expression in the epididymis. However, the epididymis has its tissue-specific microenvironment [41,42] which is different from that in the testis. The microenvironment in the epididymis is not appropriate to regulate the cells in the testis and therefore, the specific mechanism is unclear and needs further investigation.

Conclusions
We sequenced the reproductive system of male mice chronically infected with T. gondii, which provides a new direction for research into male sterility caused by Toxoplasma infection. This work provides valuable information and a comprehensive database for future studies of the interaction between T. gondii infection and the male reproductive system. Testis-specific spindle-associated factor that plays a role in spermatogenesis Spz1 1.08852 Expressed specifically in the testis (Sertoli and Leydig cells) and epididymis Play an important role in the regulation of cell proliferation and differentiation during spermatogenesis Spata18 1.29567 In testis, expressed primarily in spermatids Key regulator of mitochondrial quality that mediates the repairing or degradation of unhealthy mitochondria in response to mitochondrial damage