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–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-gamma-related antiparasitic effector mechanisms are more similar in pigs and humans compared to the same aspects of disease in mice (4, 12–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–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.

Materials and Methods

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. gondii-free 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 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^® 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™ 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–32). Approximately one microgram of total RNA per each pooled sample was used as an input for the construction of mRNA library using IlluminaTruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). The transcriptome libraries were sequenced using IlluminaHiSeq™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–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™ 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^® 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–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.

Results

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 PCR-positive, 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). At 6 dpi, 178, 476, 199, 200, and 362 DETs were detected in infected 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 RNA-seq 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).

FIGURE 1

Figure 1. (A) The number of differentially expressed transcripts (DETs) across different infected porcine tissues at 6 and 18 days postinfection. (B) Validation of RNA-seq results by qRT-PCR. Plot of gene expression (fold change) determined by the RNA-seq (X-axis) and qRT-PCR (Y-axis) of 32 selected genes (Pearson's correlation, R² = 0.715, P < 0.01). The fold-change of expression was expressed as log2 values. (C) UpSet plot showing the sets of differentially expressed transcripts from 10 different tissue samples, including the quantitative analysis of aggregate intersections between tissues. The vertical bars show the number of intersecting transcripts between tissues, denoted by the connected black circles below the histogram. The horizontal bars show the transcript set size.

FIGURE 2

Figure 2. Expression patterns of differentially expressed transcripts (DETs). (A) Cluster analysis of all DETs based on the average values of expression of two replicates. (B) 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.

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, 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 chemokine-mediated 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).

FIGURE 3

Figure 3. Stacked bar plots showing the results of GO enrichment (A) and KEGG analysis (B) of genes in the 13 clusters identified in Figure 2A.

FIGURE 4

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.

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.

FIGURE 5

Figure 5. (A) Module and traits relationships. The color-coded legend shows the correlation degree, with positive and negative numbers indicating positive and negative relationship, respectively (B). The classical multidimensional scaling plot of all transcripts.

FIGURE 6

Figure 6. Functional annotation and transcript expressional cluster of MElightyellow module. (A) GO and KEGG enrichment of transcripts in MElightyellow module. The top 10 GO terms are listed according to the enrichment P-value. (B) Expressional clustering of transcripts in MElightyellow module. Upregulated, downregulated, and missing transcripts are denoted by red, green, and gray colors, respectively. (C) Distribution of the host genes significantly correlated with T. gondii load (HGSCTG) on the pig genome. Outer circle represents the chromosome, scatter in inner circle represents the distribution of HGSCTG, the numbers adjacent to inner circle represent the number of HGSCTG that are encoded by the chromosome, red scatter shows HGSCTG distributed on the hotspot sites. (D) ROC tests between the top 30 hub gene expressions and T. gondii parasite DNA load.

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). Chromosome 4: 139–144 Mb encodes 6 HGSCTG and chromosome 7: 22–38 Mb encodes 10 HGSCTG. Therefore, chromosome 2: 0.3–13 Mb, chromosome 4: 139–144 Mb, and chromosome 7: 22–38 Mb were identified as HGSCTG genomic hotspots. The details of HGSCTG locations in the hotspots are shown in 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.

TABLE 1

Table 1. The genomic location of the host genes significantly correlated with T. gondii load (HGSCTG) in the HGSCTG hotspots.

TABLE 2

Table 2. Description of mouse or human orthologs of the top 30 hub genes of MElightyellow module.

Differentially Expressed Transcription Factors (TFs) and Their Regulatory Networks

