Transcriptome Analysis of Testes and Uterus: Reproductive Dysfunction Induced by Toxoplasma gondii in Mice

Toxoplasma gondii (T. gondii) infection in female mammals during pregnancy can result in poor pregnancy. Similarly, it can result in male reproductive disorders in male mammals. Although the testes and uterus have very different biological makeup, they are still both attacked by T. gondii resulting in reproductive dysfunctions. We hypothesized that there are significant common genes in the testes and uterus that interact with T. gondii. Finding out and studying these genes is vital to understand the infection mechanism of T. gondii and the induced disease pathogenesis. To achieve this goal, we built a mice model of acute infection with T. gondii and the testes and uterus of the mice were sequenced by RNA-Seq. A total of 291 and 679 significantly differently expressed genes (DEGs) were obtained from the testes and the uterus, respectively. In the Gene Ontology (GO) analysis, part of the DEGs in the testes and uterus were related to 35 GO functions. When compared with the KEGG database, seven pathways affecting both the testes and uterus during the course of T. gondii infection were identified. In addition, Toxoplasmosis can significantly affect the expression of Nlrp5 and Insc leading to negative outcomes in the host. On the other hand, the host regulates Gbp7, Gbp2b, and Ifit3 to defend against T. gondii infection.


Introduction
Toxoplasma gondii (T. gondii) is a "clever" obligate intracellular protozoan parasite, affecting most warm-blooded mammals [1]. In fact, about one third of the world population has been infected

Animals and Experiment Set-Up
Eight-weeks old specific-pathogen free (SPF) Kunming (KM) mice, purchased from Guangdong Medical Laboratory Animal Center, were used in this study. Male mice were randomly divided into two groups (n = 5), including infected and control groups. Two female mice were mated with one male mouse. Female mice were judged as the 0th day of pregnancy when vaginal plug occurred. Then, the pregnant mice were randomly divided into two groups (n = 5), including infected and control groups. On the 7th day of female mice pregnancy, they were infected with T. gondii. Two mice groups were injected intraperitoneally with 0.2 mL of physiological saline containing 1 × 10 2 T. gondii tachyzoite RH strain (type I) and labelled as infected groups. In our previous work, we have previously studied the optimal number of RH strain infections per mouse, in order to minimize damage by T. gondii infection [17]. Meanwhile, the control group was injected with an equal dose of physiological saline. On the 12th day of pregnancy in female mice (5th day of T. gondii infection in infected group), all the mice were euthanized by CO 2 . The mouse testes and uterus were quickly harvested and collected for subsequent experiments. One sample was taken from one mouse, and the number of samples was equal to the number of mice. The samples from infected groups were named by gender into MaleTox and FemaleTox groups, and the samples from control groups were also divided into MaleCtrl and FemaleCtrl groups according to gender. Three sets of biological replicates were implemented, for a total of 60 KM mice (male: female = 1:1).
Tachyzoites of the highly virulent RH strain of T. gondii were preserved in our laboratory (Laboratory of Parasitology, College of Veterinary Medicine, South China Agricultural University) and maintained by serial intraperitoneal passage in KM mice.

Evaluation of Sperm Quality Parameters
The epididymis was isolated from the reproductive system of 30 male mice. Afterward, it was cut into small pieces, then 1 mL PBS was added and was placed in a CO 2 incubator at 37 • C for 15 min to make a sperm suspension. Part of the sperm suspension was removed and the total number of spermatozoa was calculated using a hemocytometer. In addition, the premixes of 5% eosin and 10% nigrosine were mixed with the sperm suspension in the proportion of 1:2 and then observed on the optical microscope. Inactive sperm was stained pink, while active sperm was not stained. In total, 500 sperms were counted and the sperm survival rate was calculated.

Ultrastructural Evaluation of Testicles
For this task, two male mice were selected, one from the infected group and the other from the control one. We harvested their testes, and the half of these testes were immersed in 2.5% glutaraldehyde. The testicles were fully filled from one end with the solution. After that, the internal tissue was carefully removed after sliding the blade across albuginea testes, and then divided the size of the tissue into 1 × 1 × 3 mm. The collected samples were fixed overnight in 1% osmic acid, then dehydrated in a series of increasing concentrations of ethanol solution, followed by 100% acetone, and finally, embedded in resin. The ultrathin sections were made and double stained with uranyl acetate and lead citate. After drying, they were observed in transmission electron microscope (TECNAI, Hillsboro, OR, USA). Three parallel samples were prepared.

