Global Transcriptome Profiling of Multiple Porcine Organs Reveals Toxoplasma gondii-Induced Transcriptional Landscapes

We characterized the porcine tissue transcriptional landscapes that follow Toxoplasma gondii infection. RNAs were isolated from liver, spleen, cerebral cortex, lung, and mesenteric lymph nodes (MLNs) of T. gondii-infected and uninfected (control) pigs at days 6 and 18 postinfection, and were analyzed using next-generation sequencing (RNA-seq). T. gondii altered the expression of 178, 476, 199, 201, and 362 transcripts at 6 dpi and 217, 223, 347, 119, and 161 at 18 dpi in the infected brain, liver, lung, MLNs and spleen, respectively. The differentially expressed transcripts (DETs) were grouped into five expression patterns and 10 sub-clusters. Gene Ontology enrichment and pathway analysis revealed that immune-related genes dominated the overall transcriptomic signature and that metabolic processes, such as steroid biosynthesis, and metabolism of lipid and carboxylic acid, were downregulated in infected tissues. Co-expression network analysis identified transcriptional modules associated with host immune response to infection. These findings not only show how T. gondii infection alters porcine transcriptome in a tissue-specific manner, but also offer a gateway for testing new hypotheses regarding human response to T. gondii infection.


INTRODUCTION
Toxoplasma gondii is an obligate intracellular protozoan parasite, which infects nearly one-third of the world human population and all warm-blooded vertebrates (1,2). There are three infectious stages of this parasite: tachyzoites, bradyzoites-containing tissue cysts, and sporozoites-containing oocysts (the product of sexual reproduction in the intestine of the feline definitive host). Humans acquire infection through ingestion of raw or undercooked meat, such as pork or lamb containing cysts (3,4). Also, infection can be acquired by ingestion of food contaminated with oocysts or by exposure to soil containing oocysts (5). Despite significant research progress, our understanding of immune-related genes which are substantially involved in the pathogenesis of human toxoplasmosis remains limited. A considerable mass of research has been performed using cell culture models, which improved understanding of the pathogenesis of T. gondii infection. However, in vitro models cannot fully recapitulate in vivo processes.
Animal models can reduce the deficiencies that are inherent to in vitro models. Mice (Mus musculus) are the most widely model used to study T. gondii pathophysiology due to low cost and the availability of specific reagents (6)(7)(8). However, mice are less suitable as a model for understanding the transcriptional landscape of other mammalian species that have different transcriptional and genetic backgrounds, such as pig and humans (9). In contrast, transcriptomics and genetic technologies have shown pigs to be genetically and mechanistically relevant to study human conditions (3,10,11). Importantly, the common attributes of T. gondii infection, such as severity of infection, transplacental transmission, and interferon-gammarelated antiparasitic effector mechanisms are more similar in pigs and humans compared to the same aspects of disease in mice (4,(12)(13)(14)(15)(16)(17)(18). These facts suggest that domestic pigs (Sus scrofa domesticus) are more relevant model to the study of the pathophysiology of human toxoplasmosis.
Previous studies have provided beneficial but limited insights into the porcine response to T. gondii through using methods such as high-throughput sequencing to identify the mRNA and miRNA profiles (19)(20)(21)(22), and microarrays (23). Recent research showed that T. gondii loads vary across different porcine tissues, with high parasite loads detected in the heart and lungs during acute infection, and in the heart and brain during chronic infection, regardless of the strain of the parasite (24). Therefore, genome-wide comprehensive analysis of the differential responses of pig tissues to T. gondii infection is required to elucidate why some porcine tissues vary greatly in their response to T. gondii infection.
In this study, the transcriptomic response of five different porcine tissues [brain, lung, liver, spleen, and mesenteric lymph nodes (MLNs)] to experimental T. gondii infection was examined using RNA-sequencing (RNA-seq). Our analysis revealed the global transcriptomic changes in relation to T. gondii infection at 6 and 18 days after infection. We identified hundreds of differentially expressed transcripts (DETs), infection-specific expression patterns and porcine genes that correlated with T. gondii load in infected pig tissues.

Ethics and Biosafety Statement
The study design was reviewed and approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute (LVRI), Chinese Academy of Agricultural Sciences (CAAS). The procedures involving animals were carried out in accordance with the Animal Ethics Procedures and Guidelines of the People's Republic of China. Animals were monitored every day for the development of clinical signs of toxoplasmosis. All efforts were made to minimize suffering and to reduce the number of pigs used in the experiment. The potentially infectious clinical and laboratory waste, such as the remaining pig tissues and T. gondii oocysts, were decontaminated by autoclaving prior to disposal in accordance with the local institutional health and biosafety policy on the disposal of hazardous waste.

Animals and Parasite Challenge
Twenty-four, 14-week-old, specific-pathogen-free (SPF), outbred female white pigs were purchased from Beijing Center for SPF Swine Breeding and Management. To confirm the T. gondiifree status of pigs before being used in the experiment, pig serum was tested using modified agglutination test (MAT) as described previously (25). The 24 T. gondii-seronegative pigs were randomly assigned to eight groups (3 pigs/group), which were housed in separate units. The experimental groups included two control groups (6C_1 and 6C_2) at 6 days post infection (dpi), two control groups at 18 dpi (18C_1 and 18C_2), two infected groups at 6 dpi (6I_1 and 6I_2), and two infected groups at 18 dpi (18I_1 and 18I_2). In the infected groups, each pig was infected orally with 1,000 oocysts of T. gondii PYS strain (genotype ToxoDB#9) in 5 ml sterile Phosphate Buffered Saline (PBS, pH 7.4). Pigs in the uninfected (control) group received 5 ml sterile PBS without any oocysts. Pigs from infected and control groups were euthanized at day 6 and day 18 post infection. These two time points post infection were chosen because pig requires 6 days to develop IgM antibodies (indicative of acute infection), whereas IgG antibodies which mark chronic infection develop after 18 dpi (26). Tissue samples were collected from cerebral cortex (thereafter called brain), liver, spleen, lung and MLNs, and were stored separately at −80 • C. A total of 200 mg collected from several sites of each organ were used for RNA extraction. Confirmation of T. gondii infection in all collected tissues was performed by PCR, as described previously (27).

Determination of Normalized Parasite
Load in Pig Tissues Using Quantitative PCR (qPCR) DNA of pig tissues was extracted using a TIANamp Genomic DNA Kit (Tiangen Biotech, Beijing, China) according to the manufacturer's recommendations. T. gondii B1 gene was used to determine normalized parasite load in pig tissues, and porcine gene coding for 18S rRNA was used to normalize the T. gondii B1 DNA to host DNA. Oligonucleotides for amplification of the porcine housekeeping gene 18S rRNA were: 18S-pigF (GCCTGCTGCCTTCCTTG) and 18S-pigR (ATGGTAGTCGCCGTGCC), with an expected product size of ∼109 bp. The primers used for detection of T. gondii B1 were: B1F (TGCATAGGTTGCAGTCACTG) and B1R (TCTTTAAAGCGTTCGTGGTC) with an expected product size of ∼131 bp. The samples with an exponential-amplification curve crossing the threshold were deemed positive for T. gondii, whereas samples with no amplification curve for T. gondii B1, but amplification of the 18S rRNA gene were considered negative for T. gondii. The 2 − CT method [the method can also be displayed as 2 −(CT value of target gene in tested group − CT value of housekeeping gene in tested group) /2 −(CT value of target gene in control group − CT value of house keeping gene in control group) ] is generally used for calculation of the fold-change of the target genes in infected relative to control samples (28). However, because T. gondii gene cannot be detected in the tissue that free of T. gondii, we used the  The distribution of DETs across the chromosomes. From outer to inner circle represent chromosomes, DET bar plots of brain6, brain18, liver6, liver18, lung6, lung18, mesenteric lymph nodes6, mesenteric lymph nodes18, spleen6, and spleen18, respectively. (C) Clustering of DETs encoded by X chromosome. Red, green, and gray colors represent upregulated, downregulated and missing (undetected) data, respectively. cycle threshold (CT) value of target and housekeeping gene to calculate the relative abundance of T. gondii B1 gene normalized to pig 18S rRNA gene in each infected tissue using the equation 2 − CT . − CT = -(CT for the targeted T. gondii B1 gene -CT for the pig 18S rRNA gene). qPCR was performed in a Rotor-Gene Q system (QIAGEN, Hilden, Germany) using GoTap R qPCR Master Mix (Promega, Madison, WI, USA). The cycling conditions were 95 • C for 5 min followed by 50 cycles of 95 • C for 10 s, 60 • C for 10 s, 72 • C for 15 s; the temperatures of the melt curve analysis ranged from 72 to 95 • C to ensure the specificity of the amplification products.

RNA Extraction and Transcriptome Sequencing (RNA-Seq) Analysis
Total RNA of each collected tissue sample was extracted separately using TRIzol Reagent (Invitrogen China Ltd, Beijing, China) according to the manufacturer's instructions. The residual genomic DNA in the isolated RNA was removed using 20 units of RNase-Free DNase (Ambion, Shanghai, China). The integrity and quantity of RNA samples were tested with Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Thermo Scientific TM Nanodrop 2000 (Wilmington, DE, USA), respectively. The RNA-seq analysis was based on two biological replicates per experimental group, and each biological replicate involved pooled RNA samples from three different pigs within each group. Although sample pooling design masks the individual response of pigs, it has been considered as a cost-efficient approach (29)(30)(31)(32). Approximately one microgram of total RNA per each pooled sample was used as an input for the construction of mRNA library using IlluminaTruSeq TM RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). The transcriptome libraries were sequenced using IlluminaHiSeq TM 2000 according to the manufacturer's instructions. Adaptors and low quality sequencing reads were filtered using a quality cutoff score of Q20, which is widely used for quality control analysis (33)(34)(35)(36)(37). The resulting clean reads were mapped against the pig (Sus scrofa domesticus) reference genome (Sscrofa10.2) using SOAP aligner/SOAP2 software and genomic annotation data file (ref_Sscrofa10.2_top_level.gff3). The level of expression was calculated in units of reads per kilobase per million mapped reads (RPKM) (38). Expression analysis was performed using NOIseq R package (31). Transcripts with |log2 fold-change (FC)| ≥1 and significant value >0.8 were considered differentially expressed, as per the recommendations of the NOIseq. RNA isolation, library construction, RNA sequencing, and computational analysis were performed by BGI-Shenzhen, China.

Validation of RNA-Seq Data
Thirty-two genes identified by RNA-seq analysis, across all experimental groups, were selected for validation by quantitative reverse transcriptase PCR (qRT-PCR). RNA preparations were those used for RNA-seq experiments at the corresponding time points. cDNA was synthesized from total RNA using the PrimeScript TM II 1st Strand cDNA Synthesis Kit (Takara, Dalian, China) according to the manufacturer's instructions. All qRT-PCR experiments were performed in triplicate, with the housekeeping gene GAPDH as a control. The qRT-PCR oligonucleotide primers used in this study are described in Supplementary Table 1. qRT-PCR was performed in Rotor-Gene Q system (QIAGEN, Hilden, Germany) and using GoTap R qPCR Master Mix (Promega, Madison, WI, USA). qRT-PCR cycling conditions were as follows: 95 • C for 5 min followed by 50 cycles of 95 • C for 10 s, 60 • C for 10 s, and 72 • C for 15 s; the temperatures of the melt curve analysis ranged from 72 to 95 • C to ensure the specificity of qRT-PCR products. The 2 − CT relative expression method was used to calculate gene expression.

Gene Ontology (GO) Enrichment and KEGG Analysis
GOseq package (v1.22) in R (www.r-project.org) was used for Gene Ontology (GO) enrichment analysis, such as biological process (BP), cell component (CC), molecular function (MF). Pathway analysis was performed using Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Significantly enriched GO terms or pathways were identified using hypergeometric test followed by FDR correction method (39). The FDR corrected P < 0.05 was used as cutoff for the significantly enrichment GO terms or pathways. All differentially expressed transcripts (DETs) were clustered with log2 fold-change of DETs using Gene Cluster 3.0 and Euclidean distance. We used Upset (40) to visualize intersecting sets in order to identify the unique and common DETs across 10 tissue subsets.

Coexpression Network and Correlation Analysis
We performed coexpression analysis, which has been widely applied to identify genes involved in host-parasite interaction (41)(42)(43). The weighted gene correlation network analysis (WGCNA) R package (44) was used to establish a correlation matrix between pig mRNA expression and normalized T. gondii load in each infected tissue. WGCNA was performed as the network construction and module detection protocol in the WGCNA R package (https://horvath.genetics.ucla.edu/ html/CoexpressionNetwork/Rpackages/WGCNA/index.html). RPKMs of all transcripts were used as input data. We have chosen a soft-thresholding power of 14 because this was the lowest power needed to exceed a scale-free topology fit index of 0.75. Days post infection, infection status and normalized T. gondii load were used as input traits in the module-trait association analysis. The cluster of highly interconnected genes that share a similar expression pattern was considered as a coexpression module. Multidimensional scaling (MDS) plot was constructed to visualize pairwise relationships specified by a dissimilarity matrix, indicating dissimilarity/similarity based on gene expression data. The hub gene of coexpression module was identified based on a high module membership value or connectivity (i.e., the sum of connection strengths with the other module genes). Generally, the hub genes of specific module are located at the finger tip of MDS plot. For further testing of the predictive performance between the host gene expression and T. gondii load, pROC package was used to perform receiver operating characteristic (ROC) curve analysis, and to calculate the area under the ROC curve (AUC), a performance metric, which is generally used for the identification of potential biomarkers. The gene in coexpression module that was significantly correlated with T. gondii load, and had significant P < 0.05, RPKM > 1, AUC > 0.6, was considered as a host gene significantly correlated with T. gondii load (HGSCTG). HGSCTG genomic hotspots were defined on the basis of >5 HGSCTG per 10 Mb genomic region. Gene regulatory networks were reconstructed using the coexpression data and TRRUST database, and were visualized with Cytoscape software. The regulatory transcription factors were identified as the ones that co-express with their target genes. Finally, the relationships between enzymes and substrates were analyzed using the online DrugBank database (https://www.drugbank.ca/).

Accession Numbers
The RNA-Seq datasets described in this study have been deposited in NCBI Short Read Archive database (https://www. ncbi.nlm.nih.gov/sra) under accession numbers SRR6203124 to SRR6203163.

Confirmation of T. gondii Infection and Normalized Parasite Load in Pig Tissues
At 6 dpi, all pigs in the infected groups exhibited clinical signs, such as fever and inappetence, whereas pigs in the control groups remained clinically healthy. The brain, liver, spleen, lung, and MLNs of infected pigs were all PCRpositive, whereas no T. gondii B1 gene amplification was achieved in corresponding tissues from uninfected (control) pigs. Infected lungs showed the highest T. gondii load (Supplementary Table 2).

Transcriptomic Features and Validation of RNA-Seq Results
The RNA integrity numbers (RINs) of all RNA templates were >8.0. Also, 99% of the reads showed high quality values > Q20, and 90% of the clean reads were up to Q30 (Supplementary Figures 1-3). More than 62 million clean reads were obtained from each tissue sample and more than 32,000 transcripts were detected (Supplementary Table 3  Frontiers in Immunology | www.frontiersin.org FIGURE 4 | Differentially expressed components of neuron signaling pathways in T. gondii infected brain. Pink color represents neuron signaling pathway and cyan color represents the differentially expressed transcripts in brain at 6 or 18 dpi. The bars underneath the nodes represent the change in the expression level (left and right bars represent expressional changes at 6 and 18 dpi, respectively). Upregulated and downregulated transcripts are represented by red and green colors, respectively. The log2 fold-change of transcripts are described together with the bars. brain, liver, lung, MLNs and spleen, respectively; whereas at 18 dpi, 217, 223, 347, 119, and 161 DETs were found in the infected brain, liver, lung, MLNs and spleen, respectively ( Figure 1A). The global Pearson correlation coefficient between qRT-PCR results and RNA-seq results was 0.715 (Figure 1B), suggesting a good agreement between results obtained by the two methods, supporting the validity of the transcriptomic RNAseq data. After intersecting the differentially expressed transcript sets across tissues, there was not any DET shared across all tissues as shown in the vertical bars and the connected black circles below the histogram (Figure 1C), suggesting a lack of commonly DETs shared across the analyzed tissues. According to Euclidean distance, as shown in Figure 2A, all DETs were clustered into five distinct expression patterns, including (pattern 1) low expression in most infected tissues (downregulated in ≥ 6 tissue samples at 6 and 18 dpi), (pattern 2) high expression in infected liver only, (pattern 3) high expression at 18 dpi, (pattern 4) high expression at 6 dpi, and (pattern 5) high expression in most infected tissues (upregulated in ≥ 6 tissue samples at 6 and 18 dpi). The distribution of DETs across chromosomes is shown in Figure 2B. Three chromosomes encoded most of the DETs: chromosome 2 (110 DETs), chromosome 6 (108 DETs), and chromosome 7 (114 DETs). No DETs were found in the mitochondrial DNA or chromosome Y. However, 37 DETs were encoded by chromosome X and were clustered into two expression patterns (Figure 2C).

GO Enrichment and KEGG Analysis
We analyzed the functional enrichment and significant pathways associated with the DETs clustered in the five expression patterns. The downregulated transcript cluster in most infected tissues (pattern 1) was significantly enriched for GO terms involved in metabolic or tissue development processes, such as proteinaceous extracellular matrix, lipid metabolic process, animal organ development, PPAR signaling pathway, and metabolism of xenobiotics by cytochrome P450 (Supplementary Table 4). In the upregulated transcript clusters in infected liver (pattern 2) and most infected tissues (pattern 5), most of the enriched transcripts were related to GO terms or pathways involved in immune response. These included cytokine receptor binding, regulation of interleukin (IL)-12 production, Jak-STAT signaling pathway, nuclear factor kappa B (NF-κB) signaling pathway, Analysis is based on the pig reference genome (Sscrofa10.2) and the annotation file (Sscrofa10.2_top_level.gff3).
lymphocyte chemotaxis, IL-8 secretion, and cytokine secretion (Supplementary Tables 5, 6). The cluster containing upregulated transcripts at 6 dpi (pattern 4) was not enriched in any GO term. However, most significantly enriched GO terms in the transcript's cluster with high expression pattern at 18 dpi (pattern 3) were related to immune response, such as antigen processing and presentation pathway, lytic vacuole, response to cytokine, lymphocyte-mediated immunity, and chemokinemediated signaling pathway (Supplementary Table 7). Only five transcripts encoded by X chromosome were involved in lipid metabolic processes (Figure 2C). We also performed GO enrichment and KEGG analysis of all DETs. According to FDR corrected P-value, the top 30 significantly enriched GO terms and pathways are shown in Figure 3. Most of the transcripts that were significantly enriched in immune-related GO terms or pathways belonged to clusters 7-13, which had upregulated expression patterns (Figure 2A). However, most of the transcripts that were significantly enriched in metabolic related GO terms or pathways belonged to clusters 1-6, which had low expression patterns (Figure 2A). In the infected brains, 19 neural signaling pathways were regulated by 15 DETs, including Adora2a, P2ry13, Grin2c, Gng11, Gng13, Ppp1r1b, Pdyn, Gnal, Rac2, Atp5g1, Ndufc2, Slc5a7, Ube2l6, Tnf, and Mme (Figure 4).

Coexpression and Regulatory Network Analysis
The soft-threshold power, that exceeded a scale-free topology fit index of 0.75 for each network was 14 (Supplementary Figure 4). Therefore, number 14 was chosen as the soft power threshold for constructing WGCNA. As shown in Figure 5A, 18 modules, including gray module, denoting unassigned transcripts, were found. The module-trait association shows that MElightyellow was significantly correlated with infection status and T. gondii load (Figure 5A). The global correlation coefficient between MElightyellow and T. gondii load was 0.51 with a P-value of 7e-04. Details of the transcripts and module relationships are summarized in Supplementary Table 8. Multidimensional scaling (MDS) plot shows that the transcripts of each module were successfully categorized by coexpression analysis (Figure 5B). Most of the genes in MElightyellow module were upregulated in most infected tissues, and significantly enriched in immune-related biological processes (Figures 6A,B). The genomic locations of genes in the MElightyellow module are shown in Figure 6C.
We identified 152 genes (165 transcripts in total) as HGSCTG that were significantly correlated with T. gondii loads. The details of HGSCTG are listed in Supplementary Table 9. Chromosome 2 encodes most of the MElightyellow module genes (29 genes) and these genes were distributed on 3 spots, including head spot (chromosome 2: 0Mb-13Mb), which encodes 10 genes (8 of these were HGSCTG), middle spot (chromosome 2: 53-84 Mb), which encodes 10 genes (8 of these were HGSCTG), and end spot (chromosome 2: 124-162 Mb), which encodes 9 genes (4 of these were HGSCTG).  Table 1. Results of the ROC analysis of the top 30 hub genes of MElightyellow module are shown in Figure 6D. All the areas under ROC curve (AUC) were >0.7 and the majority of them were >0.9, indicating a high correlation between the top 30 hub gene expression and T. gondii parasite DNA load. Functions of the hub genes are described in Table 2.

Cytokine Expression
We detected 38 cytokines and 21 cytokine receptor-related transcripts that were differentially expressed, including 18 differentially expressed chemokines and seven differentially expressed chemokine receptors. Most of these were upregulated in the infected tissues (Figure 8A), including four HGSCTG (Cxcl10, Il27, Il15, and Il15ra). In infected tissues, upregulation of chemokines and chemokine receptors increases the chemotaxis of 20 immune cells, such as DC, NK, macrophage, and T cells (Figure 8B).

Comparative Toxicogenomic Analysis
We found that 45 DETs were involved in xenobiotics or drug metabolism ( Figure 8C). Most of these were downregulated, especially in the liver at 6 dpi. We also found that the DETs encode enzymes that metabolize 330 substances or drugs, such as ethanol, acetaminophen, ketoconazole, phenobarbital, and benzyl alcohol. The relationship among 330 xenobiotic substances and DETs related to metabolism are shown (Supplementary Table 10).

DISCUSSION
We compared the transcriptomes of T. gondii-infected and uninfected pigs using RNA-seq approach. Hundreds of transcripts were differentially expressed in the porcine brain, liver, lung, spleen, and MLNs at 6 and 18 dpi (Figure 1). These DETs were distributed on all the pig chromosomes, but not the mitochondrial genome and chromosome Y ( Figure 2B).
We tested whether transcriptomic changes overlap across porcine tissues. As shown in Figure 1C, there was no common DET in the 10 tissue groups. Next, we characterized the transcriptomic changes in infected tissues using gene clustering analysis, which showed that all DETs are clustered into five different expressional patterns (Figure 2A). Functional enrichment analysis of the downregulated transcripts in most infected tissues revealed that downregulation of metabolismrelated and tissue development-related transcripts is prominent in infected tissues (Supplementary Table 4 and Figure 3). This finding may have clinical relevance. During pregnancy, motherto-fetus transmission of T. gondii can occur, resulting in abortion, stillbirth, or congenital malformation (71). It remains to be determined if the downregulation of these transcripts observed in pigs can also occur in the fetus if they become congenitally infected. Mindful of the fact that successful pregnancy requires delicate immune balance to protect the fetus, the deleterious effects of T. gondii induced-inflammatory response mediated by increased expression of cytokines and cytokine receptor-related transcripts (Figure 8A and Supplementary Table 6) may cause undesirable health consequences in the fetus.
Forty-three genes involved in lipid metabolic processes, including 26 genes belonging to clusters 1-6 ( Figure 3A) showed lower expression in infected tissues (Figure 2A). Five of these were encoded by chromosome X (Figure 2C). These results suggest that, during T. gondii infection, chromosome X contributes ∼ one-eighth of the downregulated genes involved in lipid metabolism. The downregulation of metabolic terms or pathways is consistent with our previous proteomic and transcriptomic investigations in mice (72)(73)(74), suggesting that downregulation of metabolic processes may also occur in other mammalian hosts.
T. gondii influences mouse behavior via altering glutamate transporter Slc1a2 (also known as GLT-1) (75) and GABA signaling pathway (76). In our study, Slc1a2 was not differentially expressed in the brain. However, three genes (Gng11, Gng13, and Grin2c) involved in the signaling mechanism of glutamatergic and GABAergic synapse were differentially expressed. This indicates that T. gondii may alter pig behavior by interfering with glutamatergic and GABAergic synapse pathways via altering the expression of Gng11, Gng13, and Grin2c. We also found 19 neural synapse signaling pathways altered in infected brains (Figure 4), such as olfactory transduction pathway, which may alter the sense of smell of infected pigs. Chronically infected rodents exhibit behavioral changes, such as loss of aversion and even attraction to cat odors (77). In humans, infection with T. gondii has been associated with behavioral abnormalities (78) and increased risk of developing psychiatric disorders (79).
We further investigated which host genes are significantly correlated with T. gondii load (HGSCTG) using WGCNA analysis, which identified 18 coexpression modules (Figure 5 and Supplementary Table 8). By relating modules to sample traits (days after infection, infection status, and T. gondii load in tissues), MElightyellow module was significantly correlated with parasite load in infected tissues ( Figure 5A). Most of the transcripts in MElightyellow module were upregulated in infected tissues ( Figure 6B). As shown in Figure 6A, genes in MElightyellow module were significantly enriched in immune response and infection-related terms or pathways. By combining WGCNA coexpression and ROC curve analyses, we identified 152 HGSCTG (Supplementary Table 9), including three HGSCTG cytokines (Cxcl10, Il27, Il15) and one cytokine receptor (Il15ra). These four genes seem to be important for mice survival during T. gondii infection ( Table 2).
Three HGSCTG genomic hotspots ( Figure 6C) encoding 24 HGSCTG ( Table 1) were also identified. Most of these 24 genes were involved in immune processes, and some have anti-T. gondii activity, such as Gbp1, Gbp2, Gbp7, Batf2, and Tap1. As shown in Figure 6D, the AUC of the top 30 hub genes in the MElightyellow module was >0.7 and for most of these genes, the AUC was >0.9, indicating a strong correlation between the top 30 hub gene expression and T. gondii DNA load. This result shows a synergy between the results obtained by ROC and WGCNA analysis, suggesting that the identified hub genes (CD180, STX11, FAM26F, TFEC, ETV7, LOC100738612, and LOC100520491) are HGSCTG.
We detected 31 differentially expressed TFs, grouped into upregulated and downregulated clusters (Figure 7A). In these   differentially expressed TFs, Batf2, Irf7, Irf1, Tfec, and Stat1 were HGSCTG. We also found hundreds of gene targets to 15 differentially expressed TFs in the TRRUST database ( Figure 7B). The coexpression analysis revealed that 12, 7, 1, and 4 targeted genes shared similar expression pattern with Stat1, Irf1, Irf8, and Stat3, respectively (Figure 7). The Irf1 and Stat1 contribute to T. gondii control via regulating the expression of factors essential for host resistance to infection, such as TAP complex (80), Cxcl10, Ciita (81), Il27, and Jak2 (82). We also identified 38 cytokines and 21 cytokine receptor-related transcripts that were differentially expressed. Most of these were upregulated in infected tissues ( Figure 8A). In agreement with others (83), upregulation of these genes can increase chemotaxis of 20 immune cells, including DCs, NK cells, macrophages, and T cells in most infected tissues ( Figure 8B). These immune cells, which play important roles in T. gondii control (6), can be regulated by CXCl9, CXCl10, and CXCR3 signaling pathways (Figure 8B), contributing to the pig immune response to T. gondii infection. As shown in Figure 8C, 45 DETs were involved in xenobiotics or drug metabolism, most of these were downregulated in infected tissues, especially in the liver at 6 dpi. According to the DrugBank database, 330 xenobiotics or drugs were found to be metabolized by enzymes coded by these DETs (Supplementary Table 10). Acetaminophen is used to control fever, however, it can induce adverse events, such as acute liver failure (84). In our study, four downregulated genes (Cyp1a2, Cyp2e1, Cyp1a1, and Sult2a1) were involved in the pharmacokinetic of acetaminophen. The downregulation of metabolic transcripts are alleviated in infected livers at 18 dpi (Figure 3), indicating that downregulation of metabolic processes in infected liver varies by stage of infection and that acute T. gondii infection causes more inhibition of the hepatic metabolic processes. This result is consistent with previous work showing downregulation of genes involved in liver metabolism following T. gondii infection (72,73).
Downregulation of xenobiotics or drug metabolism in liver is related to inflammatory response (85). Although transcripts in the cluster showing higher expression in most infected tissues were significantly enriched in immune-related terms and pathways (Supplementary Table 5), transcripts in the high expression pattern in the liver cluster (Figure 2A) were also significantly enriched in terms related to inflammatory processes (Supplementary Table 6). This suggests that liver exhibited more inflammatory response than other tissues. Likewise, transcripts in the cluster with high expression pattern at 18 dpi were significantly enriched in GO terms related to immune responses (Supplementary Table 7). The ability of T. gondii to cause severe hepatic pathologies has been demonstrated (86,87), and it is possible that the prominent inflammatory response observed in our study contributes to the pathologies observed in infected livers.

CONCLUSION
RNA-seq was used to determine the global changes in the porcine transcriptome subsequent to T. gondii infection at 6 and 18 days post infection. Hundreds of DETs exhibited differential expression profiles in infected tissues and were clustered into five expression patterns. Infection induced downregulation of various metabolic processes in most infected tissues, especially in the liver during acute infection. The WGCNA analysis showed that, T. gondii infection causes differential expression of transcription factors, such as Irf1, Irf8, Stat1, and Stat3. We also identified 45 DETs encoding detoxifying enzymes involved in the metabolism of 300 xenobiotics or drugs. These data improve our understanding of the molecular changes that occur during T. gondii infection in pigs. Although results obtained in pigs may not be readily transferrable to humans, given the physiological and immunological similarities between pigs and humans, our findings may facilitate the understanding of how humans might respond to T. gondii infection.

ETHICS STATEMENT
The study design was reviewed and approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute (LVRI), Chinese Academy of Agricultural Sciences (CAAS). The procedures involving animals were carried out in accordance with the Animal Ethics Procedures and Guidelines of the People's Republic of China. Animals were monitored every day for the development of clinical signs of toxoplasmosis. All efforts were made to minimize suffering and to reduce the number of pigs used in the experiment.

AUTHOR CONTRIBUTIONS
X-QZ, HE, and J-JH conceived and designed the study and critically revised the manuscript. J-JH performed the experiment, analyzed the transcriptomic data, and drafted the manuscript. JM and HE helped in study design, implementation, data analysis, and manuscript revision. J-LW, F-KZ, J-XL, B-TZ, and Z-XW helped in the study implementation. All authors have read and approved the final manuscript.