We detected 31 TFs in the infected brain, liver, lung, spleen, and MLN. These 31 TFs were grouped into 2 clusters (Figure 7A). The upregulated TF cluster included Fos, Hopx, Zbtb16, Mycl, Junb, Litaf, Stat3, Tead4, Foxs1, Batf, IRF8, BATF2, IRF7, IRF1, Tfec, and Stat1. The downregulated TF cluster included Msc, Id4, Dlx5, Dlx1, Dlx2, Pbx3, Tcf21, Mycn, Sox17, Smad6, Hes1, Etv5, Gli1, Barx1, and Nr0b1. We found the target genes to 18 TFs in the TRRUST database, including Stat3, Stat1, Fos, Irf1, Mycn, Gli1, Junb, Msc, Irf8, Hes1, Zbtb16, Irf7, Nr0b1, Dlx5, Id4, Hopx, Tcf21, and Tead4. According to TRRUST database, 240 DETs were regulated by 15 differentially expressed TFs (Figure 7B). As shown in Figures 7C–F, Irf1 was co-expressed with and may function as a regulator for seven target genes (Cxcl10, Ciita, Il27, Psmb9, Socs1, Tap1, and Tap2); Stat1 was co-expressed with 12 target genes (Cxcl10, Ciita, Il27, Psmb9, Socs1, Tap1, Hsp90aa1, Jak2, Pim1, Tnfsf13b, Irf7, and Irf1); Irf8 was co-expressed with Ncf2; and Stat3 was co-expressed with Il2ra, Mcl1, Pias3, and Usp7.

FIGURE 7

Figure 7. Differentially expressed transcription factors (TFs) and their gene targets. (A) Cluster of all differentially expressed TFs. Upregulated, downregulated, and missing transcripts are denoted by red, green, and gray colors, respectively. (B) Differentially expressed TFs and their target genes. The octagon nodes represent TF and the circular nodes represent target gene of TF with the pie charts inside the nodes denoting the proportion of expressional change in infected samples. Pie chart colors represent the levels of gene expressional changes: red (log2 fold-change ≥ 1), green (log2 fold-change ≤ −1), pink (log2 fold-change between 0 and 1), light green (log2 fold-change between 0 and −1). (C–F) Co-expressed target genes of Irf1, Stat1, Stat3, and Irf8, respectively.

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).

FIGURE 8

Figure 8. Heatmaps of differentially expressed metabolic enzymes, cytokine and cytokine receptors. (A) Heatmap of the differentially expressed cytokines and cytokine receptors. (B) Immune cells regulated by differentially expressed cytokines or cytokine receptors. (C) Heatmap of differentially expressed enzymes involved in metabolism. Red, green, and gray colors represent upregulated genes, downregulated genes, missing genes, respectively. Brown color denotes non-differentially expressed chemokine receptors and immune cells. Abbreviations: Th, T helper cell; Treg, regulatory T cell; iDC, immature dendritic cell; DC, dendritic cell; pDC, plasmacytoid dendritic cell; NK, natural killer; NKT, natural killer T cell; TCM, central memory T cell; TEM, effector-memory T cell; Tfh, T follicular helper cell.

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 metabolism-related and tissue development-related transcripts is prominent in infected tissues (Supplementary Table 4 and Figure 3). This finding may have clinical relevance. During pregnancy, mother-to-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–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.

Funding

Project support was kindly provided by the International Science and Technology Cooperation Project of Gansu Provincial Key Research and Development Program (Grant No. 17JR7WA031), the National Natural Science Foundation of China (Grant No. 31230073), the Elite Program of Chinese Academy of Agricultural Sciences, and the Agricultural Science and Technology Innovation Program (ASTIP) (Grant No. CAAS-ASTIP-2016-LVRI-03).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank BGI-Shenzhen for technical assistance.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.01531/full#supplementary-material

Supplementary Figure 1. Sequencing qualities of the uninfected and infected tissues at 6 dpi. The darker the color the better the global sequencing quality.

Supplementary Figure 2. Sequencing qualities of the uninfected and infected tissues at 18 dpi. The darker the color the better the global sequencing quality.

Supplementary Figure 3. Distribution of sequencing qualities. Vertical axis represents the percentage of clean reads with sequencing quality > Q20. Horizontal axis represents the samples sequenced in the present study.

Supplementary Figure 4. Scale-free topology fit index of coexpression analysis. Scale independence map shows the relationship between soft power and scale free topology model fit of WGCNA analysis. Mean connectivity map shows the relationship between soft power and mean connectivity which summarizes the connection strengths with other genes.

Supplementary Table 1. Primers used for quantitative reverse transcriptase qRT-PCR assay to validate the RNA-seq data.

Supplementary Table 2. Normalized Toxoplasma gondii DNA load in infected pig tissues.

Supplementary Table 3. Differential expression profiles of transcripts across various tissues of the pigs.