cDNA Library Construction and Sequencing
Library construction and RNA-Seq were performed according to Illumina's standard procedures. Firstly, the total RNA was extracted from the preserved testis and uterine samples using the RNAiso Plus (TaKaRa, Dalian, China) according to the manufacturer's instruction. Then, RNA quality was assessed using a 1.0% agarose gel electrophoresis and a spectrophotometer (Biotek EPOCH, San Francisco, CA, USA). After the initial evaluation, the integrity of the RNA samples was confirmed by the Agilent Bioanalyzer 2100 and RNA 6000 Nano Lab Chip Kit (Agilent Technologies, Santa Clara, CA, USA). Qualified RNA needs to comply with RNA integrity number (RIN) greater than 7.0. To purify the product, poly (A) mRNA was isolated from approximately 10 µg of qualified RNA using poly-T oligo-attached magnetic beads (Invitrogen, Carlsbad, CA, USA). According to the recommendations of the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA), these mRNA were divided into small fragments by divalent cations under high temperature. The mRNA fragments were reverse transcribed to create the final cDNA library. Afterward, we performed the paired-end sequencing on an Illumina Hiseq 2000/2500 (Illumina, San Diego, CA, USA) following the vendor's recommended protocol.

Sequence Assembly and Acquisition of Clean Reads
The raw data obtained after the previous step were subjected to further filtering to obtain clean reads through removal of the low-quality reads from the raw data. The adaptors in the sequenced reads, which have a ratio of base N ≤ 5% and nucleotides with Q quality score ≥ 20. Q20, Q30, GC%, and the sequence repeatability of clean reads were calculated by Trinity software, as well as the clean reads were stitched.

