Acute Toxoplasma Gondii Infection in Cats Induced Tissue-Specific Transcriptional Response Dominated by Immune Signatures

RNA-sequencing was used to detect transcriptional changes in six tissues of cats, seven days after T. gondii infection. A total of 737 genes were differentially expressed (DEGs), of which 410 were up-regulated and 327 were down-regulated. The liver exhibited 151 DEGs, lung (149 DEGs), small intestine (130 DEGs), heart (123 DEGs), brain (104 DEGs), and spleen (80 DEGs)-suggesting tissue-specific transcriptional patterns. Gene ontology and KEGG analyses identified DEGs enriched in immune pathways, such as cytokine-cytokine receptor interaction, Jak-STAT signaling pathway, NOD-like receptor signaling pathway, NF-kappa B signaling pathway, MAPK signaling pathway, T cell receptor signaling pathway, and the cytosolic DNA sensing pathway. C-X-C motif chemokine 10 (CXCL10) was involved in most of the immune-related pathways. PI3K/Akt expression was down-regulated in all tissues, except the spleen. The genes for phosphatase, indoleamine 2,3-dioxygenase, Hes Family BHLH Transcription Factor 1, and guanylate-binding protein 5, playing various roles in immune defense, were co-expressed across various feline tissues. Multivariate K-means clustering analysis produced seven gene clusters featuring similar gene expression patterns specific to individual tissues, with lung tissue cluster having the largest number of DEGs. These findings suggest the presence of a broad immune defense mechanism across various tissues in cats against acute T. gondii infection.