Supplementary Table 4. Significantly enriched GO terms and pathways of DETs in the lower expression cluster in most infected tissues.

Supplementary Table 5. Significantly enriched GO terms and pathways of DETs in the higher expression cluster in most infected tissues.

Supplementary Table 6. Significantly enriched GO terms and pathways of DETs in the higher expression cluster in infected liver.

Supplementary Table 7. Significantly enriched GO terms and pathways of DETs in the higher expression cluster at 18 dpi.

Supplementary Table 8. Modules of coexpressed genes.

Supplementary Table 9. Details of the host genes significantly correlated with T. gondii load (HGSCTG).

Supplementary Table 10. Relationship between 330 xenobiotics and differentially expressed metabolism related genes.

References

1. Tenter AM, Heckeroth AR, Weiss LM. Toxoplasma gondii: from animals to humans. Int J Parasitol. (2000) 30:1217–58. doi: 10.1016/S0020-7519(00)00124-7

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Robert-Gangneux F, Darde ML. Epidemiology of and diagnostic strategies for toxoplasmosis. Clin Microbiol Rev. (2012) 25:264–96. doi: 10.1128/CMR.05013-11

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bassols A, Costa C, Eckersall PD, Osada J, Sabria J, Tibau J. The pig as an animal model for human pathologies: a proteomics perspective. Proteomics Clin Appl. (2014) 8:715–31. doi: 10.1002/prca.201300099

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Schluter D, Daubener W, Schares G, Gross U, Pleyer U, Luder C. Animals are key to human toxoplasmosis. Int J Med Microbiol. (2014) 304:917–29. doi: 10.1016/j.ijmm.2014.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hill D, Coss C, Dubey JP, Wroblewski K, Sautter M, Hosten T, et al. Identification of a sporozoite-specific antigen from Toxoplasma gondii. J Parasitol. (2011) 97:328–37. doi: 10.1645/GE-2782.1

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Sher A, Tosh K, Jankovic D. Innate recognition of Toxoplasma gondii in humans involves a mechanism distinct from that utilized by rodents. Cell Mol Immunol. (2017) 14:36–42. doi: 10.1038/cmi.2016.12

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Johnston AC, Piro A, Clough B, Siew M, Virreira Winter S, Coers J, et al. Human GBP1 does not localize to pathogen vacuoles but restricts Toxoplasma gondii. Cell Microbiol. (2016) 18:1056–64. doi: 10.1111/cmi.12579

CrossRef Full Text | Google Scholar

8. Dubey JP, Baker DG, Davis SW, Urban JJ, Shen SK. Persistence of immunity to toxoplasmosis in pigs vaccinated with a non-persistent strain of Toxoplasma gondii. Am J Vet Res. (1994) 55:982–7.

Google Scholar

9. Lin S, Lin Y, Nery JR, Urich MA, Breschi A, Davis CA, et al. Comparison of the transcriptional landscapes between human and mouse tissues. Proc Natl Acad Sci USA. (2014) 111:17224–9. doi: 10.1073/pnas.1413624111

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Lunney JK. Advances in swine biomedical model genomics. Int J Biol Sci. (2007) 3:179–84. doi: 10.7150/ijbs.3.179

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Meurens F, Summerfield A, Nauwynck H, Saif L, Gerdts V. The pig: a model for human infectious diseases. Trends Microbiol. (2012) 20:50–7. doi: 10.1016/j.tim.2011.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Dubey JP, Urban JF Jr. Diagnosis of transplacentally induced toxoplasmosis in pigs. Am J Vet Res. (1990) 51:1295–9.

PubMed Abstract | Google Scholar