Identification and Bioinformatics Analysis of DEGs
The Fragments Per Kilobase of exon model per Million mapped reads (FPKM) was used to calculate the abundance of the gene expression. The FPKM of each gene was calculated by using DEGseq software. In this study, the differently expressed genes (DEGs) of females and males were analyzed separately first, and then DEseq was used to compare the expression of various genes in the testes and uterus. The fold change was calculated based on the normalized gene expression level (FPKM-TOX/FPKM-Control). DEGs were identified as genes with "|log2 (fold change)| ≥ 1, p-value < 0.05". To understand these genes, bioinformatics analysis was applied. At first, the biological functions of the DEGs were enriched through the UniProt database, and then annotated the DEGs by the Gene Ontology (GO) (http://www.geneontology.org/) according to their biological process, molecular function, and cellular component. Finally, the pathway of DEGs was predicted through the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/kegg/). KEGG facilitates the understanding of the mechanism of different genes' coordination in the body through formation of an information network.

Verification of mRNA by Quantitative Real-Time PCR Analysis
After bioinformatics analysis, a series of significant DEGs was obtained. In order to verify whether the expression of these genes in the body is consistent with the results of transcriptome analysis, we performed Quantitative Real-Time PCR Analysis (qRT-PCR). The total RNA was extracted from the preserved testis and uterine samples using RNAiso Plus (TaKaRa, Dalian, China). Concentrations and purity were measured by the spectrophotometer (Biotek EPOCH, San Francisco, CA, USA). Approximately 1 ug total RNA was taken based on the measured concentration, then total RNA was synthesized into cDNA using a reverse transcription kit (Vazyme, Nanjing, China). Gene primers for qRT-PCR were designed by Premier 5.0 software ( Table 1). The components were mixed in advance using SYBR qPCR Master Mix (Vazyme, Nanjing, China), and these samples were placed in Roche LifeCycle 96 real-time system at 95 • C for 30 s, followed by 40 cycles of 95 • C for 10 s, and 60 • C for 30 s. Each sample was set to three independent replicates. β-actin was applied as an internal standard for normalizing the expression of gene, the blank control group was used as reference sample, which was set to 1. Additionally, relative quantification of target mRNA was performed by using the cycle threshold (Ct) 2 −∆∆Ct method [18].  ACTGGGAAAGAGATGCACCC  GCCATCATGCACAAACCCAC  147  229900  Gbp7  AGCTCATTGGGTGCAAAAATCC  ACTTACCAGAACCAGGCACTAC  73  14468  Gbp2b  GCCTCGCCTGTGTATCAGGA  GGAAAGTTGCTTCCTGTCTCC  82  15959  Ifit3  TCCTCCAAAAGGCAGCTCAG  TCCCTGTAACGGCACATGAC  133  233752  Insc  GCAGGTAGACTCGGTTCAGC  AGGTCACCTTGCGTGTCTTC  117  110557  H2-Q6  CCCCAGGTCTTATGGTGCTG  TGCCAAGTCAGGGTGATGTC  78  15018  H2-Q7  CCCAGGTCTTATGGTGCTGTC  TCAGCTCCTCCCCATTCAAC  97  15015  H2-Q4  GAAAGCCAAGGGCAATGAGC  CCGCGCTCTGGTTGTAGTAG  74  14964  H2-D1  GACAACAAGGAGTTCGTGCG  CACTCGGAACCACTGCTCTT  144  11461 β-actin AGAGAAGCTGTGCTATGTTGCT GGAACCGCTCGTTGCCAATA 128

Validations of RNA-Seq Results by Western Blotting
Western blot was used to validate RNA-Seq results, by targeting the protein of two differentially expressed genes, NACHT, LRR, and PYD domains-containing protein 5 (Nlrp5) and Interferon-induced protein with tetratricopeptide repeats 3 (Ifit3). The pre-processed protein samples were run on SDS-PAGE (80 V at 30 min, 120 V at 120 min), and then transferred into a polyvinylidene fluoride membrane (Millipore, Darmstadt, Germany) for 120 min at 200 mA. The membrane was blocked in 0.1% Tween 20-PBS containing 5% skimmed milk powder for 2 h and then incubated with β-actin mouse monoclonal antibodies (1:10,000 dilution, Proteintech Group, Chicago, USA) and specific mouse anti-Nlrp5 antibodies (1:1000, Abcam, Massachusetts, USA) for 1 h. Another membrane was incubated with β-actin mouse monoclonal antibodies and anti-Ifit3 antibodies (1:1000, Abcam, Massachusetts, USA). The membranes were probed with goat anti-mouse IgG conjugated with horseradish peroxidase (HRP) (Tiangen Biotech, Beijing, China) at 1:2500 dilution after washing three times. Finally, the membranes were visualized with BeyoECL Plus (Beyotime, Shanghai, China) and the image was taken and analyzed by the Tanon 5200 Western blotting detection system (Tanon, Shanghai, China).

Assessment of Reproductive Damage
Epididymis is the sperm storage room of mammals. Statistics of the total number of sperms in the epididymis are an important indicator of sperm production efficiency. In addition, calculating the sperm viability can reflect the fertility of male mice. Therefore, we chose the total number of sperm in the epididymis and the sperm survival rate as the sperm quality parameters. The sperm quality parameters of each male mouse in the control group (n = 15) and the infected group (n = 15) were counted three times and averaged to reduce the error. As a result, the total number of sperms in the control group was in a range of 6.98 ± 0.71 × 10 7 to 12.62 ± 0.43 × 10 7 , whereas in the infection group, the range was 3.86 ± 0.55 × 10 7 to 11.22 ± 0.33 × 10 7 (p < 0.05). The average sperm survival rates of the control group and the infected group were 71.92% and 46.28%, respectively, which was also statistically significant, p < 0.05 ( Figure 1). The ultramarine structure of testes was observed by the transmission electron microscope (TEM) (Figure 2). The testes of mice were damaged after acute infection with T. gondii. In control group, ultrathin sections of the testes showed that the Sertoli cells were clearly separated from spermatogonia, and the cells were closely arranged. However, in the infected group, the cell borders were blurred, and the BTB composed of tight junctions between the Sertoli cells was damaged. In addition, the tight junction gaps were significantly increased. This phenomenon is a strong evidence of parasites breaking through the BTB, entering the testes from the blood and causing direct damage to it.

Identification of Expressed Transcripts
We divided KM mice into four groups for RNA-Seq analysis including two experimental groups and two control groups. These groups formed from male mice infected with T. gondii (MaleTox), male

Identification of Expressed Transcripts
We divided KM mice into four groups for RNA-Seq analysis including two experimental groups and two control groups. These groups formed from male mice infected with T. gondii (MaleTox), male

Identification of Expressed Transcripts
We divided KM mice into four groups for RNA-Seq analysis including two experimental groups and two control groups. These groups formed from male mice infected with T. gondii (MaleTox), male mice injected with saline (MaleCtrl), pregnant female mice infected with T. gondii (FemaleTox), and pregnant females injected with saline (FemaleCtrl). Three sets of biological replicates are required. Total RNA harvested from the testes and uterus of all mice groups was subjected to high-throughput sequencing on Illumina Hiseq 2000/2500, respectively. 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 2). After filtration process, 31,114 (male groups) and 28,718 (female groups) DEGs were obtained. In order to screen the significant DEGs, it is necessary to conform to "p-value < 0.05, |log 2 (fold change)| ≥ 1". MA-Plots assessed the relationship of log 2 (FPKM) and log 2 (fold change) to realize the display of data distribution (see Supplementary file 1). Compared with the control group, the male mice infected group had 291 significantly DEGs, of which 78 and 213 DEGs were labelled as upregulated and downregulated, respectively. On the other hand, the female mice infected group had 679 significantly DEGs, of which 309 and 370 DEGs were labelled as upregulated and downregulated, respectively. The overall distribution of DEGs can be seen by plotting the volcano maps ( Figure 3), highlighting the genes with fold change ≥ 1 and p-value < 0.05. Among these significant DEGs, we found that 13 genes were co-expressed in the male group and the female group ( Figure 4a). mice injected with saline (MaleCtrl), pregnant female mice infected with T. gondii (FemaleTox), and pregnant females injected with saline (FemaleCtrl). Three sets of biological replicates are required. Total RNA harvested from the testes and uterus of all mice groups was subjected to high-throughput sequencing on Illumina Hiseq 2000/2500, respectively. 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 2). After filtration process, 31,114 (male groups) and 28,718 (female groups) DEGs were obtained. In order to screen the significant DEGs, it is necessary to conform to "p-value < 0.05, |log2 (fold change)| ≥ 1". MA-Plots assessed the relationship of log2 (FPKM) and log2 (fold change) to realize the display of data distribution (see Supplementary file 1). Compared with the control group, the male mice infected group had 291 significantly DEGs, of which 78 and 213 DEGs were labelled as upregulated and downregulated, respectively. On the other hand, the female mice infected group had 679 significantly DEGs, of which 309 and 370 DEGs were labelled as upregulated and downregulated, respectively. The overall distribution of DEGs can be seen by plotting the volcano maps (Figure 3), highlighting the genes with fold change ≥ 1 and p-value < 0.05. Among these significant DEGs, we found that 13 genes were coexpressed in the male group and the female group (Figure 4a).

