Transcriptional pathways associated with the slow growth phenotype of transformed Anaplasma marginale

Insert analysis

The OPTIMIZER web server was used to perform predictions and adjustments in codon usage for the insert present in AmTR [15]. The ‘one amino acid–one codon’ method, a guided random method based on a Monte Carlo algorithm, was used. The A. marginale str. St. Maries codon usage table was extracted from Artemis [16].

Cell culture and RNA-seq

The accession number for this RNA-seq study is: SRP014580. Three T75 flasks of ISE6 cells were infected with either wild type A. marginale or AmTR, representing three biological replicates. Both inoculums were originally derived from the St. Maries strain of A. marginale. The cell lines were maintained in ISE6 cells cultured at 34°C as previously described [9, 17, 18]. When passage 27 of either wild type or AmTR in a T75 cell culture flask infected greater than 80% of ISE6 cells, as determined by examination of Giemsa stained cytospin preparations, cell cultures were mechanically lysed by passaging 50 times through a bent 27 g needle attached to a 6 mL luer lock syringe, then filtered through a 2 μm syringe filter (Whatman) which allows passage of A. marginale but not cellular debris. The filtrate was pelleted by centrifugation, and then washed in PBS. The final pellet was resuspended in 200 μL PBS before adding 1 mL of TRIzol (Invitrogen). Total RNA was isolated from wild type or AmTR-infected ISE6 cell culture using TRIzol (Invitrogen), per manufacturer directions. Eukaryotic sequences were negatively selected through hybridization using the MICROBEnrich kit (Ambion). The Duplex‒Specific thermostable Nuclease (DSN) normalization protocol was applied to all samples. Samples were sequenced with Illumina technology with 100 bp reads, each sample was run on 1/3 of a lane. Data was processed using CLC Genomics Workbench (CLC Bio). Mapping parameters were adjusted to map a maximum number of reads to the reference bacterial genomes. The distribution of the expression values for all samples was analyzed and compared. Normalization by quantiles was applied to adjust the distributions for further comparison. Fold changes with respect to RPKM values were calculated [19]. Two different tests were applied to evaluate the statistical significance of fold changes: Kal and Baggerly’s statistical tests on proportions [20, 21].

Comparative transcriptional analysis

Expression was measured using quantitative reverse transcription PCR (qPCR) using the SYBR Green ER RT-PCR Kit (Invitrogen). Briefly, 200 ng of RNA were processed with the SuperScript® VILO™ cDNA Synthesis Kit (Invitrogen) to obtain cDNA. The steady state, single copy gene msp5 was used to calibrate the qPCR. Relative expression ratios were calculated by a mathematical model, which includes efficiency correction of individual transcripts through the REST software [22] with modifications as previously described [23]. This procedure was used for verification of transcriptional differences found by RNA-seq. The variance of the normalized expression values of the genes flanking the insert was evaluated through ANOVA paired with an F-test.

K-means clustering

We used K-means clustering to assign the mean transcriptional value of each gene to the cluster whose center is nearest. Lloyd’s algorithm was used for these experiments [24]. Euclidean distance was used as distance metric; five partitions were used to generate the clusters. For each gene, the mean gene expression value over all input samples was subtracted. Normalized expression values were used for clustering.

Gene Set enrichment analysis

The Gene Set Enrichment Analysis (GSEA) test described by Tian et al. [25] was used. Briefly, the test implicitly calculates and uses ANOVA statistic for multiple group experiments for each feature, as measures of association. For each category, the test statistics for the features in that category are summed and a category based test statistic is calculated as this sum divided by the square root of the number of features in the category. In order to annotate the A. marginale genome with Gene Ontology we used the Comprehensive Microbial Resource [26]. Briefly, we matched feature IDs in the A. marginale genome to synonyms and gene product form IDs from the Comprehensive Microbial Resource. Additionally, we complemented this source with manual curation of the annotation file, including genes reported in the Kyoto Encyclopedia of Gene and Genomes (KEGG) [27].

Insert characteristics

The nature of the construct used for transformation can be a source for phenotypic alterations. The effect of insert size on bacterial fitness (measured by growth rate), has been previously shown for bacteria with intermediate genome sizes [11, 12]. Inserts that ranged from 2.9 to 12.6 kb were used to transform either E. coli or Bacillus subtilis. In both studies, bacteria transformed using the largest inserts displayed slower growth rates. These larger inserts represent a small proportion of the chromosome they are inserted in: 0.2 and 0.3% of the total genome size in E. coli and B. subitlis, respectively. The insert used for transformation of A. marginale represents 0.4% of its single chromosome. Obligate intracellular pathogens in the Rickettsia and Coxiella genera have been transformed with inserts that vary in size, ranging from 0.001 to 0.5% of their chromosomes [28–30]. As these transformants displayed different growth phenotypes, an association between insert size and slow growth seems unlikely.