13. Dubey JP. Toxoplasmosis in pigs–the last 20 years. Vet Parasitol. (2009) 164:89–103. doi: 10.1016/j.vetpar.2009.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Shiono Y, Mun HS, He N, Nakazaki Y, Fang H, Furuya M, et al. Maternal-fetal transmission of Toxoplasma gondii in interferon-gamma deficient pregnant mice. Parasitol Int. (2007) 56:141–8. doi: 10.1016/j.parint.2007.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jungersen G, Jensen L, Riber U, Heegaard PM, Petersen E, Poulsen JS, et al. Pathogenicity of selected Toxoplasma gondii isolates in young pigs. Int J Parasitol. (1999) 29:1307–19. doi: 10.1016/S0020-7519(99)00078-8

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Senegas A, Villard O, Neuville A, Marcellin L, Pfaff AW, Steinmetz T, et al. Toxoplasma gondii-induced foetal resorption in mice involves interferon-gamma-induced apoptosis and spiral artery dilation at the maternofoetal interface. Int J Parasitol. (2009) 39:481–7. doi: 10.1016/j.ijpara.2008.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wallon M, Peyron F, Cornu C, Vinault S, Abrahamowicz M, Kopp CB, et al. Congenital Toxoplasma infection: monthly prenatal screening decreases transmission rate and improves clinical outcome at age 3 years. Clin Infect Dis. (2013) 56:1223–31. doi: 10.1093/cid/cit032

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Nau J, Eller SK, Wenning J, Spekker-Bosker KH, Schroten H, Schwerk C, et al. Experimental porcine Toxoplasma gondii infection as a representative model for human toxoplasmosis. Mediators Inflamm. (2017) 2017:3260289. doi: 10.1155/2017/3260289

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hou Z, Liu D, Su S, Wang L, Zhao Z, Ma Y, et al. Comparison of splenocyte microRNA expression profiles of pigs during acute and chronic toxoplasmosis. BMC Genomics. (2019) 20:97. doi: 10.1186/s12864-019-5458-y

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li S, Yang J, Wang L, Du F, Zhao J, Fang R. Expression profile of microRNAs in porcine alveolar macrophages after Toxoplasma gondii infection. Parasit Vectors. (2019) 12:65. doi: 10.1186/s13071-019-3297-y

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Zhou CX, Elsheikha HM, Zhou DH, Liu Q, Zhu XQ, Suo X. Dual identification and analysis of differentially expressed transcripts of porcine PK-15 cells and Toxoplasma gondii during in vitro infection. Front Microbiol. (2016) 7:721. doi: 10.3389/fmicb.2016.00721

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhou CX, Zhou DH, Liu GX, Suo X, Zhu XQ. Transcriptomic analysis of porcine PBMCs infected with Toxoplasma gondii RH strain. Acta Trop. (2016) 154:82–8. doi: 10.1016/j.actatropica.2015.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Okomo-Adhiambo M, Beattie C, Rink A. cDNA microarray analysis of host-pathogen interactions in a porcine in vitro model for Toxoplasma gondii infection. Infect Immun. (2006) 74:4254–65. doi: 10.1128/IAI.00386-05

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Gisbert Algaba I, Verhaegen B, Jennes M, Rahman M, Coucke W, Cox E, et al. Pork as a source of transmission of Toxoplasma gondii to humans: a parasite burden study in pig tissues after infection with different strains of Toxoplasma gondii as a function of time and different parasite stages. Int J Parasitol. (2018) 48:555–60. doi: 10.1016/j.ijpara.2017.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Dubey JP. Validation of the specificity of the modified agglutination test for toxoplasmosis in pigs. Vet Parasitol. (1997) 71:307–10. doi: 10.1016/S0304-4017(97)00016-2

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lind P, Haugegaard J, Wingstrand A, Henriksen SA. The time course of the specific antibody response by various ELISAs in pigs experimentally infected with Toxoplasma gondii. Vet Parasitol. (1997) 71:1–15. doi: 10.1016/S0304-4017(97)00010-1

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jiang HH, Huang SY, Zhou DH, Zhang XX, Su C, Deng SZ, et al. Genetic characterization of Toxoplasma gondii from pigs from different localities in China by PCR-RFLP. Parasites Vectors. (2013) 6:227. doi: 10.1186/1756-3305-6-227

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2^−ΔΔCT method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

CrossRef Full Text | Google Scholar

29. Biswas S, Agrawal YN, Mucyn TS, Dangl JL, Jones CD. Biological averaging in RNA-seq. Quant Methods. arXiv:1309.0670 [q-bio.QM]. (2013).

Google Scholar