Gene Ontology of Differentially Expressed Genes
The exploration of biological functions of DEGs from other biological processes, cellular components, and molecular functions were carried out by Gene Ontology (GO) enrichment analysis. The distribution diagram was used to describe the top 10 GO terms of each biological function of the male group and female group ( Figure 5). In the male group, the main categories of biological processes include transport, immune system process, oxidation-reduction process, and cellular response to interferon. As for cellular components, extra cellular exosome, cell surface, and endoplasmic reticulum are closely related to the cells. Furthermore, in molecular functions, GTP binding, oxidoreductase activity, and GTPase activity occupy a very important part. On the other hand, the female group showed that the biological processes are involved mainly in transport, cell differentiation, and cell adhesion. In respect to the cellular components, most of them were involved in membrane, cytoplasm, and extracellular exosome. Whereas, for molecular functions, it includes protein, calcium ion, and metal ion binding, as well as hydrolase and oxidoreductase activity. Among the GO functions, 35 were involved in both the testes and uterus, representing 11, 9, and 15 involved in biological processes, molecular functions, and cellular components, respectively. The identity of the individual genes within the KEGG pathways and GO terms are shown in Supplementary file 2. (a)

Gene Ontology of Differentially Expressed Genes
The exploration of biological functions of DEGs from other biological processes, cellular components, and molecular functions were carried out by Gene Ontology (GO) enrichment analysis. The distribution diagram was used to describe the top 10 GO terms of each biological function of the male group and female group ( Figure 5). In the male group, the main categories of biological processes include transport, immune system process, oxidation-reduction process, and cellular response to interferon. As for cellular components, extra cellular exosome, cell surface, and endoplasmic reticulum are closely related to the cells. Furthermore, in molecular functions, GTP binding, oxidoreductase activity, and GTPase activity occupy a very important part. On the other hand, the female group showed that the biological processes are involved mainly in transport, cell differentiation, and cell adhesion. In respect to the cellular components, most of them were involved in membrane, cytoplasm, and extracellular exosome. Whereas, for molecular functions, it includes protein, calcium ion, and metal ion binding, as well as hydrolase and oxidoreductase activity. Among the GO functions, 35 were involved in both the testes and uterus, representing 11, 9, and 15 involved in biological processes, molecular functions, and cellular components, respectively. The identity of the individual genes within the KEGG pathways and GO terms are shown in Supplementary file 2.