INTRODUCTION
The intracellular protozoan Toxoplasma gondii infects almost all warm-blooded animals and approximately one-third of the world's human population, causing toxoplasmosis, a serious illness with fatal consequences in immune-compromised individuals and the unborn fetus (1)(2)(3). The emergence of drug-resistant parasite strains (4) together with the adverse effects of the currently-available drug therapies (5), and the inability to clear chronic infection highlight the need for improved treatment strategies to combat toxoplasmosis in humans and animals.
GRAPHICAL ABSTRACT | The top 10 Differentially Expressed Genes (DEGs) per examined cat tissue. The gene name and symbol, the log2 fold change and P-value for each gene are shown.
T. gondii has a complex lifecycle, wherein the parasite undergoes asexual reproduction in the intermediate host and sexual reproduction in the definitive host (members of the Felidae family). Cats acquire infection by ingesting prey tissue containing parasite cysts or, more rarely, oocysts. Prenatal infection may also occur in cats and humans (2,6). Infected cats discharge oocysts (the product of parasite sexual reproduction in the cat's intestine) in their feces. The ability to accommodate both sexual and asexual reproduction of T. gondii, makes cats a significant source of infection to humans and animals (1,7). T. gondii is a food-borne pathogen acquired via oral infection; however, it has a dispersive nature and can disseminate throughout the cat's body to infect multiple tissues (8)(9)(10)(11)(12)(13)(14).
Recognizing early transcriptional signatures of infection, while knowing the factors that determine tissue susceptibility to T. gondii infection in cats, would assist in planning preventative measures against environmental contamination with oocysts. Previous investigations have provided insights into the transcriptome of many intermediate host species during T. gondii infection (15)(16)(17)(18)(19)(20)(21). However, knowledge of the mechanisms that underpin the feline host transcriptional response to T. gondii remains poorly understood; and no genome-wide expression profiling of multiple tissues from cat, has been reported. Single-tissue gene expression can provide information on how a specific tissue responds to infection; however, understanding the patterns of gene expression across various tissues may advance our understanding of T. gondii molecular-pathogenesis events occurring during acute infection in cats and reveal the mechanisms by which the definitive host counters a complex infection such as T. gondii. RNA sequencing (RNA-Seq) of the whole transcriptome has proved a powerful and versatile tool for global gene expression analysis (22)(23)(24), enabling a comparison of transcriptomes between the merozoite and tachyzoite stages of T. gondii infecting the cat's intestine (25).
Here, we hypothesized that tissues of T. gondii-infected cats exhibit characteristic transcriptional signatures which are dominated by a number of genes, and may be exclusive to a particular tissue with variance across tissues from the same individual. Differential gene expression and gene clustering analyses were performed on six tissues, individually or combined, in cats infected with T. gondii. Our studies provided novel information about the transcriptomic landscape of the major tissue types in cats during acute T. gondii infection and revealed tissue responsive signatures during acute T. gondii infection.

Ethics Approval
All efforts were made to minimize suffering and to reduce the number of animals in the experiment. All work with live T. gondii was performed at biosafety level 2 and the animal experimental protocols were approved by The Animal Administration and Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (Protocol Permit Number: LVRIAEC-2014-001).

Animal Husbandry
Twelve domestic cats (Felis catus; 2 to 3-months-old; Chinese Li Hua breed) were purchased from a local breeder and housed individually in a controlled environment. Six cats each from two litters-were randomly allocated into four groups (2x infected and 2x control) with three cats per group. All of the cats were tested negative for feline immunodeficiency virus and feline leukemia virus using SNAP FIV/FeLV Combo Test (IDEXX, Westbrook, US), and feline calicivirus and feline parvovirus by ELISA KIT (NuoYuan, Shanghai, China) prior to the experiment. Also, cats tested seronegative for the presence of specific anti-T. gondii antibodies using indirect fluorescent antibody test (IFAT). Cats were supplied with a commercial diet (Royal Canine Inc., St. Charles, MO, USA) and water ad libitum during the 3 weeks prior to experimentation. During the experiment cats were fed a maintenance diet, based on their daily energy requirement.

Parasite Strain, Infection and Sample Collection
T. gondii Pru strain (genotype II) was maintained by passage through Kunming mice (26). The number of T. gondii cysts (determined using an optical microscope) was adjusted to 100 cysts per 1ml phosphate buffered saline (PBS, pH 7.4). Each experimental animal was infected with 100 cysts in 1 ml PBS by intragrastric inoculation. The enteroepithelial sexual cycle of T. gondii is completed within 3-10 days post ingestion of T. gondii cysts. This period can be extended to >18 days if cats are infected by oocysts. Control cats were sham-infected with 1 ml PBS only. Six different tissues (brain, heart, liver, lung, small intestine, and spleen) were collected from cats seven days after infection, in order to allow sufficient time for infection to establish. Tissue collection was performed as a terminal procedure from cats under deep isoflurane anesthesia, and unresponsiveness to all stimuli. Collected tissue samples were rinsed in saline, flash frozen in liquid nitrogen and stored at −80 • C until processing.

Confirmation of Infection
Total genomic DNA was extracted from all harvested tissues using TIANamp Genomic DNA kit (TianGen TM , Beijing, China).
A PCR assay targeting T. gondii B1 gene was used to detect T. gondii infection in all cat tissues (27,28). The PCR products were subjected to electrophoresis on ethidium bromide-stained 2% agarose-Tris-acetate-EDTA gels, and the banding pattern was visualized by UV transillumination. All of the electrophoresed PCR products were run with positive and negative control. T. gondii genotype was confirmed by PCR-restriction fragment length polymorphism (RFLP) analysis of the positive amplicons (29).

Read Alignment and Transcript Assembly
Index of the reference genome was built using Hisat2 (v2.1.0) (30). Trimmed paired-end reads were aligned to the reference genome using Hisat2 for expression estimation. StringTie (v1.3.3) was then used to assemble the read alignments into known transcripts for each sample (31). In addition, StringTie also produces a gene abundance table (FPKM and TPM), which was used for clustering analysis.

Tissue-Specific Differential Expression Analysis
A combined read count table at the gene level for all the samples was generated using a python script available from StringTie. The Bioconductor package edgeR (v3.18.1) (32) was used to identify the differentially expressed (DE) genes per-tissue condition (infected vs. uninfected; two biological replicates per condition). Genes with a 5% false discovery rate (FDR < 0.05) and log fold change (logFC) ≥ 1 were considered differentially expressed.

Gene Ontology (GO) and Pathway Analyses
The Bioconductor package GOstats (v2.42.0) (33) was used to test for over-representation of GO terms using a hypergeometric test (hyperGTest). Orthologues for cat gene-sets were found using Ensembl BioMart against human data; then the cat gene orthologues were used for gene ontology analysis. GO terms with a corrected P < 0.05 were considered significantly enriched. Pathway analysis was performed using bioconductor package "pathview" (v1.16.5), which implements Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The significance level of enrichment of KEGG pathways was identified using FDR < 0.05 and a corrected P < 0.05.

K-Means Clustering
The 21,890 genes identified through HiSat2 and StringTie were subjected to clustering analysis. In detail: the FPKM values of the genes from both replicates of each tissue and each treatment (e.g., rep 1 and rep 2 of infected brain tissues) were averaged; the log of the ratio between infected and non-infected conditions for each tissue was calculated; then, two subsequent k-means clustering analyses were performed [MeV_4_8 v10.2 (Multi Experiment viewer)] using Euclidean distance and k = 10 number of clusters. A stringent expression cutoff was applied to each sub-group of genes to discard background noise. Specifically, in each of the six generated clusters where genes showed increased expression levels concentrated on a single infected tissue, only those genes with FPKM values >1 in the specific infected condition were considered as expressed and selected for further analysis. To visualize gene expression in each cluster, bean plots (which represent the actual distribution of the individual data sets) were produced using BoxPlotR: a web-tool for generation of box plots.

Quantitative Real Time (RT)-PCR
Total RNA was extracted from various cat tissues using the TRIzol method (Invitrogen) and reverse-transcripted to single strand cDNA using the GoScript TM Reverse Transcription System (Promega, MI, USA). GoTaq R qPCR Master Mix (Promega, MI, USA) was used to perform RT-PCR reactions in a QIAGEN's real-time PCR cycler (Rotor-Gene Q). Amplification reactions were performed under the following conditions: 2 min at 95 • C, 40 cycles of 95 • C for 15s, 55 • C for 30s, and 72 • C for 30s. All quantitative measurements were carried out in triplicate and normalized to the housekeeping gene glyceraldehyde-3phosphatedehydrogenase (GAPDH) for every reaction (34). Twelve significant DEGs were selected to validate the sequencing data. Primers used for RT-PCR are listed in Table 1. The mRNA fold change was calculated using the following equations (35):

Presence of T. Gondii Infection in Cats
Parasite dissemination from the intestine to other tissues, both close to the inoculation site (small intestine, liver, and spleen) and to distantly placed organs (brain, heart, and lungs) was confirmed seven days post infection. T. gondii B1 gene-based PCR analysis detected the parasite DNA in the brain, heart, liver, lung, small intestine, and spleen of all infected cats ( Figure S1). The parasite load did not seem to vary across the cat tissues. However, it is possible that the parasite load may vary over the course of infection. T. gondii DNA was not detected in any tissue of the control cats. All positive amplicons that were characterized by RFLP produced a restriction fragment pattern that correlated with T. gondii genotype II.

General Features of the Transcriptome Data
Transcriptome-sequencing analysis generated ∼143 million sequence reads (125 bp in length) from 24 libraries. After quality control analysis and the removal of low quality reads, ∼170 Gb clean reads were obtained, with an average of 7 Gb clean reads/tissue. Less than 95% of the clean reads had Phred-like quality scores at the Q20 level and GC content of about 50% ( Table S1). The clean reads were mapped to the genome of Felis catus. The majority of the clean reads were distributed in the exon region, with fewer in the intergenic region and the intron region. Approximately 80% of the clean reads were unique (Table S2); and subsequent analyses were based on these uniquely mapped reads.

Infection Induced Significant Alterations in Gene Expression
We investigated the distribution of gene expression values across the six tissues by fragment Per Kilobase of exon per Million mapped reads (FPKM) ( (151 genes) compared to lung (149 genes), small intestine (130 genes), heart (123 genes), brain (104 genes), and spleen (80 genes) (Figure 1). DEGs of each tissue are listed in Table S4. We also used quantitative real-time PCR to validate the expression levels of representative genes, across cat tissues, detected by RNA-sequencing analysis (Figure 2).

Gene Ontology (GO) Enrichment and Functional Annotation Analyses
GO enrichment analysis was used to identify the significantly enriched GO terms in all DEGs using the GOseq R package. The enriched GO terms of each tissue are shown in Table S5.
All DEGs were mapped to terms in the KEGG database. DEGs and KEGG pathways related to immune response were highly represented and are summarized in Table S6.

Signatures of Gene Co-expression
Differences were detected between the transcriptomes of cat tissues, with significant variations in gene expression between infected and uninfected tissues (Table S4). Gene co-expression analysis of DEGs indicated that phosphatase and indoleamine 2,3-dioxygenase (IDO) were co-expressed in five tissues (brain, heart, liver, small intestine, and spleen); while the HES1 (Hes Family BHLH Transcription Factor 1) was co-expressed in brain, heart, liver, and small intestine. Guanylate-binding protein 5 was detected in heart, liver, lung, and spleen (Figure 3, Table S7). Expression patterns and pathways associated with co-expressed DEGs across the various cat tissues are shown ( Figure S2).

Global Gene Regulation During T. Gondii Infection
We employed cluster analysis, which was based on the partition of genes into clusters according to the log-ratios of gene expression between infected and uninfected tissues using a two-tier K-means algorithm. This analysis identified a set of 10 gene groups each characterized by a unique pattern. The patterns corresponding to individual clusters are visualized as line plots (Figure 4). Seven of these 10 clusters showed significant expression levels: spleen (cluster 2); brain (cluster 3); lung (cluster 7); liver (cluster 8), and heart (cluster 9), included up-regulated genes. Small intestine tissue produced two clusters with opposite regulation: cluster 4 (up-regulated genes) and cluster 10 (downregulated genes). Two clusters (5 and 6) did not exhibit gene expression patterns specific to any tissue. To isolate the most representative genes in each cluster, a second clustering analysis was performed, focusing on genes contained in each of the seven clusters identified by the initial cluster analysis. Ultimately, seven groups of tissue-specific genes were identified, and are shown in heat maps and graphical formats, based on patterns of expression of individual genes across various tissues (Figures 5A,B).

Refined Gene Clusters After FPKM Filtering
Here, a stringent expression cutoff was applied to each subgroup of genes to remove background noise. In each of the six clusters where genes showed increased expression concentrated on a single infected tissue, only those genes with FPKM values >1 in the infected condition were selected for further analysis. Five   of these clusters contained up-regulated genes in the spleen (6 genes), brain (13 genes), lung (75 genes), liver (1 gene), and heart (6 genes). However, the small intestine produced one cluster of 19 up-regulated genes and one cluster of 21 down-regulated genes. GO analysis of genes within these seven clusters provided new insights into biological processes, cellular processes and molecular functions regulated in cats during T. gondii infection. The five most highly associated biological process terms with the enrichment P-values of each cluster are shown ( Figure 5C).
KEGG analysis of the DEGs within these seven final clusters was performed to gain more insight into the biological processes influenced during infection. We found no similarity in gene annotation in any cluster in any of the tissues examined. In small intestine cluster, KEGG analysis identified networks associated with hematopoietic cell lineage, metabolism of xenobiotics by cytochrome-P, chemical carcinogenesis and cytokine-cytokine receptor interaction. Lung cluster showed high expression levels of genes involved in nicotine addiction, biosynthesis of unsaturated fatty acids, long-term potentiation, glycosaminoglycan biosynthesis -heparan sulfate/heparin, cell adhesion molecules, insulin secretion, GABAergic synapse, fatty acid metabolism, amyotrophic lateral sclerosis, and adrenergic signaling in cardiomyocytes. The heart cluster showed high expression levels in complement and coagulation cascades and platelet activation. Small intestine cluster showed high expression levels of genes involved in metabolic metabolism, drug metabolism and fatty acid and amino acid metabolism.
Genes with different expression patterns across the seven clusters are shown in Table S8. All of the analyzed clusters contained genes annotated to GO terms from the three ontologies at roughly equivalent levels. However, we found considerable bias in some tissues in regards to the representation of clustered genes among the three ontologies ( Figure S3). The high proportion of genes annotated to biological process, reflects a trend toward clustering of metabolic and immune pathways. The most striking example was in the lung cluster, where 253 terms were found in the biological process.

Tissue-Specific Changes in Gene Expression
Even though the number of genes in each cluster was small, large-magnitude gene expression changes between infected and uninfected was observed in most of the tissues (Figure 6). The liver gene cluster had the highest magnitude (7.68 fold), followed by lung (5.44 fold), brain (3.58 fold), small intestine cluster 4 (3.55 fold), spleen (3.09 fold), and heart (2.98 fold); whereas small intestine tissue cluster 10 had a magnitude of −3.22 fold. Apart from the liver with only on gene, lung, brain, small intestine, and heart showed large magnitude transcriptional changes occurring in a small number of genes. This suggests that changes induced by T. gondii infection seem to rely on a small number of genes, but with large transcriptional changes. The lung was the tissue showing the greatest number of genes with the largest magnitude of transcriptional change. Among all highly expressed genes, the FIGURE 4 | k-Means clustering patterns of the differentially expressed genes across cat tissues. All DEGs were clustered into 10 groups using k-means clustering method and visualized with TM4 software. The pink line shows average expression z-scores to visualize the dominant expression trend of each cluster. Each line in the figure represents an expression value of the corresponding gene. The horizontal axis indicates the type of cat tissues and the vertical axis is the log2 expression ratio. Fold expression changes between various tissues from infected and uninfected cats were calculated using the log2 ratios. The numbers of genes for each cluster are indicated. The clusters included DEGs whose expression was either up regulated (clusters 2, 3, 4, 7, 8, and 9) or down regulated (cluster 10). There were no changes in gene expression observed in clusters 5 and 6. DEGs in cluster 1 was not considered due to lack of specificity to a single tissue. metallothionein 3 (MT3) gene showed the highest expression in the infected lung, and it was also one of the genes with the highest log ratio between infected and uninfected tissues. Metallothionein-3 protein encoded by MT3 gene binds metals both in natural (such as Zn, Cu, Se) and xenobiotic (such as Ag, Cd) conditions, conferring a protective role against metal toxicity and oxidative stress.

DISCUSSION
The influence of T. gondii infection on the transcriptome of different tissue of cats is largely unknown. In this study, we examined the transcriptomes of liver, lung, small intestine, heart, brain, and spleen tissues of cats, seven days post infection. Intertissue transcriptome comparison revealed that the levels of DEGs were highest in the liver (151 genes) and lungs (149 genes), indicating a significant gene transcriptional response to infection in the liver and lung compared to other feline tissues. This result supports previous reports, wherein liver and lung seemed to be the most likely tissues to be involved in clinical cases of toxoplasmosis with a rapid fatal outcome (11-13, 34, 36).
It is plausible that the timing and duration of infection can influence host tissue transcriptional response to the parasite, as previously indicated by the difference in the transcriptomes of mouse brain between acute and chronic T. gondii infection (17), and this could alter cytokine responses of the host to infection.
The transcriptional landscape of infected cat tissues was dominated by an immune gene expression signature, wherein, cytokine-cytokine receptor interaction, Jak-STAT signaling pathway, NOD-like receptor signaling pathway, NF-kappa B signaling pathway, MAPK signaling pathway, T cell receptor signaling pathway and the cytosolic DNA sensing pathway, were amongst the up-regulated immune pathways in almost all tissues. The importance of Toll-like receptor signaling in controlling T. gondii infection has been established (37). Many of the upregulated genes (e.g., CXCL10, SOCS3, MAPK13, CXCL9, CD2, CSF2RA, PI4K2B, IGF1, PFKFB1, MMP7, FZD8, TNFSF10, and RelB) identified in the livers are involved in immune-related pathways. Some of these genes are involved in the development of Natural Killer (NK) and adaptive T cell responses, leading to the production of Interferon gamma (IFN-γ) and resistance to infection (38). T. gondii is very FIGURE 5 | Hierarchical and K-mean clustering analysis of the DEGs within clusters. Differentially expressed, co-regulated genes in each cluster were grouped using k-means clustering. Average cluster size varied considerably among tissues with lung containing the largest cluster with 75 genes. The smallest cluster was found in (Continued) Frontiers in Immunology | www.frontiersin.org efficient in manipulating the host's immune defense (39,40); and previous studies have indicated that host immune response is a key determinant of clinical outcome following infection with T. gondii (16-20, 41, 42). Although spleen had the lowest number of DEGs (80 genes) compared to other tissues, KEGG analysis identified multiple immune signaling pathways that were influenced by infection. Most of the immune pathways were upregulated in spleen of infected cats; such as Toll-like receptor signaling pathway, cytosolic DNA-sensing pathway, RIG-I-like receptor signaling, TNF signaling pathway, FoxO signaling pathway, chemokine signaling pathway, and PI3K/Akt signaling pathway. Of note, C-X-C motif chemokine 10 (CXCL10) was involved in most of the immune-related pathways. Our previous study demonstrated that CXCL10 was up-regulated in pig peripheral blood mononuclear cells during early T. gondii infection (20,21). CXCL10 is a chemokine secreted from cells stimulated with IFN-γ, and plays an important role in chemo-attraction of immune cells (43). This result indicates that T. gondii influences chemokine gene expression in spleen during early T. gondii infection.
T. gondii is a food-borne pathogen and infection is initially established in the small intestine by consuming prey containing parasite cysts or oocyst-contaminated water (44); leading to enteropathy (45). Therefore, successful infection of the host intestine is essential for subsequent parasite dissemination to different tissues. Others have showed that in vitro infection of rat intestinal epithelial cells can trigger an inflammatory response characterized by Tumor Necrosis Factor alpha (TNFα) signaling via NF-κB (46). However, feline host factors that are influenced by intestinal infection remain largely unknown. In our study, the global transcriptomic changes in the intestine of the definitive feline host, have been investigated. A total of 130 DEGs in the small intestine of infected cats were identified; of which 75 were down-regulated and 55 up-regulated. Some of the DEGs were involved in immune processes and signaling pathways, including LOC100302541, TNFSF18, CCL20, TNFRSF6B, IFNG, SOCS3, ICAM1, CD36, and FGF19. Other immune pathways upregulated in liver and spleen tissues were also up-regulated in the small intestine.
An earlier study investigating the transcriptome of mouse brain during T. gondii infection reported an increased expression of genes involved in immune responses and cell activation (16). Also, host immune and inflammatory response was the major feature of genes affected by T. gondii infection of mouse peritoneal cells at five days post infection (47). We also showed that the expression of genes and signaling pathways involved in host immune response and cell fate, such as PI3K-Akt signaling pathway, Hippo signaling pathway and MAPK signaling pathway, was altered in the cat brain. Transcriptional signatures observed in the cat brain tissue showed also that Notch signaling pathway is involved in T. gondii neuro-pathogenesis. Previous T. gondii host/pathogen KEGG pathway interactome analysis suggested the involvement of six genes of the Notch pathway in psychiatric/neurological disorders (45,48). Notch signaling interacts with other signaling pathways, including phosphatidylinositol-3-kinase (PI3K)/serine/threonine kinase (Akt) and NF-κB to regulate cell fate. PI3K/Akt signaling pathway regulates diverse cellular activities related to cell growth, metabolism, migration, and apoptosis (49). Notch signaling plays an important role in various facets of T. gondii pathogenesis (50). The possibility that PI3K/Akt signaling pathway participates in promoting T. gondii survival and proliferation (51) and in mediating cell survival and blockage of apoptotic responses during T. gondii infection (50), suggests a link between T. gondii and this signal pathway. This finding supports an earlier observation wherein many T. gondii strains were found able to regulate genes enriched for processes involved in cell cycle regulation in murine macrophages (15).
T. gondii exploits heterotrimeric Gi-protein-mediated signaling to activate PI3K, leading to phosphorylation of downstream serine/threonine kinase AKT and extracellular signal-regulated protein kinases 1/2 (ERK1/2), and the inhibition of apoptosis (52). In our study, down-regulation of PI3K/AKT signaling in all tissues, except spleen, was detected and this can enhance apoptosis of host tissue and limits parasite growth. Phosphatase (a negative regulator of PI3K/Akt signaling, by converting PIP3 to PIP2, an opposite action to PI3K) can interfere with a number of cellular functions, such as cell proliferation and cell-cycle progression (53), was co-expressed in five cat tissues. These results indicate diverse roles played by PI3K/Akt signaling in T. gondii-host interaction. Further investigation into how the PI3K/Akt pathway interacts with other signaling mediators is required. NF-κB and mitogen-activated protein kinase (MAPK) signaling pathways are involved in the innate immune response to T. gondii (54). Altered MAPK signaling has been also implicated in toxoplasmosis of mice (55,56) and humans (57). Considering that NF-κB and MAPK are downstream effectors of Akt, it is of significance to clarify whether this pathway influences the fate of infected cells via the regulation of NF-κB and MAPK.

Gene co-association Across Tissues
Co-expression analysis of DEGs indicated that phosphatase gene expression overlapped in five tissues; including brain, heart, liver, small intestine, and spleen. Transcripts for the inner-membrane complex (IMC) protein phosphatase have been involved in gene expression and cell division (46). Indoleamine 2,3-dioxygenase (IDO) was also co-expressed in five tissues: heart, liver, lung, small intestine, and spleen, but not in the brain tissue. IDO, the rate limiting catabolic enzyme in the degradation pathway of tryptophan, initiates the production of tryptophan degradation products, which exert important immuno-regulatory functions. IDO, through T-cell functions and other mechanisms (58), modulates pathophysiological processes, such as antimicrobial and antioxidant activities, and immune-regulation. IFNγ produced in response to T. gondii infection induces IDO enzyme to degrade L-tryptophan, an amino acid for which T. gondii is auxotrophic (59,60). IFNγinduced L-tryptophan starvation was believed to trigger T. gondii clearance via noncanonical, ubiquitin-mediated autophagy (61). Guanylate-binding protein (GBP) was detected in the heart, liver, lung, and spleen. GBP, IFN-γ-inducible effector, is a member of the GTPase family, and plays an essential role in mediating host defense against T. gondii (62). Human guanylate-binding protein 1 (hGBP1) functions against T. gondii infection in human MSCs (hMSCs). T. gondii replication can be significantly inhibited by the recruitment of hGBP1 to the parasitophorous vacuole (PV) membrane in IFN-γ-stimulated hMSCs (63). In mice, the recruitment of mGBP2 to T. gondii-containing PV was essential for controlling T. gondii replication (64).

Patterns of Gene Clusters
Genes clustered according to their pattern of expressions across the six tissues examined, led to the identification of seven clusters featuring expressions differentiated across tissues and infected/non-infected conditions. We also detected tissuespecific variations in the percentage of clustered genes, and in the properties of gene clusters including functional annotation and magnitude of gene expression. The differences in each tissue's response to infection may imply that some tissuespecific defense mechanisms exist in order to maintain the balance between enhancing the host's response to infection and promoting the parasite's survival. The genes within these clusters after a further processing step imposing a stringent expression cutoff to avoid any background expression noise were analyzed for differential transcriptional changes between infected and uninfected samples (magnitude). Interestingly, the lung besides showing the greater number of genes showed also the largest magnitude of transcriptional changes between infected and noninfected conditions, suggesting that changes induced by infection in lungs seem to rely on genes, which had large transcriptional changes.

CONCLUSION
We used two complementary approaches to characterize alterations in tissue-specific gene expression in cats infected with T. gondii. Our results revealed considerable transcriptional differences between cat body tissues, and between infected and healthy cats. We identified significant tissue-specific differences in gene expression, and in gene cluster content and functional annotations. The differences in gene expression and gene clusters may result from tissue-specific differences in the defense processes that shape the host-pathogen interaction. Our data also underlined the importance of immune and inflammatory response to T. gondii infection, regardless of tissue types. Genes and pathways discovered in this study, should serve as a basis for further understanding the cellular and molecular basis of cat response to T. gondii. These results may assist the selection of biomarkers useful for developing new diagnostic tools or therapeutic interventions to control toxoplasmosis in cats.

ACCESSION NUMBER(S)
The RNA-seq data reported in this study have been deposited in GenBank's Short Read Archive (SRA) database under BioProject number PRJNA296557.

AUTHOR CONTRIBUTIONS
X-QZ, HE, and WC designed the experiment. WC supervised the experimental infection. WC, F-KZ, C-XZ, J-JH, and X-XZ performed the experiments. WC, HE, TD, FK, and RE contributed reagents, materials, and analysis tools. TD, FK, and RE developed computational algorithms and performed the bioinformatics analysis. WC, HE, and X-QZ wrote the paper. All authors commented on the manuscript.

ACKNOWLEDGMENTS
We thank Hind Mamdowh Saad for excellent drawing of the Graphical Abstract.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu. 2018.02403/full#supplementary-material Figure S1 | Agarose gel electrophoresis of PCR amplicons after nested amplification of Toxoplasma gondii B1 gene-specific fragment from cat tissue DNA. ∼96-bp products of B1 gene were amplified from ∼193-bp B1 PCR products, originally generated from genomic DNA extracted from different cat tissues, using the nested primers 5 ′ -TGCATAGGTTGCAGTCACTG-3 ′ and 5 ′ -GGCGACCAATCTGCGAATACACC-3 ′ . Samples were analyzed by electrophoresis through 2% (wt/vol) agarose gels. Gels were stained with ethidium bromide and DNA was visualized under UV. Lanes: M, DNA ladder marker (TAKARA, China); 1, positive control; 2, negative control without DNA template; 3-8, positive PCR products from brain, heart, liver, lung, spleen and small intestine of infected cats; 9-14, negative results of samples obtained from the equivalent tissues of uninfected cats. The numbers to the left refer to the size (bp) of marker DNA fragments.