30. Hill JT, Demarest BL, Bisgrove BW, Gorsi B, Su YC, Yost HJ. MMAPPR: mutation mapping analysis pipeline for pooled RNA-seq. Genome Res. (2013) 23:687–97. doi: 10.1101/gr.146936.112

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tarazona S, Furio-Tari P, Turra D, Pietro AD, Nueda MJ, Ferrer A, et al. Data quality aware analysis of differential expression in RNA-seq with NOISeq R/Bioc package. Nucleic Acids Res. (2015) 43:e140. doi: 10.1093/nar/gkv711

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Williams AG, Thomas S, Wyman SK, Holloway AK. RNA-seq data: challenges in and recommendations for experimental design and analysis. Curr Protoc Hum Genet. (2014) 83: 1–20. doi: 10.1002/0471142905.hg1113s83

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Cong W, Dottorini T, Khan F, Emes RD, Zhang FK, Zhou CX, et al. Acute Toxoplasma Gondii infection in cats induced tissue-specific transcriptional response dominated by immune signatures. Front Immunol. (2018) 9:2403. doi: 10.3389/fimmu.2018.02403

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Williams CR, Baccarella A, Parrish JZ, Kim CC. Trimming of sequence reads alters RNA-Seq gene expression estimates. BMC Bioinformatics. (2016) 17:103. doi: 10.1186/s12859-016-0956-2

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Mahdavi Mashaki K, Garg V, Nasrollahnezhad Ghomi AA, Kudapa H, Chitikineni A, Zaynali Nezhad K, et al. RNA-Seq analysis revealed genes associated with drought stress response in kabuli chickpea (Cicer arietinum L.). PLoS ONE. (2018) 13:e0199774. doi: 10.1371/journal.pone.0199774

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Kim B-M, Ahn SH, Choi NR, Heo J, Kim H, Kwon KH, et al. Transcriptome profiles of Daphnia magna across to the different water chemistry of surface water of the Korean Demilitarized Zone. Toxicol Environ Health Sci. (2017) 9:188–98. doi: 10.1007/s13530-017-0320-6

CrossRef Full Text | Google Scholar