Gene Ontology of Differentially Expressed Genes
The exploration of biological functions of DEGs from other biological processes, cellular components, and molecular functions were carried out by Gene Ontology (GO) enrichment analysis. The distribution diagram was used to describe the top 10 GO terms of each biological function of the male group and female group ( Figure 5). In the male group, the main categories of biological processes include transport, immune system process, oxidation-reduction process, and cellular response to interferon. As for cellular components, extra cellular exosome, cell surface, and endoplasmic reticulum are closely related to the cells. Furthermore, in molecular functions, GTP binding, oxidoreductase activity, and GTPase activity occupy a very important part. On the other hand, the female group showed that the biological processes are involved mainly in transport, cell differentiation, and cell adhesion. In respect to the cellular components, most of them were involved in membrane, cytoplasm, and extracellular exosome. Whereas, for molecular functions, it includes protein, calcium ion, and metal ion binding, as well as hydrolase and oxidoreductase activity. Among the GO functions, 35 were involved in both the testes and uterus, representing 11, 9, and 15 involved in biological processes, molecular functions, and cellular components, respectively. The identity of the individual genes within the KEGG pathways and GO terms are shown in Supplementary file 2.

KEGG Enrichment Analysis
Enrichment analysis of DEGs between the testes and uterus was performed using KEGG database. The male group was significantly involved in 25 pathways, while the female group was involved in 29 pathways. The top 10 enriched pathways in the male group include Metabolic pathways, Phagosome, Antigen processing and presentation, Oxidative phosphorylation, Parkinson's disease, Herpes simplex infection, HTLV-I infection, Graft-versus-host disease, Cell adhesion molecules (CAMs), and Allograft rejection. Whereas, the top 10 enriched pathways in the female group include Metabolic pathways, Biosynthesis of antibiotics, HTLV-I infection, Cell adhesion molecules (CAMs), Focal adhesion, Hippo signaling pathway, cAMP signaling pathway, Wnt signaling pathway, Transcriptional misregulation in cancer, and ECM-receptor interactions ( Figure 6). By comparing these pathways, we found that seven pathways are commonly enriched in male and female groups (Figure 4b), which is briefly described in Table 3. An autoimmune disease of endocrine system characterized by minimal or absent production of insulin by the pancreas.

KEGG Enrichment Analysis
Enrichment analysis of DEGs between the testes and uterus was performed using KEGG database. The male group was significantly involved in 25 pathways, while the female group was involved in 29 pathways. The top 10 enriched pathways in the male group include Metabolic pathways, Phagosome, Antigen processing and presentation, Oxidative phosphorylation, Parkinson's disease, Herpes simplex infection, HTLV-I infection, Graft-versus-host disease, Cell adhesion molecules (CAMs), and Allograft rejection. Whereas, the top 10 enriched pathways in the female group include Metabolic pathways, Biosynthesis of antibiotics, HTLV-I infection, Cell adhesion molecules (CAMs), Focal adhesion, Hippo signaling pathway, cAMP signaling pathway, Wnt signaling pathway, Transcriptional misregulation in cancer, and ECM-receptor interactions ( Figure 6). By comparing these pathways, we found that seven pathways are commonly enriched in male and female groups (Figure 4b), which is briefly described in Table 3.   Table 4 shows the co-expressed genes and fold changes in the male and female groups. We used qRT-PCR to detect the differential expression trends of these genes to verify the reliability of RNA-Seq results. Through the preliminary analysis of the data, the co-DEGs in the testes and uterus were obtained, and then the relevant literature was searched to find out the genes related to pathogen invasion or barrier function among these co-DEGs. There are some genes, although they are co-DEGs between the testes and uterus, but they were not selected for verification. Here, these results are similar to those of RNA-Seq (Figure 7), and PCR amplification efficiency and primer specificity are shown by amplification curve and melting curve, respectively (see Supplementary file 3).   Table 4 shows the co-expressed genes and fold changes in the male and female groups. We used qRT-PCR to detect the differential expression trends of these genes to verify the reliability of RNA-Seq results. Through the preliminary analysis of the data, the co-DEGs in the testes and uterus were obtained, and then the relevant literature was searched to find out the genes related to pathogen invasion or barrier function among these co-DEGs. There are some genes, although they are co-DEGs between the testes and uterus, but they were not selected for verification. Here, these results are similar to those of RNA-Seq (Figure 7), and PCR amplification efficiency and primer specificity are shown by amplification curve and melting curve, respectively (see Supplementary file 3).