Foreign genes and their unfamiliar codon usage can affect bacteria through slow translation efficiency, difficultly folding the foreign protein, and protein degradation [31, 32]. The issue of codon usage is highlighted for bacteria with unusual genomic G + C content, such as the obligate intracellular pathogens in the order Rickettsiales [33, 34]. In order to analyze the relationship of the insert codon usage with that of A. marginale, we used OPTIMIZER, a web server that optimizes a DNA sequence using pre-computed codon usage tables [15]. Comparison of the insert with the recipient A. marginale revealed a Codon Adaptation Index (CAI) of 0.802. CAI values range from 0, when the codon usage of a sequence and that of the reference set are very different, to 1 when both codon usages are the same. These results suggest that codon usage is an unlikely explanation for reduced growth rate in AmTR.

Whole genome transcriptional profiling

Since the genes immediately adjacent to the insertion site were not significantly altered, we next identified the most differentially transcribed genes across all replicates. The top ten down and up regulated genes are shown in Figure 1. These genes were all significantly differentially transcribed (p < 1E-05) (Additional file 1 & Figure 1). Four of the top differentially transcribed genes (AM360, AM357, rnhB and fpba), one of the genes flanking the insertion site (ndk) and a gene involved in translation (rpmE) were evaluated with qPCR to validate the RNA-seq data (Additional file 2).

The RNA-seq analysis revealed genes that were significantly differentially transcribed, including genes with predicted functions, such as rnhB, fpba, tRNA-Phe-1, xerD, guaB, purL, dnaE, tRNA-Arg-3 and tRNA-Ser-3. Gene rnhB, which encodes RNase H, an enzyme that specifically cleaves the RNA strand of RNA/DNA hybrids and plays a role in removing RNA primers from Okazaki fragments [37], was significantly down regulated in AmTR. Although it is known that this gene is dispensable in several bacterial models, knockout mutants of rnhB show growth rates that are about half of that of the wild type [38]. Another gene with predicted function that was found to be significantly down regulated in AmTR was fpba, which encodes an iron binding protein. The most common strategy of iron acquisition in bacteria is the production of high-affinity ferric iron binding molecules that can sequester iron from the host’s iron-binding proteins. Knock-out of iron acquisition mechanisms has been shown to attenuate pathogenic bacteria and impair their growth [39, 40]. One example of this can be found in Streptococcus pyogenes, where inactivation of an iron transporter affected growth in vivo [41]. The transcriptional status of fpba and rnhB, together with evidence for reduced growth in other bacterial pathogens with impaired iron transport systems or RNase H enzymes, is consistent with the phenotype exhibited by AmTR.

Differentially transcribed groups of genes and pathways

Most phenotypes are the result of more complex genetic changes than divergence at a single gene or locus, or even a few loci [42]. Evaluation of the transcriptional status of groups of genes revealed that genes involved in translation exhibited a trend: of 61 genes, 50 were found to be significantly down regulated in AmTR (Additional file 3). A first step toward addressing the challenge of statistically evaluating genes as groups is the use of clustering techniques. We used K-means clustering in order to validate a trend we were able to identify through observation and find further relationships. Most genes (882) fell into Cluster 1, which is comprised of genes that were either slightly up or down regulated between wild type and AmTR (Figure 2 & Additional file 1). The second cluster included 34 genes that were clearly down regulated in AmTR. Twenty of these genes are involved in translation and translation elongation. Cluster 3 included three genes that were down regulated in AmTR. The difference in average expression values for these genes was much greater than the differences seen in Cluster 2. This cluster included an additional two genes involved in translation. Cluster 4 incorporated two genes that were significantly up regulated in AmTR: AM291 and guaB. Finally, Cluster 5 included 31 genes that were significantly up regulated in AmTR. This cluster included genes involved in translation as well as nucleotide biosynthetic pathways: carB, purL and pyrG.

Clusters 2 and 3 were of particular interest as they identified genes that were down regulated in AmTR compared to wild type. A large proportion (45%) of the genes involved in translation and translation elongation were included in these two clusters, validating our observation. Interestingly, it has been previously shown that mutations and transcriptional alterations in genes involved in translation provoke slow growth of E. coli [43], and three of these, rpoB, rpsL and rpsE, were also found to be altered in AmTR.

Additionally GSEA allowed for identification of other significantly altered pathways such as purine biosynthesis. This is consistent with the finding that two genes involved in the purine biosynthesis pathway are among the top differentially transcribed genes: guaB and purL. Several previous reports have described the importance of nucleotide biosynthesis in bacteria. For instance, certain auxotroph mutants of Salmonella [47], Staphylococcus aureus [48] and Steptococcus pneumoniae [49] have shown impaired growth. These studies suggest that purines and selected amino acids are scarce in different environments and essential to pathogen growth. Impairment in nucleotide biosynthesis would further underlie the slow growth phenotype.