37. Rosani U, Gerdol M. A bioinformatics approach reveals seven nearly-complete RNA-virus genomes in bivalve RNA-seq data. Virus Res. (2017) 239:33–42. doi: 10.1016/j.virusres.2016.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. (2008) 5:621–8. doi: 10.1038/nmeth.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Statist Soc B. (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

40. Lex A, Gehlenborg N, Strobelt H, Vuillemot R, Pfister H. UpSet: visualization of intersecting sets. IEEE Trans Vis Comput Graph. (2014) 20:1983–92. doi: 10.1109/TVCG.2014.2346248

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Chepelev I, Wei G, Wangsa D, Tang Q, Zhao K. Characterization of genome-wide enhancer-promoter interactions reveals co-expression of interacting genes and modes of higher order chromatin organization. Cell Res. (2012) 22:490–503. doi: 10.1038/cr.2012.15

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Kommadath A, Bao H, Arantes AS, Plastow GS, Tuggle CK, Bearson SM, et al. Gene co-expression network analysis identifies porcine genes associated with variation in Salmonella shedding. BMC genomics. (2014) 15:452. doi: 10.1186/1471-2164-15-452

PubMed Abstract | CrossRef Full Text | Google Scholar

43. De Bodt S, Proost S, Vandepoele K, Rouze P, Van de Peer Y. Predicting protein-protein interactions in Arabidopsis thaliana through integration of orthology, gene ontology and co-expression. BMC Genomics. (2009) 10:288. doi: 10.1186/1471-2164-10-288

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Lieberman LA, Banica M, Reiner SL, Hunter CA. STAT1 plays a critical role in the regulation of antimicrobial effector mechanisms, but not in the development of Th1-type responses during toxoplasmosis. J Immunol. (2004) 172:457–63. doi: 10.4049/jimmunol.172.1.457

CrossRef Full Text | Google Scholar

46. Luder CG, Algner M, Lang C, Bleicher N, Gross U. Reduced expression of the inducible nitric oxide synthase after infection with Toxoplasma gondii facilitates parasite replication in activated murine macrophages. Int J Parasitol. (2003) 33:833–44. doi: 10.1016/S0020-7519(03)00092-4

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Tussiwand R, Lee WL, Murphy TL, Mashayekhi M, Kc W, Albring JC, et al. Compensatory dendritic cell development mediated by BATF-IRF interactions. Nature. (2012) 490:502–7. doi: 10.1038/nature11531

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Norose K, Kikumura A, Luster AD, Hunter CA, Harris TH. CXCL10 is required to maintain T-cell populations and to control parasite replication during chronic ocular toxoplasmosis. Invest Ophthalmol Vis Sci. (2011) 52:389–98. doi: 10.1167/iovs.10-5819

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Saeij JP, Frickel EM. Exposing Toxoplasma gondii hiding inside the vacuole: a role for GBPs, autophagy and host cell death. Curr Opin Microbiol. (2017) 40:72–80. doi: 10.1016/j.mib.2017.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Degrandi D, Konermann C, Beuter-Gunia C, Kresse A, Würthner J, Kurig S, et al. Extensive characterization of IFN-induced GTPases mGBP1 to mGBP10 involved in host defense. J Immunol. (2007) 179:7729–40. doi: 10.4049/jimmunol.179.11.7729

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Khan IA, Moretto M, Wei XQ, Williams M, Schwartzman JD, Liew FY. Treatment with soluble interleukin-15Ralpha exacerbates intracellular parasitic infection by blocking the development of memory CD8+ T cell response. J Exp Med. (2002) 195:1463–70. doi: 10.1084/jem.20011915

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Goldszmid RS, Bafica A, Jankovic D, Feng CG, Caspar P, Winkler-Pickett R, et al. TAP-1 indirectly regulates CD4+ T cell priming in Toxoplasma gondii infection by controlling NK cell IFN-gamma production. J Exp Med. (2007) 204:2591–602. doi: 10.1084/jem.20070634

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Hall AO, Beiting DP, Tato C, John B, Oldenhove G, Lombana CG, et al. The cytokines interleukin 27 and interferon-gamma promote distinct Treg cell populations required to limit infection-induced pathology. Immunity. (2012) 37:511–23. doi: 10.1016/j.immuni.2012.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Stumhofer JS, Laurence A, Wilson EH, Huang E, Tato CM, Johnson LM, et al. Interleukin 27 negatively regulates the development of interleukin 17-producing T helper cells during chronic inflammation of the central nervous system. Nat Immunol. (2006) 7:937–45. doi: 10.1038/ni1376

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Clough B, Wright JD, Pereira PM, Hirst EM, Johnston AC, Henriques R, et al. K63-linked ubiquitination targets Toxoplasma gondii for endo-lysosomal destruction in IFNgamma-stimulated human cells. PLoS Pathog. (2016) 12:e1006027. doi: 10.1371/journal.ppat.1006027

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Ahn YH, Park S, Choi JJ, Park BK, Rhee KH, Kang E, et al. Secreted tryptophanyl-tRNA synthetase as a primary defence system against infection. Nat Microbiol. (2016) 2:16191. doi: 10.1038/nmicrobiol.2016.191

CrossRef Full Text | Google Scholar

57. Murofushi Y, Villena J, Morie K, Kanmani P, Tohno M, Shimazu T, et al. The toll-like receptor family protein RP105/MD1 complex is involved in the immunoregulatory effect of exopolysaccharides from Lactobacillus plantarum N14. Mol Immunol. (2015) 64:63–75. doi: 10.1016/j.molimm.2014.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Neve EP, Svensson K, Fuxe J, Pettersson RF. VIPL, a VIP36-like membrane protein with a putative function in the export of glycoproteins from the endoplasmic reticulum. Exp Cell Res. (2003) 288:70–83. doi: 10.1016/S0014-4827(03)00161-7

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Anderson AC, Anderson DE, Bregoli L, Hastings WD, Kassam N, Lei C, et al. Promotion of tissue inflammation by the immune receptor Tim-3 expressed on innate immune cells. Science. (2007) 318:1141–3. doi: 10.1126/science.1148536

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Shindou H, Hishikawa D, Nakanishiu H, Harayama T, Ishii S, Taguchi R, et al. A single enzyme catalyzes both platelet-activating factor production and membrane biogenesis of inflammatory cells. Cloning and characterization of acetyl-CoA: lyso-PAF acetyltransferase. J Biol Chem. (2007) 282:6532–9. doi: 10.1074/jbc.M609641200

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Moessinger C, Kuerschner L, Spandl J, Shevchenko A, Thiele C. Human lysophosphatidylcholine acyltransferases 1 and 2 are located in lipid droplets where they catalyze the formation of phosphatidylcholine. J Biol Chem. (2011) 286:21330–9. doi: 10.1074/jbc.M110.202424

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Teufel M, Saudek V, Ledig JP, Bernhardt A, Boularand S, Carreau A, et al. Sequence identification and characterization of human carnosinase and a closely related non-specific dipeptidase. J Biol Chem. (2003) 278:6521–31. doi: 10.1074/jbc.M209764200

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Kogl T, Muller J, Jessen B, Schmitt-Graeff A, Janka G, Ehl S, et al. Hemophagocytic lymphohistiocytosis in syntaxin-11-deficient mice: T-cell exhaustion limits fatal disease. Blood. (2013) 121:604–13. doi: 10.1182/blood-2012-07-441139

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Malik U, Javed A. FAM26F: an enigmatic protein having a complex role in the immune system. Int Rev Immunol. (2016) 35:1–11. doi: 10.1080/08830185.2016.1206098

CrossRef Full Text | Google Scholar

65. Szklarczyk R, Wanschers BF, Cuypers TD, Esseling JJ, Riemersma M, van den Brand MA, et al. Iterative orthology prediction uncovers new mitochondrial proteins and identifies C12orf62 as the human ortholog of COX14, a protein involved in the assembly of cytochrome c oxidase. Genome Biol. (2012) 13:R12. doi: 10.1186/gb-2012-13-2-r12

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Rehli M, Sulzbacher S, Pape S, Ravasi T, Wells CA, Heinz S, et al. Transcription factor Tfec contributes to the IL-4-inducible expression of a small group of genes in mouse macrophages including the granulocyte colony-stimulating factor receptor. J Immunol. (2005) 174:7111–22. doi: 10.4049/jimmunol.174.11.7111

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Holster T, Pakkanen O, Soininen R, Sormunen R, Nokelainen M, Kivirikko KI, et al. Loss of assembly of the main basement membrane collagen, type IV, but not fibril-forming collagens and embryonic death in collagen prolyl 4-hydroxylase I null mice. J Biol Chem. (2007) 282:2512–9. doi: 10.1074/jbc.M606608200

CrossRef Full Text | Google Scholar

68. Jung DH, Kim KH, Byeon HE, Park HJ, Park B, Rhee DK, et al. Involvement of ATF3 in the negative regulation of iNOS expression and NO production in activated macrophages. Immunol Res. (2015) 62:35–45. doi: 10.1007/s12026-015-8633-5

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Cardone M, Kandilci A, Carella C, Nilsson JA, Brennan JA, Sirma S, et al. The novel ETS factor TEL2 cooperates with Myc in B lymphomagenesis. Mol Cell Biol. (2005) 25:2395–405. doi: 10.1128/MCB.25.6.2395-2405.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Hata S, Doi N, Kitamura F, Sorimachi H. Stomach-specific calpain, nCL-2/calpain 8, is active without calpain regulatory subunit and oligomerizes through C2-like domains. J Biol Chem. (2007) 282:27847–56. doi: 10.1074/jbc.M703168200

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Goldstein EJ, Montoya JG, Remington JS. Management of Toxoplasma gondii infection during pregnancy. Clin Infect Dis. (2008) 47:554–66. doi: 10.1086/590149

CrossRef Full Text | Google Scholar

72. He JJ, Ma J, Elsheikha HM, Song HQ, Huang SY, Zhu XQ. Transcriptomic analysis of mouse liver reveals a potential hepato-enteric pathogenic mechanism in acute Toxoplasma gondii infection. Parasites Vectors. (2016) 9:427. doi: 10.1186/s13071-016-1716-x

PubMed Abstract | CrossRef Full Text | Google Scholar

73. He JJ, Ma J, Elsheikha HM, Song HQ, Zhou DH, Zhu XQ. Proteomic profiling of mouse liver following acute Toxoplasma gondii infection. PLoS ONE. (2016) 11:e0152022. doi: 10.1371/journal.pone.0152022

PubMed Abstract | CrossRef Full Text | Google Scholar

74. He JJ, Ma J, Song HQ, Zhou DH, Wang JL, Huang SY, et al. Transcriptomic analysis of global changes in cytokine expression in mouse spleens following acute Toxoplasma gondii infection. Parasitol Res. (2016) 115:703–12. doi: 10.1007/s00436-015-4792-5

PubMed Abstract | CrossRef Full Text | Google Scholar

75. David CN, Frias ES, Szu JI, Vieira PA, Hubbard JA, Lovelace J, et al. GLT-1-dependent disruption of CNS glutamate homeostasis and neuronal function by the protozoan parasite Toxoplasma gondii. PLoS Pathog. (2016) 12:e1005643. doi: 10.1371/journal.ppat.1005643

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Mendez OA, Koshy AA. Toxoplasma gondii: entry, association, and physiological influence on the central nervous system. PLoS Pathog. (2017) 13:e1006351. doi: 10.1371/journal.ppat.1006351

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Vyas A, Kim SK, Giacomini N, Boothroyd JC, Sapolsky RM. Behavioral changes induced by Toxoplasma infection of rodents are highly specific to aversion of cat odors. Proc Natl Acad Sci USA. (2007) 104:6442–7. doi: 10.1073/pnas.0608310104

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Flegr J. Effects of Toxoplasma on human behavior. Schizophr Bull. (2007) 33:757–60. doi: 10.1093/schbul/sbl074

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Sugden K, Moffitt TE, Pinto L, Poulton R, Williams BS, Caspi A. Is Toxoplasma gondii infection related to brain and behavior impairments in humans? Evidence from a population-representative birth cohort. PLoS ONE. (2016) 11:e0148435. doi: 10.1371/journal.pone.0148435

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Etheridge RD, Alaganan A, Tang K, Lou HJ, Turk BE, Sibley LD. The Toxoplasma pseudokinase ROP5 forms complexes with ROP18 and ROP17 kinases that synergize to control acute virulence in mice. Cell Host Microbe. (2014) 15:537–50. doi: 10.1016/j.chom.2014.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Luder CG, Lang C, Giraldo-Velasquez M, Algner M, Gerdes J, Gross U. Toxoplasma gondii inhibits MHC class II expression in neural antigen-presenting cells by down-regulating the class II transactivator CIITA. J Neuroimmunol. (2003) 134:12–24. doi: 10.1016/S0165-5728(02)00320-X

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Ceravolo IP, Chaves AC, Bonjardim CA, Sibley D, Romanha AJ, Gazzinelli RT. Replication of Toxoplasma gondii, but not Trypanosoma cruzi, is regulated in human fibroblasts activated with gamma interferon: requirement of a functional JAK/STAT pathway. Infect Immun. (1999) 67:2233–40.

Google Scholar

83. Griffith JW, Sokol CL, Luster AD. Chemokines and chemokine receptors: positioning cells for host defense and immunity. Annu Rev Immunol. (2014) 32:659–702. doi: 10.1146/annurev-immunol-032713-120145

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Larson AM, Polson J, Fontana RJ, Davern TJ, Lalani E, Hynan LS, et al. Acetaminophen-induced acute liver failure: results of a United States multicenter, prospective study. Hepatology. (2005) 42:1364–72. doi: 10.1002/hep.20948

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Renton KW. Regulation of drug metabolism and disposition during inflammation and infection. Expert Opin Drug Metab Toxicol. (2005) 1:629–40. doi: 10.1517/17425255.1.4.629

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Ustun S, Aksoy U, Dagci H, Ersoz G. Incidence of toxoplasmosis in patients with cirrhosis. World J Gastroenterol. (2004) 10:452–4. doi: 10.3748/wjg.v10.i3.452

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Shapira Y, Agmon-Levin N, Renaudineau Y, Porat-Katz BS, Barzilai O, Ram M, et al. Serum markers of infections in patients with primary biliary cirrhosis: evidence of infection burden. Exp Mol Pathol. (2012) 93:386–90. doi: 10.1016/j.yexmp.2012.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar