Identification of putative cis-regulatory elements in Cryptosporidium parvum by de novo pattern finding

Genes with conserved upstream motifs have similar expression profiles

The groups of genes selected in this study are involved in parasite-specific functions as well as housekeeping-type activities (Table 1). Parasite-specific gene families included genes encoding cryptosporidial oocyst wall proteins (COWPs) and large secretory proteins (Cp LSPs), genes known to show concerted post-infection expression patterns [4, 11]. The housekeeping genes used in the analysis included genes involved in nucleotide salvage, DNA replication and glycolysis. We show that each of these groups of genes share different conserved upstream motifs. No common, general motif was conserved across all groups, barring AT-rich stretches, which were not statistically significant given the AT-richness of the genome.

Algorithms for pattern finding can report several hundred motifs ordered by their statistical significance as compared to a background model. Since statistical significance alone is not a sufficient indicator of biological significance, we applied a rule-based approach to identify candidate motifs that warrant further investigation. Candidate motifs were required to be within the top ten motifs predicted by both algorithms; display high information content and preferably show multiple occurrences within each upstream region. Information content as determined by MEME depends on the frequencies of the bases in a given column compared to the overall frequencies of those bases in the group of sequences. The more conserved the position is and the more rare the conserved nucleotides are, the higher the information content is. We provided MEME with a background file containing all inter-coding regions in the genome against which information content was calculated for each motif. (See additional files 2 and 3 for the scores, E-values and information content of reported motifs).

Each group of genes that shared conserved, upstream motifs was examined for correlated expression profiles via comparative real-time PCR at 6 different post-infection time points. The housekeeping groups of genes were further resolved into three sub-sets based on their expression profiles. For each sub-set, a second iteration of pattern finding was performed to determine if conserved motifs within each sub-set existed. We find a correlation between genes that contain distinct, conserved, upstream motifs and their corresponding expression profiles over a 72 h post-infection period.

Genes encoding the COWP family

The cryptosporidial oocyst wall proteins comprise a multi-gene family that demonstrate a defined pattern of expression during in vitro development, with expression levels peaking at 48 h through 72 h post-infection [11]. Genes encoding members of this family are scattered throughout the genome and not clustered in a tandem array. The most significant motif found in the upstream regions of these genes as determined by both algorithms was a 12 bp motif (Figure 2a). This motif is present one or two times in the upstream region of all COWP genes, and when present in pairs, the motifs are often within 50 – 100 bp of each other. The promoter regions of this gene family are not alignable outside of the conserved motifs identified, indicating that the conserved motif is not simply a function of recent gene & promoter duplication.

Genes encoding the Cp LSP family

The large secretory proteins comprise a gene family that shows genomic co-localization in a cluster on chromosome 7. They are also co-expressed during the life cycle [4]. Figure 2b shows a single DNA motif found upstream of each of these genes, with a well-conserved sequence. This motif occurs 2–3 times in all of the upstream regions, and is often located within -350 bp from the translational start. As is the case with the COWP gene family, the promoter regions of this gene family are also not alignable outside of the conserved motifs identified, indicating that the conserved motif is not simply a function of recent gene & promoter duplication.

Genes involved in nucleotide metabolism

Cryptosporidium parvum possesses highly streamlined nucleotide metabolic pathways, relying on the host cell for the salvage of both purine and pyrimidine residues. These pathways also contain genes that have been transferred into the nuclear genome of C. parvum from bacteria and plants via intracellular or horizontal gene transfer. The essential functions of these genes and their distinct evolutionary origin make them important drug-targets in developing anti-cryptosporidial chemotherapy [12]. We examined ten genes involved in nucleotide salvage and modification to identify significant motifs common to their upstream regions. We could not find a significant motif reported by both algorithms to be present upstream of all of the genes. MEME alone reported an 8 bp AT-rich motif present at least once in all the sequences at varying positions from the translational start (Figure 3a). Based on their real-time PCR expression profiles, the genes were divided into 3 sub-sets (sub-set 1, 2 and 3). Sub-set 1 contained three enzymes involved in the transport and modification of purines (AT, IMPDH and GMPS), and also one pyrimidine-modifying enzyme (CTPS). These genes were characterized by high expression levels at 2 h and 12 h post-infection and the most significant motif specific to this sub-set was a 10 bp motif shown (Figure 3a, sub-set 1). This motif was often found multiple times in the upstream regions and almost always found on the reverse strand. Three remaining pyrimidine-modifying enzymes (RDPR, dCMPD and DHFR-TS) had expression levels that peaked at 48 h post-infection and dropped subsequently. They comprise sub-set 2. These genes contained a 12 bp motif in their upstream regions, also seen at varying positions from the translational start (Figure 3a, sub-set 2). The three kinases (AK, TK and UK) involved in nucleotide salvage were grouped together in sub-set 3 based on their high expression levels at 48 h and 72 h post-infection. They were found to contain a conserved 14-bp AT-rich motif in their upstream regions (Figure 3a, sub-set 3).

Genes involved in DNA replication

Analysis of the C. parvum genome reveals that the organism possesses a reduced complement of genes involved in DNA replication [4]. We chose to study genes involved in DNA replication expecting that they would be co-regulated in a time-dependent manner associated with the life cycle [13, 14]. The most significant motif identified by both MEME and AlignACE was a single G-rich motif present upstream of all of the genes occurring multiple times in some of the upstream regions (Figure 3b). These genes could be resolved into 3 sub-sets based on their comparative RT-PCR expression profiles. Three genes peaking at 2 h post-infection were classified into sub-set 1. A 14 bp motif with a core conserved 5'-CGCCAA-3' sequence was found occurring once upstream of these three genes (Figure 3b, sub-set 1). At 6 h post-infection, a few genes coding for MCM-like proteins and the single-stranded binding protein RP-A were found to peak in expression levels. These were classified into sub-set 2. The most significant motif found specific to this sub-set was an 11 bp motif occurring one or two times in the upstream regions of this sub-set (Figure 3b, sub-set 2). Most of the MCM-like proteins peaked at 48 h post-infection and were classified into sub-set 3. These genes were found to contain a 13 bp motif, with a relatively less-conserved sequence, present one or more times in their upstream regions (Figure 3b, sub-set 3).

Genes involved in glycolysis

Glycolysis is considered to be the main source of energy in the Coccidia, and especially so in C. parvum due to the lack of evidence of a mitochondrion and a functional respiratory chain [4, 15, 16]. Ten genes associated with glycolysis were considered for this study, and a single motif was found to be over-represented in all of their upstream regions (Figure 3c). This 10 bp motif contains a core 5'-GGCG-3' sequence and is present multiple times in some of the upstream regions. No outstanding pattern with respect to the orientation or position relative to the start of translation is apparent. Comparative real-time PCR experiments resolved the glycolytic genes into three sub-sets based on their expression profiles during development. The genes that peaked at 6 h- 12 h post-infection were included in sub-set 1. The most significant motif found upstream of these genes was a 14 bp motif occurring uniquely, or many times and almost always on the opposite strand (Figure 3c, sub-set 1). Sub-set 2 is comprised of 3 genes exhibiting weakly correlated expression profiles, two of which peak at 48 h post-infection. An 11 bp C-rich motif was found to occur one time in each of their upstream regions (Figure 3c, sub-set 2). We are encouraged that a conserved motif was found for this group, however, we remain unconvinced about the validity of this group based on their expression profiles alone since there is a discrepancy in the profiles. Sub-set 3 consists of 3 genes with expression levels peaking at 72 h post-infection. A 9 bp motif was found to be common to their upstream regions. This motif was very well conserved at the sequence level (TGC [A/G] [T/G]G [C/G]GA) (Figure 3c, sub-set 3) and was found occurring once upstream each gene.

Genome-wide occurrences of candidate motifs