Verification Results by Western Blotting
After analyzing the raw data and searching relevant literature, Nlrp5 (UniprotKB: Q9R1M5, 134 kDa) and Ifit3 (UniprotKB: Q64345, 56 kDa) were selected as the research targets to explore their role in T. gondii infection. Therefore, we verified the expression changes of protein by Western blotting (Figure 8). The results were consistent with the sequencing of the transcriptome, where T. gondii infection suppressed Nlrp5 expression and increased expression of Ifit3.

Verification Results by Western Blotting
After analyzing the raw data and searching relevant literature, Nlrp5 (UniprotKB: Q9R1M5, 134 kDa) and Ifit3 (UniprotKB: Q64345, 56 kDa) were selected as the research targets to explore their role in T. gondii infection. Therefore, we verified the expression changes of protein by Western blotting (Figure 8). The results were consistent with the sequencing of the transcriptome, where T. gondii infection suppressed Nlrp5 expression and increased expression of Ifit3.

Verification Results by Western Blotting
After analyzing the raw data and searching relevant literature, Nlrp5 (UniprotKB: Q9R1M5, 134 kDa) and Ifit3 (UniprotKB: Q64345, 56 kDa) were selected as the research targets to explore their role in T. gondii infection. Therefore, we verified the expression changes of protein by Western blotting (Figure 8). The results were consistent with the sequencing of the transcriptome, where T. gondii infection suppressed Nlrp5 expression and increased expression of Ifit3.

