Original Research ARTICLE
Dual Identification and Analysis of Differentially Expressed Transcripts of Porcine PK-15 Cells and Toxoplasma gondii during in vitro Infection
- 1National Animal Protozoa Laboratory and College of Veterinary Medicine, China Agricultural University, Beijing, China
- 2State Key Laboratory of Veterinary Etiological Biology, Key Laboratory of Veterinary Parasitology of Gansu Province, Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Lanzhou, China
- 3Faculty of Medicine and Health Sciences, School of Veterinary Medicine and Science, University of Nottingham, Loughborough, UK
- 4Jiangsu Co-innovation Center for the Prevention and Control of Important Animal Infectious Diseases and Zoonoses, Yangzhou University College of Veterinary Medicine, Yangzhou, China
Toxoplasma gondii is responsible for causing toxoplasmosis, one of the most prevalent zoonotic parasitoses worldwide. The mechanisms that mediate T. gondii infection of pigs (the most common source of human infection) and renal tissues are still unknown. To identify the critical alterations that take place in the transcriptome of both porcine kidney (PK-15) cells and T. gondii following infection, infected cell samples were collected at 1, 3, 6, 9, 12, 18, and 24 h post infection and RNA-Seq data were acquired using Illumina Deep Sequencing. Differential Expression of Genes (DEGs) analysis was performed to study the concomitant gene-specific temporal patterns of induction of mRNA expression of PK-15 cells and T. gondii. High sequence coverage enabled us to thoroughly characterize T. gondii transcriptome and identify the activated molecular pathways in host cells. More than 6G clean bases/sample, including >40 million clean reads were obtained. These were aligned to the reference genome of T. gondii and wild boar (Sus scrofa). DEGs involved in metabolic activities of T. gondii showed time-dependent down-regulation. However, DEGs involved in immune or disease related pathways of PK-15 cells peaked at 6 h PI, and were highly enriched as evidenced by KEGG analysis. Protein-protein interaction analysis revealed that TGME49_120110 (PCNA), TGME49_049180 (DHFR-TS), TGME49_055320, and TGME49_002300 (ITPase) are the four hub genes with most interactions with T. gondii at the onset of infection. These results reveal altered profiles of gene expressed by PK-15 cells and T. gondii during infection and provide the groundwork for future virulence studies to uncover the mechanisms of T. gondii interaction with porcine renal tissue by functional analysis of these DEGs.
Infection with the apicomplexan protozoan parasite Toxoplasma gondii can cause severe disease and even death in immune-compromised individuals and in congenital infections (Elsheikha, 2008). However, association between parasite genotype, host genetic background and severity of the disease may occur in healthy subjects (Bela et al., 2012; McLeod et al., 2012; Xiao and Yolken, 2015). T. gondii strains from Europe and North America have been reported to belong to three (type I, II, and III) main evolutionary lineages (Howe and Sibley, 1995; Khan et al., 2011; McLeod et al., 2012). Genotypes not belonging to the three lineages have been detected in South America (Pena et al., 2008). Recent studies revealed even more genetic diversities, which seem to be driven via genetic recombination events that occur during the sexual phase of the life cycle in the gut epithelium of the definitive felid host (Minot et al., 2012).
What makes T. gondii so special compared to other apicomplexan protozoa is its ability to infect any nucleated cell types in virtually all warm-blooded animals (Dubey et al., 1998; Dubremetz, 1998; Schlüter et al., 2014; Yarovinsky, 2014). Successful infection of T. gondii tachyzoites depends on their ability to contend with immune responses mounted by the hosts they infect. T. gondii organisms have a remarkable ability to manipulate host cell biological process to their own advantage and to evade both innate and adaptive host immune defenses (Hunter and Sibley, 2012; Coombes and Hunter, 2015). Hence, if T. gondii tachyzoites are to replicate intracellularly and survive long enough to effectively establish infection they need to exploit host processes that are beneficial to their metabolic, anti-apoptotic and pathogenetic functions. The host cells, on the other hand, employ several strategies, such as repair and stress pathway, to adapt to and mitigate the damage caused by the parasite (Gazzinelli et al., 2014).
Interaction between T. gondii and host cells is complex due to the numerous parasite and host factors mediating this interaction. T. gondii invades host cells and establishes a parasitiphorous vacuole (PV) in the host cell cytoplasm in which it resides (Pissuwan et al., 2007). Within the PV, T. gondii undergoes endodyogeny, assembling two daughter parasitic cells within each parental parasite in every mitotic cell cycle (~7–10 h post infection; Hu et al., 2002). The main events of T. gondii cell cycle, G1, single S and mitotic (S/M) phases, and cytokinesis (C) phase (Butler et al., 2014), are associated with a coordinated program of gene expression as revealed by microarray and RNA-Sequencing analyses (Behnke et al., 2010; Gaji et al., 2011). During this intracellular replication cycle, specific host cellular state can influence T. gondii gene expression (Radke et al., 2006) and T. gondii secreted effectors could modulate host gene expression leading to an altered host cell microenvironment to which it subsequently responds. The correlations between host and T. gondii gene expression clusters (Melo et al., 2013) clearly support the presence of coordinated cross-talks that mediate interaction between T. gondii and host cells.
The remarkable capability of T. gondii to exploit surrogate host cells have spurred extensive research on the relationship between this parasite and its host using various technologies, including proteomics (Nelson et al., 2008; Zhou et al., 2011, 2013), metabolomics (Tymoshenko et al., 2015) and transcriptomics (Melo et al., 2013; He et al., 2016). Specifically, high-throughput RNA-Sequencing can provide a powerful method for identifying novel RNAs and for studying the expression profiles of genes in different biological samples (Melo et al., 2013), such that hundreds of differentially expressed genes (DEGs) can be identified and quantified in a single analysis. The high sensitivity of RNA-seq can facilitate the precise identification of the DEGs within particular cells or tissues (Sansom et al., 2014).
The pathogenesis of T. gondii in various body tissues has been extensively studied, but little is known about how T. gondii interacts with host renal tissue. Also, although pig is a natural host of T. gondii and represents a major source of zoonotic infection (Tenter et al., 2000; Steinparzer et al., 2015), knowledge of the pathophysiology of toxoplasmosis in this intermediate host species has been very limited. Further, little is known about DEGs and the role they play during T. gondii infection. To address these needs, the present study was performed to investigate the concurrent transcriptomic changes of T. gondii and porcine kidney (PK-15) host cells spanning the first 24 h of infection using next-generation mRNA sequencing (RNA-seq). Quantitative real-time PCR was also performed to validate the transcriptome sequencing data. The obtained transcriptomic sequences were mapped to reference genomes, and were analyzed to identify unique regulatory pathways that are associated with T. gondii during the first 24 h of infection. Understanding of these pathways may provide information likely to be critical for the development of rationally designed T. gondii therapeutics and vaccines.
Materials and Methods
Porcine kidney (PK-15) cells were purchased from National Platform of Experimental Cell Resources for Sci-Tech (Beijing, China) and cultured in Dulbecco's modified Eagle medium (DMEM, HyClone, China) supplemented with 10% heat-inactivated fetal bovine serum (FBS, Gibco, USA), 100 units/mL penicillin and 100 μg/mL streptomycin (Invitrogen), and incubated at 37°C in a humidified 5% CO2 atmosphere.
In vitro Parasite Cultivation and Infection of PK-15 Cells
Type I virulent T. gondii cultures of RH strain tachyzoites were maintained in our lab in immortalized PK-15 epithelial cells. Briefly, monolayers of PK-15 cells, at 80% confluence in 25-cm2 culture flasks, were infected with T. gondii tachyzoites at a multiplicity of infection (MOI) of 5. Tachyzoites were harvested from host feeder cells when 70–80% of infected PK-15 epithelial cells had lyzed. The egressed tachyzoites and the remaining infected cells were harvested with a cell scraper, passed through a syringe and 22-gauge needle to rupture remaining intact host cells, washed once in cold sterile phosphate buffered saline (PBS), and centrifuged twice at 5000 × g for 10 min to eliminate cell debris. The obtained parasite pellet was washed twice with cold PBS with centrifugation at 350 × g for 10 min (Qu et al., 2013). The final purified tachyzoites were re-suspended in PBS. The enriched parasite preparation was counted using a hemocytometer and the concentration of tachyzoites was adjusted by adding fresh medium to achieve a desired MOI of 5, i.e., five tachyzoites per cell. Infected PK-15 cells were incubated for 1 h, followed by removal of the culture medium, washing three times with PBS to remove unattached tachyzoites and addition of fresh culture medium. Infected cultures were incubated and at 1, 3, 6, 9, 12, 18, or 24 h PI samples from infected culture and controls were collected and processed for sequencing analysis. Mock-infected cells and purified parasites were used as controls. Infection experiment with T. gondii was performed under biosafety level 2 (BSL-2) conditions. All infections were performed at 37°C in an atmosphere of 5% CO2.
RNA Isolation, cDNA Library Preparation and Illumina Deep Sequencing
RNA was extracted from samples at various times post infection by using the PrimeScript™ 1st Strand cDNA Synthesis Kit (Takara Bio, Shiga, Japan) according to the manufacturer's instructions. Total RNA was treated with PQI DNase (Promega, MI, USA) to eliminate the residual DNA contamination. RNA degradation and contamination were monitored on 1% agarose gels. The purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and the integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Finally RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA).
The cDNA synthesis was performed with 3 μg of total RNA using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's recommendations. Briefly, mRNA was selected by using poly-T oligo-attached magnetic beads and then fragmented using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). Reverse transcription was performed using random hexamer primer and M-MuLV Reverse Transcriptase with the mRNA fragments as templates. Double stranded cDNA fragments were obtained using DNA Polymerase I and RNase H. Remaining overhangs were blunt-ended and the cDNA fragments were ploy-adenylated. NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. The library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA) to obtain cDNA fragments of preferentially 150~200 bp in length. Purified DNA fragments were then enriched using Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer in a 15-cycle PCR reaction. PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
Before sequencing, the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. Sequencing runs were performed on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated.
Raw sequencing reads (raw data) were processed through in-house perl scripts and the read counts were extracted. The raw reads were then filtered to remove low-quality sequences (Q < 20), empty reads, and reads containing adapter or poly-N to gain clean reads. All subsequent analyses were based on the clean reads with high quality. Reference genomes and gene model annotation files were downloaded (ftp://ftp.ensembl.org/pub/release-76/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa10.2.dna.toplevel.fa.gz;ftp://ftp.ensemblgenomes.org/pub/release-23/protists//fasta/toxoplasma_gondii/dna/Toxoplasma_gondii.ToxoDB-7.1.23.dna.toplevel.fa.gz). Index of the reference genome was built using Bowtie v2.2.3 and paired-end clean reads were aligned to the reference genome sequence using TopHat v2.0.12 (Zhang et al., 2014). Picard-tools v1.96 and samtools v0.1.18 were used to sort mark duplicated reads and reorder the bam alignment results of each sample (Liu et al., 2015a). GATK2 (v3.2) software was used to perform SNP calling (DePristo et al., 2011). The Cufflinks v2.1.1 Reference Annotation Based Transcript (RABT) assembly method was used to construct and identify both known and novel transcripts from TopHat alignment results. Alternative splicing (AS) events were classified to 12 basic types by the software Asprofile v1.0 (Florea et al., 2013). The number of AS events in each sample was estimated separately. To annotate all unigenes obtained, we searched them against the protein databases NR, Swiss-Prot, KEGG, COG and Toxo DB using blastx (BLAST, the basic local alignment search tool; Liu et al., 2015b).
Differential Gene Expression Analysis
To detect differentially expressed genes (DEGs), the number of clean reads assigned to a gene was counted using HTSeq v0.6.1. The read counts were then normalized into the values of fragments per kilobase of exon per million fragments mapped (FPKM) (Trapnell et al., 2010). DEGSeq R package (1.20.0) was utilized to perform differential expression analysis of two conditions (Robinson et al., 2010). Corrected p-value (q-value) < 0.005 and |log2 (fold change)| > 1 were set as the threshold for significantly DEGs. Subsequently, the enrichment analysis of Gene Ontology (GO) and KEGG pathways was performed based on these DEGs by using GOseq R package and KOBAS 2.0 software (http://kobas.cbi.pku.edu.cn/home.do) (Xie et al., 2011). Finally, predicted Protein-Protein Interactions (PPI) analysis was performed based on the STRING database (Szklarczyk et al., 2015).
Verification by Quantitative Real Time RT-PCR
Total RNA was extracted using TRIzol method (Invitrogen) and then was reverse-transcripted to single strand cDNA using GoScript™ Reverse Transcription System (Promega, MI, USA) according to the manufacturer's instructions. GoTaq® qPCR Master Mix (Promega, MI, USA) was used according to the manufacturer's protocols to perform RT-PCR reactions on QIAGEN's real-time PCR cycler, the Rotor-Gene Q. The amplification reactions were performed with the following conditions: 2 min at 95°C, 40 cycles of 95°C for 15 s, 55°C for 30 s, 72°C for 30 s. All quantitative measurements were carried out in triplicate and normalized to the internal glyceraldehyde-3-phosphatedehydrogenase (GAPDH) control for every reaction (Schmittgen and Livak, 2008). Results were expressed as mean value ± standard deviation (SD). More than 20 significantly differentially expressed genes were selected to validate the sequencing data and sequences of forward and reverse primers used in this study are listed in Table S1. The mRNA fold change was calculated by the following equations: ΔCT=ΔCT(target)−ΔCT(GAPDH); ΔΔCT = ΔCT(infected)−ΔCT(control); mRNA fold change = 2−ΔΔCT (Xu et al., 2012).
Illumina Sequencing and Read Mapping
We used RNA-sequencing (RNA-seq) of poly-adenylated RNA to capture the dynamic temporal changes of PK-15 cells and T. gondii RH strain transcriptomes at 1, 3, 6, 9, 12, 18, or 24 h PI. Nine sequencing libraries were prepared and sequenced with the Illumina paired-end method. More than 4.3 × 107 raw reads were generated in each group (Table 1). After removing low quality reads, clean reads were obtained and >95% of the clean reads had Phred-like quality scores at the Q20 level (an error probability of 0.01) and the GC-contents were about 50%.
We then mapped the sequenced clean reads to T. gondii and wild boar (Sus scrofa) genome, respectively. The majority of the clean reads were distributed in the exon region, followed by intergenic region and the intron region. Clean reads that were not mapped to the Sus scrofa genome were mapped to the genome sequence of T. gondii. About 90% of the total clean reads that were mapped to the Sus scrofa genome were unique. However, above 99% of the total clean reads were uniquely mapped to the T. gondii genome, while less than 1% of the clean reads were mapped to both reference genomes. As shown in Figure S1, the increased total reads and the unique mapped reads to T. gondii genome was proportional to time after infection. In contrast, total and uniquely mapped reads to the Sus scrofa genome decreased as time after infection increases. All subsequent analyses were based on the uniquely mapped reads.
Discovery and Characterization of SNPs and InDels in Transcriptomes
RNA-Seq is an efficient way to comprehensively identify single nucleotide polymorphisms (SNPs) and Insertions/Deletions (InDels) events from the expressed genes. Putative SNPs were predicted from T. gondii infected and control samples based on read depth and quality score of alignment results. Among the detected SNPs, A/G, G/A, C/T, and T/C transitions, were the most abundant and equally distributed (Figures 1A,B). For SNPs transversion, A/T and T/A were the smallest types. Alternative splicing (AS) is an important regulatory mechanism underlying increased cellular and functional complexity in eukaryotes and plays an important role in the generation of proteomic and biological diversities (Pan et al., 2008). Among the 12 basic types of AS events, transcription start site (TSS) and transcription terminal site (TTS) were the two major events (Figures 1C,D). Additionally, the number of putative SNPs and InDels identified from T. gondii transcriptomes increased proportional to time after infection. The AS events showed a similar trend.
Figure 1. Statistics of single nucleotide polymorphisms (SNPs) and alternative splicing (AS) events identified from the RNA-Seq data. (A) Classification of putative SNPs and INDELs identified from the transcriptomes of T. gondii. The x-axis represents types of SNP events. The y-axis represents number of SNP events. (B) Classification of putative SNPs and INDELs identified from the transcriptomes of PK-15 cells. The x-axis represents types of SNP events. The y-axis represents number of SNP events. (C) Numbers of the AS events in T. gondii. The x-axis represents types of AS events. The y-axis represents number of AS events. (D) Numbers of AS events in PK-15 cells. The x-axis represents types of AS events. The y-axis represents number of AS events. The AS events include transcription start site (TSS), transcription terminal site (TTS), skipped exon (SKIP_ON,SKIP_OFF pair) (SKIP), approximate SKIP (XSKIP_ON,XSKIP_OFF pair) (XSKIP), multi-exon SKIP (MSKIP_ON,MSKIP_OFF pair) (MSKIP), approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair) (XMSKIP), intron retention (IR_ON, IR_OFF pair) (IR), approximate IR (XIR_ON, XIR_OFF pair) (XIR), multi-IR (MIR_ON, MIR_OFF pair) (MIR), approximate MIR (XMIR_ON, XMIR_OFF pair) (XMIR), alternative exon ends (5′, 3′, or both) (AE), and approximate AE(XAE).
Dynamics of Gene Expression
Gene expression is calculated as FPKM and showed a comparable median distribution across different sequencing libraries (Figure S2). To characterize the expression patterns of mRNA activated during T. gondii infection, heat maps were constructed using FPKM values to determine overall transcriptome difference among different infection times PI. More than 80% of T. gondii genes were found to be expressed (FPKM>1). Heat map representation of differently expressed transcripts revealed a similarity in the trancriptomes within the early (T1, T3), the middle (T6, T9, T12) and the late (T18, T24) stages PI, respectively (Figure 2A). Parasites at stage T0 are just egressed from the host cells, which explains the clustering of the trascriptome at stage T0 with T18 and T24. Meanwhile, >50% of the PK-15 cell genes in the genome were found to be expressed (FPKM>1). Heat map based on the log2 FPKM values of differently expressed transcripts in PK-15 cells showed a similar cluster profile, except PK-15 cell transcriptome at T0, T1, and T3, which clustered together (Figure 2B). For differential gene expression analysis, pair wise comparisons of datasets from T. gondii during cell infection vs. egressed parasites (T0-T) and infected host cells vs. controls (T0-C) were performed, respectively. As shown in Figure 2C, more than half of the parasite's DEGs are down-regulated and the number of DEGs decreased in a time-dependent manner. However, the number of identified host DEGs fluctuated and peaked at T6, which was near the end of the first parasite cell cycle (Figure 2D). We further compared the DEGs in the parasites and host cells during the first T. gondii cell cycle (T1~T9), and the numbers of DEGs among these samples are visualized in the Venn diagrams. As shown in Figure 3A, DEGs belonging to a specific infection stage decrease in a time-dependent manner and there are 12 T. gondii DEGs that varied at all four infection stages. Meanwhile, 8 host DEGs varied during the process (Figure 3B). The validity of the differential expression of these DEGs that varied during the first four time points PI was confirmed by qRT-PCR analysis. The results of the qRT-PCR analysis agreed with the Illumina sequencing data, which are shown in Tables S2, S3.
Figure 2. Gene expression profiles during T. gondii infection. (A) Hierarchical clustering of differentially expressed genes of T. gondii during infection. Expression values (FPKM) were log2-transformed and subsequently median-centered by gene. Rows were hierarchically clustered based on average linkage using Pearson correlation coefficients as the distance measure. The expression levels are visualized using gradient color scheme, the scale from least abundant to highest range is from -2.0 to 2.0. Green color indicates low expression, and red color indicates high expression of the detected genes. The left vertical axis represents sample ID. The horizontal axis shows clusters of samples and the above vertical axis shows clusters of DEGs. (B) Hierarchical clustering of differentially expressed genes of T. gondii-infected PK-15 cells. (C) Numbers of DEGs in the parasites during infection process. (D) Numbers of DEGs in the host PK-15 cells during infection process. Up, up-regulated DEGs; Down, down-regulated DEGs; Total, total DEGs.
Figure 3. Venn diagram showing the overlap of DEGs clustered into four comparison groups (T1, T3, T6, and T9) represented by four ellipses. Numbers outside the ecllipses correspond to groups unique to a given infection stage. Numbers in the overlapping parts of different ellipses correspond to the number of DEGs in common from those comparison groups. (A) Veen diagram of identified DEGs in T. gondii. (B) Veen diagram of identified DEGs in PK-15 cells. Numbers in parentheses indicate the total number of DEGs in each data set.
Gene Ontology Enrichment Analysis
A total of 1436 T. gondii genes were considered to be significant when we compared T1 to T0-T. Of these, 524 genes were up-regulated and 912 were down-regulated. Gene GO is an international standardized gene functional classification system, which was applied to search for significantly enriched GO terms in DEGs. To obtain a comprehensive view of DEGs, we carried out the GO analysis using GOseq R package. Go term enrichment analysis was carried out to evaluate significantly over-represented GO terms of these 1436 DEGs and has resulted in a total of 2169 GO terms, including 1194 biological process terms, 304 cellular component terms and 671 molecular function terms. Each DEG could be assigned to more than one GO terms and the number of enriched GO terms was higher for down-regulated genes than for up-regulated genes. The top 40 enriched GO terms are shown in Figure 4A. The majority of the DEGs were classified under the “molecular function”. Within the three main categories of the GO classification, the “cellular process,” “cell,” and “binding” terms were most prevalent. We also noted that a high percentage of DEGs were classified under “metabolic process” and “catalytic activity,” while only a few DEGs were classified under “gene expression” and “biological regulation.” As infection progresses the number of DEGs decrease. A total of 288 DEGs were identified between T6 and T0-T, and these genes were categorized into 1100 GO terms. “Binding”, “membrane” and “cellular process” were the highest enrichment GO terms in the three main categories (Figure 4B).
Figure 4. DEGs GO cluster distribution. DEGs were classified into three main categories: biological process, cellular component and molecular function. Shown are the identified functions and the corresponding numbers of DEGs for each GO category. (A) T1/T0-T DEGs GO cluster distribution. (B) T6/T0-T DEGs GO cluster distribution. (C) T1/T0-C DEGs GO cluster distribution. (D) T6/T0-C DEGs GO cluster distribution. The y-axis indicates the number of genes in a specific category. X-axis indicates different GO terms. Top 40 enriched GO terms are listed.
Meanwhile, a total of 51 DEGs in the host cells were found between group T1 and T0-C. These genes were annotated to 522 GO terms, and “binding,” “intracellular” and “cellular processes” were the highest enrichment terms in the three main categories (Figure 4C). We then noted that in the host PK-15 cells the total number of DEGs between T6 and T0-C was 396 (205 up-regulated vs. 191 down-regulated). Go term enrichment analysis resulted in a total of 1491 GO terms, including 792 biological process terms, 240 cellular component terms and 459 molecular function terms. In contrast with the comparison between T1 and T0-C, the number of enriched GO terms in T6/T0-C was lower for down-regulated genes than for up-regulated genes, and the majority of the DEGs were classified under “molecular function,” especially under “binding” and “catalytic activity” (Figure 4D).
Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis
DEGs were mapped to the reference pathway in the KEGG database in order to identify the parasite biological pathway operating during infection process. KEGG pathways analysis of T. gondii DEGs showed that only small number of DEGs was annotated to the KEGG pathways and most of which were related to metabolism. During the intracellular infection cycle of T. gondii a total of 41 metabolic pathways were predicted (Table S4). The “pyrimidine metabolism,” “carbon metabolism,” “biosynthesis of secondary metabolites,” “purine metabolism,” and “biosynthesis of amino acids” were highly enriched. Of the 1436 DEGs between T1 and T0-T, 203 genes had a KO ID and were associated with 68 pathways (23 up-regulated genes categorized into 24 pathways, and 180 down-regulated genes categorized into 68 pathways). The top 20 enriched pathways were shown in Figure 5A and all the enriched pathway were closely related to metabolic activities except “sulfur relay system” and “base excision repair.” Additionally, “metabolic pathway” included 99 DEGs, account for half of those genes with KO ID. The cytoscape software was used to visualize the protein-protein interaction networks of the products of DEGs. As shown in Figure 6, TGME49_120110 (PCNA), TGME49_049180 (DHFR-TS), TGME49_055320 (hypothetical protein), and TGME49_002300 (ITPase) are the four hubs with most interactions.
Figure 5. Statistics of KEGG pathway enrichment. The x-axis shows the enrichment factor. The y-axis corresponds to KEGG Pathway. The color of the dot represents q value and size of the dot represents the number of DEGs mapped to the reference pathways. (A) Top 20 enriched pathways in T. gondii at T1. (B) Top 20 enriched pathways in PK-15 cell at T6.
Figure 6. Protein-protein interaction (PPI) networks of the DEGs identified in the comparison between T1 and T0-T. The PPI network with a high combined score >0.9 was prepared using the STRING 10 database program. DEGs are represented as round nodes. The red node indicates upregulated or green node indicates downregulated EDGs. The node size indicates high interaction degree (large) or low degree (small). Proteins that are associated to each other are linked by an edge. The color of the edge indicates the combined interaction score.
To further investigate the biochemical pathways of DEGs in the host cells during T. gondii infection, all host cell DEGs were mapped to terms in the KEGG database and compared with the whole transcriptome background. In order to get a better understanding of the defense system in host cells during infection we identified DEGs and KEGG pathways related to immune response or disease activities. As shown in Table 2, a total of 27 enriched immune-relevant pathways from the KEGG database were predicted. There was large number of DEGs at the infection stage T6 and T9, especially at stage T6. Of the 396 DEGs between T6 and T0-C, 157 genes had a KO ID and were associated with 194 pathways (87 up-regulated genes involved in 142 pathways, and 70 down-regulated genes involved in 134 pathways). The top 20 enriched pathways are shown in Figure 5B, with DEGs involved in disease-related “small cell lung cancer,” “renal cell carcinoma,” “prostate cancer,” “pancreatic cancer,” and “bladder cancer” being highly enriched. Additionally, among theses DEGs, some immune-related molecules involved in “TGF-β signaling pathway,” “p53 signaling pathway,” and “rap1 signaling pathway” were identified through KEGG enrichment. Additionally, variations of log2 fold changes of 95 genes that are involved in immune-related pathways are present in Figure 7A. The majority of these genes are active throughout the entire infection course except at stage T1. Meanwhile, the alterations of 99 genes that are involved in disease-related pathways are shown in Figure 7B, which shows a similar profile. To confirm the differential gene expression from the transcriptome data, real-time PCR was performed to validate the expression of the selected SLC2A2 and SERPINB2. As shown in Figures 7C,D, the expression patterns of these two DEGs were consistent with the high-throughput sequencing data, indicating that the transcriptome sequencing data were reliable.
Figure 7. Heat maps of the fold changes of immune (A) or disease (B) related genes with significantly changing expression levels across the entire infection course. DEGs have been divided into three parts: the left part indicates genes that are upregulated at most infection times; the middle part indicates genes that are upregulated or downregulated at all infection times; the right part indicates genes that are downregulated at most infection times. The range of expression is represented by a color grade ranging from low (blue) to high (yellow). The x-axis indicates DEGs and y-axis means different samples. (C) Relative abundance of SLC2A2 as determined by RT-PCR and HiSeq analysis. X-axis indicates different samples. Y-axis indicates log2 fold change values. (D) Relative abundance of SERPINB2 as determined by RT-PCR and HiSeq analysis. X-axis indicates different samples. Y-axis indicates log2 fold change values.
Despite being recognized as the most common source of food-borne toxoplasmosis, very little is known about the mechanisms of T. gondii infection in pigs. Transcriptomic studies can provide useful information on the underlying pathogenic mechanisms and interactions following the course of T. gondii infection. Indeed, gene expression patterns of T. gondii and host cells have already been investigated (He et al., 2016; Zhou et al., 2016). However, apart from one study that has taken a combined host and parasite transcriptional profiling approach (Melo et al., 2013), previous studies have focused on either the host cells or the parasite, which has largely caused knowledge gap in the host-pathogen interactome. To address this need and to elucidate the temporally dynamic gene expression patterns of both host cells and T. gondii, transcriptome sequencing and differential gene expression analysis of T. gondii infecting swine PK-15 cells were performed through the course of an in vitro infection.
Our combined transcriptome analysis has identified distinct mRNA profiles in host PK-15 cells in response to T. gondii infection and has elucidated the molecular pathways and gene transcripts involved in innate defenses of PK-15 cells infected with T. gondii. Parasite infection induced substantial changes in the host cell transcriptome and revealed time-specific transcripts. Despite continuous infection for 24 h, the abundance of gene transcripts exhibited a dramatic change within the first hour PI, which decreased as infection time increased (Figure 2). These findings are consistent with a previous study, which also reported a time-dependent reprogramming of host transcriptomes in response to T. gondii infection (Zhou et al., 2016).
To determine the DEGs profile of T. gondii, cDNA were prepared from samples collected at different time points PI and sequenced using Illumina technology. Hierarchical clustering to group genes according to their level of expression (log-transformed FPKM values) provided an overall view of the parasite transcriptomic landscape. This analysis revealed that early infection is associated with the most DEGs, which are more likely to play a role in cell invasion and PV formation. Surprisingly, early infection (T1/T0-C) did not induce significant changes in gene expression of host cells, indicating that T. gondii invasion occurs in a subtle manner to maneuver the host defense surveillance. However, host cell DEGs peaked at stage T6, a period near the end of the first parasite cell cycle.
The KEGG pathway-based analysis was performed to better understand the DEGs functions and interactions. Among the 68 KEGG pathways to which obtained DEGs between group T1 and T0-T were assigned, “metabolic pathways” represented the largest category, followed by “biosynthesis of secondary metabolites,” “pyrimidine metabolism,” and “carbon metabolism.” Of the 68 KEGG pathways, 41 metabolic pathways were identified, implying that T. gondii tachtzoites can affect the host by altering the host cell microenvironment for their own benefit. PPI analysis revealed that DHFR-TS is a hub gene, which plays an important role in parasite invasion, and it is recognized as a crucial drug target due to its involvement in the folate metabolism and nucleotide synthesis (Sharma et al., 2013). Another hub gene, PCNA, codes for an evolutionarily conserved protein in all eukaryotic species and plays key roles in DNA replication, chromatin remodeling, DNA repair, and cell cycle regulation (Strzalka and Ziemienowicz, 2011).
The activation of anti-T. gondii immune response requires precise regulation of host cellular mRNAs expression, and any aberrations in the host cellular mRNAs expression during T. gondii infection may impair fine mechanism controlling host anti-T. gondii immunity. Upon entry into host cell, T. gondii did not induce significant DEGs in host cells until stage T6. The DEGs identified were involved in essential biological pathways and differential expression analysis showed host immune genes, such as JUNB, CCL5, HSP90B1, FASN, TFDP1, MAP3K6, LAMC2, MALT1, RASSF5, ICAM-1, CCNE2, GADD45G, SFN, CISH_TV2, TGFB1, LAMA3, HSP70.2, GADD45B, ACKR3, and SRRM2 to be up-regulated throughout the infection cycle. In contrast, genes like MAP3K3, IGFBP3, ITGB5, IGF1R, CACNA2D1, GUCY1B3, PRKC1, GUCY1A3, CSF1, MCSF_ALPHA, EPHA1, MCL1, FGFR1IIC, ERBB3, NOD1, THBS1, PLCD1, and CDC25B were down-regulated.
We further analyzed the alterations of disease related genes in PK-15 cells. CYP1A1 and SERPINB2 displayed high expression levels along with T. gondii infection, and these two genes are involved in chemical carcinogenesis and amoebiasis, respectively (Faust and Guillen, 2012; Modi et al., 2012). Also, SLC2A2 showed decreased expression levels; the product of this gene plays an important role in glycogen storage disease (Santer et al., 2002). These results further substantiate the dynamic host-pathogen cross-talks that occur during T. gondii infection.
In summary, the present study provided the first combined RNA-seq-based transcriptomic analysis of T. gondii and host cells during the first 24 h of infection and the first in-depth analysis of the molecular mechanisms T. gondii uses to establish itself within porcine kidney cells. T. gondii infection induced significant changes in porcine PK-15 cell's transcriptome concurrent with considerable reprogramming of metabolic processes and signaling regulatory pathways in T. gondii. Immune- and disease-related DEGs in PK-15 cells were induced during the entire infection course (24 h), except at the beginning of infection, supporting the crucial role played by the host immune defenses against T. gondii infection. Parasite DEGs involved in metabolic pathways showed significant changes at the onset of the infection. These findings indicate the gene expression of PK-15 cells engaging with T. gondii is tightly linked. Future studies should employ a systems biology functional analysis pipeline to characterize not only the functions of the identified genes, but also the linked proteins and metabolites during T. gondii-host interactions. These would in turn shed light on key processes that underpin the inter-kingdom signaling communication between this highly zoonotic apicomplexan parasite and the mammalian host. Ultimately, this may lead to identification of new diagnostic biomarkers of T. gondii infection or novel therapeutic targets.
XZ, XS conceived and designed the experiments. CZ, QL performed the experiments. CZ analyzed the data. DZ contributed reagents/materials/analysis tools. CZ, HE wrote the paper. XZ contributed to the revision. All authors read and approved the final manuscript.
Data Availability Statement
The data sets supporting the results of this article are available in the NCBI SRA repository with the accession number SRP063957.
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.
Project support was provided by the National Natural Science Foundation of China (Grant No. 31230073) to XZ. We thank Novogene Bioinformatics Technology Co., Ltd (Beijing, China) for Illumina transcriptome sequencing and preliminary data analysis.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fmicb.2016.00721
Figure S1. Percent of reads uniquely mapped to reference genomes compared to total mapped. Tg reads that were mapped to T. gondii genome; C, reads that were mapped to pig genome. X axis, time post infection (h); y axis, percent of mapped reads.
Figure S2. The comparisons of gene expression levels among samples examined at different time points post infection (PI) and controls based on FPKM values. (A) Box plots compare the levels of actively expressed genes of parasites at different time points PI. (B) Box plots compare the levels of actively expressed genes of host PK-15 cells at different time points PI. X axis indicates different samples; y axis indicates log10 (FPKM+1) values in RNA-Seq libraries. T0-T indicates reads library from T. gondii control sample. T0-C indicates read library from PK-15 cell control sample.
Table S1. List of primers used in real-time quantitative RT-PCR analysis.
Table S2. Real-time qRT-PCR of DEGs in the parasites that vary at all four infection stages (T1, T3, T6, and T9).
Table S3. Real-time qRT-PCR of DEGs in the host PK-15 cells that vary at all four infection stages (T1, T3, T6, and T9).
Table S4. Statistics of the number of parasite DEGs involved in metabolic activities.
ACKR3, atypical chemokine receptor 3; CACNA2D1, calcium voltage-gated channel auxiliary subunit alpha 2 delta 1; CCL5, C-C motif chemokine ligand 5; CCNE2,cyclin E2; CDC25B, cell division cycle 25B; CISH_TV2, cytokine inducible SH2-containing protein; CSF1, colony stimulating factor 1; CYP1A1, cytochrome P450 family 1 subfamily A member 1; DHFRTS, bifunctional dihydrofolate reductase-thymidylate synthase; EPHA1, EPH receptor A1; ERBB3, erb-b2 receptor tyrosine kinase 3; FASN, fatty acid synthase; FGFR1IIC, fibroblast growth factor receptor 1; GADD45B, growth arrest and DNA damage inducible beta; GADD45G, growth arrest and DNA damage inducible gamma; GUCY1A3, guanylate cyclase 1,soluble, alpha 3;GUCY1B3, guanylate cyclase 1,soluble, beta 3; HSP70.2, heat shock 70 kDa protein 1B;HSP90B1, heat shock protein 90kDa beta family member 1; ICAM-1, intercellular adhesion molecule 1; IGFBP3, insulin like growth factor binding protein 3; IGF1R, insulin like growth factor 1 receptor; ITGB5, integrin subunit beta 5; ITPase, inosine triphosphate pyrophosphatase; JUNB, jun B proto-oncogene; LAMA3, laminin subunit alpha 3; LAMC2, laminin subunit gamma 2; MALT1, mucosa-associated lymphoid tissue lymphoma translocation protein 1; MAP3K3, mitogen-activated protein kinase kinase kinase 3; MAP3K6, mitogen-activated protein kinase kinase kinase 6; MCL1, myeloid cell leukemia 1; MCSF_ALPHA, macrophage colony stimulating factor alpha; NOD1, nucleotide binding oligomerization domain containing 1; PCNA, proliferating cell nuclear antigen; PLCD1, phospholipase C delta 1; PRKC1, protein kinase C1; RASSF5, Ras association domain family member 5; SERPINB2, serpin family B member 2; SFN, stratifin; SLC2A2, solute carrier family 2 member 2; SRRM2, serine/arginine repetitive matrix 2; TFDP1, transcription factor Dp-1; TGFB1, transforming growth factor beta 1 and THBS1, thrombospondin 1.
Behnke, M. S., Wootton, J. C., Lehmann, M. M., Radke, J. B., Lucas, O., Nawas, J., et al. (2010). Coordinated progression through two subtranscriptomes underlies the tachyzoite cycle of Toxoplasma gondii. PLoS ONE 5:e12354. doi: 10.1371/journal.pone.0012354
Bela, S. R., Dutra, M. S., Mui, E., Montpetit, A., Oliveira, F. S., Oliveira, S. C., et al. (2012). Impaired innate immunity in mice deficient in interleukin-1 receptor-associated kinase 4 leads to defective type 1 T cell responses, B cell expansion, and enhanced susceptibility to infection with Toxoplasma gondii. Infect. Immun. 80, 4298–4308. doi: 10.1128/IAI.00328-12
Butler, C. L., Lucas, O., Wuchty, S., Xue, B., Uversky, V. N., and White, M. (2014). Identifying novel cell cycle proteins in apicomplexa parasites through co-expression decision analysis. PLoS ONE 9:e97625. doi: 10.1371/journal.pone.0097625
DePristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J. R., Hartl, C., et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. 43, 491–498. doi: 10.1038/ng.806
Florea, L., Song, L., and Salzberg, S. L. (2013). Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res. 2:188. doi: 10.12688/f1000research.2-188.v2
Gaji, R. Y., Behnke, M. S., Lehmann, M. M., White, M. W., and Carruthers, V. B. (2011). Cell cycle-dependent, intercellular transmission of Toxoplasma gondii is accompanied by marked changes in parasite gene expression. Mol. Microbiol. 79, 192–204. doi: 10.1111/j.1365-2958.2010.07441.x
Gazzinelli, R. T., Mendonça-Neto, R., Lilue, J., Howard, J., and Sher, A. (2014). Innate resistance against Toxoplasma gondii: an evolutionary tale of mice, cats, and men. Cell Host Microbe 15, 132–138. doi: 10.1016/j.chom.2014.01.004
He, J. J., Ma, J., Song, H. Q., Zhou, D. H., Wang, J. L., Huang, S. Y., et al. (2016). Transcriptomic analysis of global changes in cytokine expression in mouse spleens following acute Toxoplasma gondii infection. Parasitol. Res. 115, 703–712. doi: 10.1007/s00436-015-4792-5
Hu, K., Mann, T., Striepen, B., Beckers, C. J., Roos, D. S., and Murray, J. M. (2002). Daughter cell assembly in the protozoan parasiteToxoplasma gondii. Mol. Biol. Cell 13, 593–606. doi: 10.1091/mbc.01-06-0309
Khan, A., Dubey, J. P., Su, C., Ajioka, J. W., Rosenthal, B. M., and Sibley, L. D. (2011). Genetic analyses of atypical Toxoplasma gondii strains reveal a fourth clonal lineage in North America. Int. J. Parasitol. 41, 645–655. doi: 10.1016/j.ijpara.2011.01.005
Liu, J., Wang, S., Qin, T., Li, N., Niu, Y., Li, D., et al. (2015b). Whole transcriptome analysis of Penicillium digitatum strains treatmented with prochloraz reveals their drug-resistant mechanisms. BMC Genomics 16, 855. doi: 10.1186/s12864-015-2043-x
McLeod, R., Boyer, K. M., Lee, D., Mui, E., Wroblewski, K., Karrison, T., et al. (2012). Prematurity and severity are associated with Toxoplasma gondii alleles (NCCCTS, 1981-2009). Clin. Infect. Dis. 54, 1595–1605. doi: 10.1093/cid/cis258
Melo, M. B., Nguyen, Q. P., Cordeiro, C., Hassan, M. A., Yang, N., McKell, R., et al. (2013). Transcriptional analysis of murine macrophages infected with different Toxoplasma strains identifies novel regulation of host signaling pathways. PLoS Pathog. 9:e1003779. doi: 10.1371/journal.ppat.1003779
Minot, S., Melo, M. B., Li, F., Lu, D., Niedelman, W., Levine, S. S., et al. (2012). Admixture and recombination among Toxoplasma gondii lineages explain global genome diversity. Proc. Natl. Acad. Sci. U.S.A. 109, 13458–13463. doi: 10.1073/pnas.1117047109
Modi, B. G., Neustadter, J., Binda, E., Lewis, J., Filler, R. B., Roberts, S. J., et al. (2012). Langerhans cells facilitate epithelial DNA damage and squamous cell carcinoma. Science 335, 104–108. doi: 10.1126/science.1211600
Nelson, M. M., Jones, A. R., Carmen, J. C., Sinai, A. P., Burchmore, R., and Wastling, J. M. (2008). Modulation of the host cell proteome by the intracellular apicomplexan parasite Toxoplasma gondii. Infect. Immun. 76, 828–844. doi: 10.1128/IAI.01115-07
Pan, Q., Shai, O., Lee, L. J., Frey, B. J., and Blencowe, B. J. (2008). Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing. Nat. Genet. 40, 1413–1415. doi: 10.1038/ng.259
Pena, H. F., Gennari, S. M., Dubey, J. P., and Su, C. (2008). Population structure and mouse-virulence of Toxoplasma gondii in Brazil. Int. J. Parasitol. 38, 561–569. doi: 10.1016/j.ijpara.2007.09.004
Pissuwan, D., Valenzuela, S. M., Miller, C. M., and Cortie, M. B. (2007). A golden bullet? Selective targeting of Toxoplasma gondii tachyzoites using antibody-functionalized gold nanorods. Nano Lett. 7, 3808–3812. doi: 10.1021/nl072377
Qu, D., Zhou, H., Han, J., Tao, S., Zheng, B., Chi, N., et al. (2013). Development of reverse transcription loop-mediated isothermal amplification (RT-LAMP) as a diagnostic tool of Toxoplasma gondii in pork. Vet. Parasitol. 192, 98–103. doi: 10.1016/j.vetpar.2012.10.010
Radke, J. R., Donald, R. G., Eibs, A., Jerome, M. E., Behnke, M. S., Liberator, P., et al. (2006). Changes in the expression of human cell division autoantigen-1 influence Toxoplasma gondii growth and development. PLoS Pathog. 2:e105. doi: 10.1371/journal.ppat.0020105
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Sansom, S. N., Shikama-Dorn, N., Zhanybekova, S., Nusspaumer, G., Macaulay, I. C., Deadman, M. E., et al. (2014). Population and single-cell genomics reveal the Aire dependency, relief from Polycomb silencing, and distribution of self-antigen expression in thymic epithelia. Genome Res. 24, 1918–1931. doi: 10.1101/gr.171645.113
Santer, R., Groth, S., Kinner, M., Dombrowski, A., Berry, G. T., Brodehl, J., et al. (2002). The mutation spectrum of the facilitative glucose transporter gene SLC2A2 (GLUT2) in patients with Fanconi-Bickel syndrome. Hum. Genet. 110, 21–29. doi: 10.1007/s00439-001-0638-6
Sharma, H., Landau, M. J., Vargo, M. A., Spasov, K. A., and Anderson, K. S. (2013). First three-dimensional structure of Toxoplasma gondii thymidylate synthase-dihydrofolate reductase: insights for catalysis, interdomain interactions, and substrate channeling. Biochemistry 52, 7305–7317. doi: 10.1021/bi400576t
Steinparzer, R., Reisp, K., Grünberger, B., Köfer, J., Schmoll, F., and Sattler, T. (2015). Comparison of different commercial serological tests for the detection of Toxoplasma gondii antibodies in serum of naturally exposed pigs. Zoonoses Public Health 62, 119–124. doi: 10.1111/zph.12122
Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi: 10.1093/nar/gku1003
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Tymoshenko, S., Oppenheim, R. D., Agren, R., Nielsen, J., Soldati-Favre, D., and Hatzimanikatis, V. (2015). Metabolic needs and capabilities of Toxoplasma gondii through combined computational and experimental analysis. PLoS Comput. Biol. 11:e1004261. doi: 10.1371/journal.pcbi.1004261
Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483
Xu, X., Liu, T., Zhang, A., Huo, X., Luo, Q., Chen, Z., et al. (2012). Reactive oxygen species-triggered trophoblast apoptosis is initiated by endoplasmic reticulum stress via activation of caspase-12, CHOP, and the JNK pathway in Toxoplasma gondii infection in mice. Infect. Immun. 80, 2121–2132. doi: 10.1128/IAI.06295-11
Zhang, N., Zhang, H. J., Zhao, B., Sun, Q. Q., Cao, Y. Y., Li, R., et al. (2014). The RNA-seq approach to discriminate gene expression profiles in response to melatonin on cucumber lateral root formation. J. Pineal. Res. 56, 39–50. doi: 10.1111/jpi.12095
Zhou, C. X., Zhou, D. H., Liu, G. X., Suo, X., and Zhu, X. Q. (2016). Transcriptomic analysis of porcine PBMCs infected with Toxoplasma gondii RH strain. Acta Trop. 154, 82–88. doi: 10.1016/j.actatropica.2015.11.009
Zhou, D. H., Yuan, Z. G., Zhao, F. R., Li, H. L., Zhou, Y., Lin, R. Q., et al. (2011). Modulation of mouse macrophage proteome induced by Toxoplasma gondii tachyzoites in vivo. Parasitol. Res. 109, 1637–1646. doi: 10.1007/s00436-011-2435-z
Keywords: Toxoplasma gondii, PK-15 cells, host-pathogen interaction, transcriptome, KEGG, protein-protein interaction
Citation: Zhou C-X, Elsheikha HM, Zhou D-H, Liu Q, Zhu X-Q and Suo X (2016) Dual Identification and Analysis of Differentially Expressed Transcripts of Porcine PK-15 Cells and Toxoplasma gondii during in vitro Infection. Front. Microbiol. 7:721. doi: 10.3389/fmicb.2016.00721
Received: 03 February 2016; Accepted: 29 April 2016;
Published: 13 May 2016.
Edited by:Jörg Linde, Leibniz-Institute for Natural Product Research and Infection Biology -Hans-Knoell-Institute, Germany
Reviewed by:Biswapriya Biswavas Misra, University of Florida, USA
Samer Angelone-Alasaad, Zurich University, Switzerland
Copyright © 2016 Zhou, Elsheikha, Zhou, Liu, Zhu and Suo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.