The candidate motifs reported in this study were selected based on rules as described earlier. In all the cases except one, these motifs are also found elsewhere in the genome, upstream of other genes. This is not surprising given that these are short sequences, with several degenerate positions. Expression profiles are not yet available for every gene to test if the other genes containing these motifs upstream also display similar expression profiles. (See additional file 3).

Comparative studies in C. hominis and other apicomplexans

Comparative analyses of upstream regions in Toxoplasma gondii, a more distant apicomplexan does not reveal the presence of the same conserved motifs in the groups studied (data not shown). This is not surprising considering the evolutionary distance between these species, and indicates that other apicomplexans may not serve as appropriate model systems for exploring the role of these C. parvum motifs.

Gene prediction and retrieval of intergenic regions

We used GLIMMER [30] to predict genes on the C. parvum genome, wrote scripts in PERL to extract gene-coordinate information and created the intergenic regions database (3396 sequences). To exclude the possibility of including coding regions in this set, a BLASTX was performed against known annotated C. parvum proteins using the set of intergenic regions as the query. 1000 sequences that contained portions of 100% identity to fragments of C. parvum protein sequences were trimmed to remove the protein coding regions. The upstream region for a gene refers to the entire intergenic region upstream of the translational start.

Organization of genes into functional groups

In the absence of expression information, we classified genes into putatively co-regulated groups based on their function. To identify genes belonging to each pathway/group, we made use of existing annotation, and BLASTP searches using orthologues from the related apicomplexans Plasmodium falciparum and Toxoplasma gondii.

Identification of conserved motifs in the upstream regions

We applied two pattern finding algorithms MEME and AlignACE to identify de novo patterns in the upstream regions. We used a background model based on the entire set of intergenic regions (3396 sequences) in C. parvum to train these algorithms. To identify patterns, the length range was set between 8 to 20 bp, and three different modes of occurrence were specified. The top 10 non-overlapping results from each algorithm were examined and compared, and the best motifs predicted by both algorithms were selected. We used WebLogo [31] to create sequence logos to represent the best motifs found in each search.

C. parvum culture and RNA isolation

C. Parvum infected cultures – Human ileocecal adenocarcinoma cells (HCT-8, ATCC CCL-244; American Type Culture Collection, Rockville, MD.) were plated on 10 cm plates and cultured to approximately 70% confluency in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS), sodium pyruvate, and antibiotics/antimycotic solution (100U penicillin G/ml, 100 μg streptomycin/ml and 0.25 μg amphotericin B/ml). C. parvum oocysts (Iowa strain) harvested from calves were purchased commercially (Pleasant Hill farms), stored at 4°C and used for in vitro infections prior to three months of age as previously described [11]. Briefly, oocysts were surface sterilized by treatment with a 1:3 dilution of Clorox bleach (1 ml/3 × 107oocysts) on ice for 7 minutes, washed repeatedly with Hank's buffered saline solution (HBSS), and stored in Cp media [RPMI 1640 media containing 10% fetal bovine serum, 15 mM Hepes, 50 mM glucose, 0.1u bovine insulin/ml, 35 μg ascorbic acid/ml, 1.0 μg folic acid/ml, 4.0 μg 4-aminobenzoic acid/ml, 2.0 μg calcium pantothenate/ml, 100U of penicillin G/ml, 100 μg of streptomycin/ml and 0.25 μg of amphotericin B/ml (pH 7.4)] at 4°C overnight. HCT-8 cultures were switched to Cp media approximately 18 hours prior to infection. Oocysts were warmed to room temperature for 30 minutes, and added to HCT-8 monolayers at a 1:1 ratio. Cells were incubated in a humidified incubator at 37°C in an atmosphere containing 5% CO₂. Following a 2 h excystation period at 37°C, cells were washed repeatedly with warm HBSS and incubated at 37°C in fresh Cp media until harvested. Infection was estimated to be between 70%-90% depending on the batch and storage period of oocysts. Total RNA was harvested in TRIzol reagent (Invitrogen) at 2, 6, 12, 24, 48, and 72 hours post infection and purified by following manufacture's instructions. Mock-infected cultures, cultures treated identically with the exception of infection, were harvested at exact time points as C. parvum infected cultures. Three independent time-courses were plated, infected, and harvested for this study.