Discussion
The placental barrier is a special physiological structure that is formed only during the pregnancy of a female mammal. The zygote of mice developed into blastocysts after 3.5 days, and trophoblasts around the blastocyst proliferated close to the endometrium to become syncytial trophoblast cells [19], and the placental barrier began to form. In early pregnancy, the barrier is weak because the placenta is not fully differentiated, and the uterus is tightly integrated with the embryo. From the second trimester, the structure of the placental barrier is fully completed and the barrier function is well developed [20]. In order to achieve a transcriptomic analysis of the effects of T. gondii infection on the placental barrier, we chose to infect T. gondii in the second trimester of pregnancy (7th day) in mice. Simultaneously, male mice were also infected with T. gondii experimentally. In both females and males, we found that T. gondii go around the reproductive barrier causing direct damage to the mouse reproductive functions.
Using RNA-seq data analysis, we obtained 28,718 and 31,114 DEGs in FemaleTox group and MaleTox group, respectively. We elaborate the common points of the placental barrier and BTB in case of acute infection with T. gondii in respect to DEGs and pathways. In DEGs co-expressed by the uterus and testes, H2-D1, H2-Q4, H2-Q6, and H2-Q7 were all upregulated. These four genes are the major histocompatibility complex (MHC) located on mouse chromosome 17. The main role of these genes is to present foreign antigens to the body's immune system, restricting mutual recognition, and inducing immune responses [21]. In particular, the H2 gene is usually upregulated when foreign microorganisms invade the body [22], indicating that T. gondii has launched an attack on the mouse's reproductive system. These four genes are involved in many pathways such as Allograft rejection, Graft-versus-host disease, Autoimmune thyroid disease, HTLV-I infection, and Cell adhesion molecules (CAMs), and may increase the risk of allogeneic rejection, autoimmune disease, and viral infection.
Nlrp5, also known as Mater, is a conserved gene mainly expressed in the testes and uterus of adult mammals and related to animal reproductive functions [23]. From our experimental results, we found that Nlrp5, as a gene co-expressed by FemaleTox group and MaleTox group, is downregulated compared with the control group, and it is downregulated in the uterus by a factor of more than 304 times (log 2 (fold change) = 8.46). Tong et al. [24] considered that Nlrp5 is a crucial gene for early embryonic development. Knock-out of the Nlrp5 gene in female mice would stagnate the early embryonic development at the two-cell stage, thereby causing infertility. Male mice that lack the Nlrp5 gene are also infertile. Impaired Nlrp5 function that causes reduced fertility is ubiquitous in other mammals, such as monkeys [25], bovine [26], and pigs [27]. However, in our experiments, female mice were infected with the parasites in the second trimester of pregnancy, and embryonic development had passed through the two-cell stage. The results showed that Nlrp5 was still significantly downregulated. The reason may be that Nlrp5 can not only affect the early embryonic development, but also the subsequent fetal development. Interestingly, with the aging of animals, the content of Nlrp5 in the oocytes will also decrease, indicating that the decline in the fertility of older females may be related to the decrease or loss of Nlrp5 [28]. Moreover, some studies have pointed out that Nlrp5 mediates the mitochondrial function in mouse embryos, thus, the lack of Nlrp5 in the mother's body will automatically activate the mitochondrial pool prematurely leading to embryonic mitochondrial damage [29]. It is expected that Nlrp5 may have more complex functions during the maternal pregnancy. Although Nlrp5 is expressed in the testes, this seems to be not related to the reports showing the function of this gene in male animals. Up to now, studies are focusing mainly on the effect of Nlrp5 on the maternal pregnancy. Our experimental results showed that as a pathogen, the invasion process of T. gondii results in downregulation of Nlrp5, especially in the female mice. In this study, we tested the fertility of mice and the results showed that T. gondii can also reduce the fertility of male mice. At present, there is little report on the relationship between T. gondii and Nlrp5. According to the results of the current study, Nlrp5 may be a new target for studying the reproductive damage caused by T. gondii, which breaking through the reproductive system barrier.
Insc is a weight that controls the balance between symmetric cell division (SCD) and asymmetric cell division (ACD). Tissue homeostasis requires an ordered cell division. SCD enhances the proliferation of cells, meanwhile, ACD promotes the tissue differentiation and ensures the cell diversity [30]. In our present study, Insc showed a very different expression in the two different reproductive systems after the invasion of T. gondii. The gene is upregulated five times, (log 2 (flod change) = 2.35) in the MaleTox group and downregulated 11 times (log 2 (flod change) = −3.46) in the FemaleTox group (Table 3). Since Insc interacts with Par3 and LGN, and participates in ACD, Insc has long been considered as the molecular connection between the apically localized complex (Par3/Par6/aPKC) and the microtubule-associated complex (NuMA/LGN/Gai). It promotes spindle-polarity coupling during the mitosis process leading to an unequal cell division [31,32]. Recent studies have found that high expression of Insc will drive the dividing mode of cells to SCD, resulting in an increase in the cell number and volume [33]. In both cases of the upregulation or downregulation of Insc, it may cause the tissue to lose its differentiation function but it can only increase the surface area. It is suggested that the change in the expression of Insc may lead to several pathological reactions in the host and adverse maternal pregnancy. After infection with T. gondii, the opposite change in the expression of Insc in male and female mice may be explained by the different mechanisms of T. gondii infection due to the wide difference in the male and female reproductive systems. It is recommended to combine immunological and pathological experiments for further investigations.
Interferon-inducible Guanylate-binding protein (Gbp, 65-72 kDa) can protect organisms from being infected and facilitate the proinflammatory signal transduction by promoting the host cells' death [34]. There are seven types of Gbps reported in human genome (hGbp1 to hGbp7), whereas 11 types of Gbps in murine genome (mGbp 1-11). During infection of exogenous pathogenic microorganism, many types of murine Gbps [35][36][37] and human Gbp2 and Gbp5 [38,39] can trigger the host cells' pyroptosis by activating inflammasome complex containing caspase-1 or caspase-4. After infection of T. gondii, mGbp2B (also named Gbp1) can activate inflammasome complex signaling cascades, which leads to Nlrp1, Nlrp3, Nlrp4, and pro-caspase-1 degradation resulting in converting the cell death type from pyroptosis to apoptosis [40]. Mouse guanylate-binding protein7 (mGbp7) exhibits high affinity for GTP, nevertheless; it is inferior to the members of small GTPase families such as Ras and Gα. Additionally, it is one of the members with highest affinity in mGbp counterparts [41,42], which enables mGbp7 to hydrolyze GTP more efficiently to promote the host defense capability to some extent. In addition, the vimineous C-terminal tail of mGbp7, composed of 49 amino acid residues, can play a similar role to the C-terminal CaaX motif of other mGbps to function as anchorage by being inserted into intramembrane [42]. Although mGbp7 and other mGbps (mGbp1, -2, -3, -6, and -9) can effectively recruit T. gondii [43], the CT tail on mGbp7 hardly migrates into the parasites' parasitophorous vacuole [42]. This interesting phenomenon suggests that there may be an unknown mechanism of CT tail, which plays an important role in membrane protein localization, membrane anchoring, or interaction.
The interferon-induced proteins with tetratricopeptide repeats (Ifits) are a conserved family of antiviral RNA-binding proteins in vertebrates, and play a significant role in antiviral immune responses [44]. On the one hand, Ifit3, as an intermediate molecule, stabilizes the normal expression of Ifit1 in the host cells helping them to recognize the virus invasion and make a rapid inhibitory response [45]. On the other hand, the tetrapeptide repeat motif of Ifit3 interacts with the N terminus of TBK1 making a bridging protein that connects TBK1 to MAVS on the mitochondria, and is important for protein transport, translational initiation, cell migration and proliferation, antiviral signaling, and virus replication. Recent studies have found that the protection mechanism of Ifit3 for the host is not only against viruses, but is also effective against tumors such as breast cancer [46], and oral squamous cell carcinoma [47], as well as bacteria and parasite infections such as pneumococcus [48], plasmodium [49]. Among the immune proteins induced by Interferon gamma (IFN-γ), Gbps is one of the most abundant proteins. IFN-γ plays a vital role in antimicrobial defense. The lack of IFN-γ receptor signaling in the mammals results in loss of function or mutation of IFN-γ receptor signaling, which can significantly reduce the defense ability of invading foreign microorganisms and make the host extremely susceptible to T. gondii [50]. T. gondii can induce the death of infected cells and promote its own transmission in the host. This situation includes the mechanism of apoptosis. The inhibition of apoptosis varies depending on the type of host and the presence of IFN-γ. Ifits and Gbps proteins belong to IFN-α and IFN-γ stimulated genes (ISGs). Interestingly, researchers are accustomed to using Ifit3 and Gbps as a marker to distinguish between type I interferons (IFN-α and IFN-β) and type II interferons (IFN-γ) [51]. The upregulated expression of marker proteins just shows that the reproductive system can initiate the defense effect of interferon to resist the invasion of T. gondii. Gbps is an interferon-induced protein that can resist the attack of pathogens. Several reports demonstrated the shape of co-expression relationship between Gbps and Ifit3 [52][53][54]. Philippe et al. suggested that ISGs may induce different proteins to cooperate in the control of the production of specific type of IFN [55]. The current mainstream view revealed that the host's autoimmune defense is induced by IFN-γ after T. gondii infection [56]. Type I interferon is mainly involved in the immune system regulation and control of antiviral factors to prevent viral infection. INF-α induces NK cells to produce IFN-γ [57], which is an upstream event of IFN-γ against parasitic infection. The upregulation of Ifit3 may be essential to promote the IFN-α induced expression of IFN-γ in the NK cells, and ultimately enhances the host's functional ability to resist the parasitic infection. On the other hand, research studies investigating the relationship between Ifit3 and T. gondii is scanty. The research on the relationship between Ifit3 and parasites may help us to understand the defense mechanism of the interferon of the host against T. gondii infection.

Conclusions
Results of this study confirm that T. gondii can result in serious reproductive damage in both genders such as decreased sperm count, decreased sperm survival rate, and poor pregnancy in pregnant mice. We also provided transcriptome data of differentially expressed genes in testes and uterus after T. gondii infection and provided data to support subsequent research on parasite invasion to the reproductive system. From the investigated data, we screened a series of genes that can respond to T. gondii infection in the reproductive system of both sexes. They may be the key members to lend a hand to T. gondii to cross through the tissue barrier and cause the subsequent damage to the tissues.