Comparative real-time PCR

To investigate gene expression during C. parvum in vitro development, gene-specific primers (see additional file 1 for primer sequences) were designed and used in comparative real-time PCR analysis. First strand cDNA was made using manufacturer (Invitrogen) protocols. Briefly, 2 μg of total RNA that had been previously DNased following manufacturer instructions (Turbo DNase, Ambion) was mixed with 0.5 μg of random hexamer and RNase-free water (to bring up the volume to 12 μl), heated at 70°C for 10 min, and cooled on ice. To this mixture was added a 7 μl aliquot, consisting of 4 μl of 5X first strand buffer, 2 μl of 100 mM dithiothreitol, and 1 μl of 10 mM dNTP mixture. The reaction was equilibrated at 42°C for 2 min on an iCycler (BioRad), after which 1 μl (200 U) of SuperScript II RT was applied. The reaction mixture was then incubated at 42°C for 50 min, heated at 70°C for 15 min, and held at 4°C. Identical reactions were set up without the addition of reverse transcriptase to test for the presence of contaminating genomic DNA using primers to C. parvum rRNA and 50 cycles of PCR under standard conditions. cDNA made from RNA containing no detectable product in the above PCR reaction was used for comparative real-time PCR.

Comparative real-time PCR was performed using a Stratagene Mx3000 P real-time instrument in a 96 well format. Due to the sensitivity of the machine the amount of ROX normalizing dye needed to be reduced in the reactions. Therefore, 20 μl reactions were set-up using a modified master mix consisting of 10 μl 2X SYBR master mix (1 part SYBR green master mix containing ROX dye to 5 parts SYBR green master mix without ROX dye), 2 η mol of each primer and water to 13 μl. cDNA was diluted 1:150 and 7 μl of template was added to each reaction. After an initial denaturation at 94°C for 2 min, the reaction mixture underwent 42 cycles of denaturation at 94°C for 30 sec, annealing at 58 or 59°C for 20 sec, and extension at 68°C for 30 sec. Fluorescence was read after the end of each annealing cycle. Following the end of amplification cycles, a melting curve was run. This cycle started with an initial denaturation step at 94°C followed by annealing at 56°C. Melting was performed by increasing the temperature in single degree increments until the temperature reached 94°C. Fluorescence of SYBR green was read at each increase and the data was plotted onto a graph using the Mx3000 P software. cDNA made from RNA harvested from three independent timecourses was run in duplicate reactions and the average Ct value of each duplicate was determined using the Stratagene Mx3000 P software. As the number of developing C. parvum life stages within infected cells changes over time, primers specific for C. parvum 18S rRNA were used to normalize the amount of cDNA product of the target genes to that of C. parvum rRNA in the same sample. Due to the fact that rRNA is much more abundant than any specific mRNAs, the cDNA was diluted an additional 40 times and reactions were set up as above using three replicates of each time-course. Average C. parvum rRNA Ct values for each time-course was determined as above. An additional single time point of each time-course was run using C. parvum rRNA primers on each plate containing target genes to test the consistency of runs from plate to plate.

To determine comparative gene expression of a target gene, the average Ct values of all time points (expressed in log scale) of a single biological replicate were linearized and the ratio of expression for each time point was determined by dividing each by the product of the time point with the lowest expression (highest Ct). The ratio of expression for each time point of C. parvum rRNA expression was determined exactly as above. Relative target gene expression was determined by normalizing the ratio of the target gene expression to the ratio of C. parvum rRNA expression. Values from three biological replicates were imported into Prism (GraphPad Software). The mean and standard errors for each time point was determined using Prism's statistical package and the resulting graph was normalized to 100% maximum expression.