Comparative genomics and transcriptomics of trait-gene association
© Pierlé et al.; licensee BioMed Central Ltd. 2012
Received: 7 August 2012
Accepted: 16 November 2012
Published: 26 November 2012
The Order Rickettsiales includes important tick-borne pathogens, from Rickettsia rickettsii, which causes Rocky Mountain spotted fever, to Anaplasma marginale, the most prevalent vector-borne pathogen of cattle. Although most pathogens in this Order are transmitted by arthropod vectors, little is known about the microbial determinants of transmission. A. marginale provides unique tools for studying the determinants of transmission, with multiple strain sequences available that display distinct and reproducible transmission phenotypes. The closed core A. marginale genome suggests that any phenotypic differences are due to single nucleotide polymorphisms (SNPs). We combined DNA/RNA comparative genomic approaches using strains with different tick transmission phenotypes and identified genes that segregate with transmissibility.
Comparison of seven strains with different transmission phenotypes generated a list of SNPs affecting 18 genes and nine promoters. Transcriptional analysis found two candidate genes downstream from promoter SNPs that were differentially transcribed. To corroborate the comparative genomics approach we used three RNA-seq platforms to analyze the transcriptomes from two A. marginale strains with different transmission phenotypes. RNA-seq analysis confirmed the comparative genomics data and found 10 additional genes whose transcription between strains with distinct transmission efficiencies was significantly different. Six regions of the genome that contained no annotation were found to be transcriptionally active, and two of these newly identified transcripts were differentially transcribed.
This approach identified 30 genes and two novel transcripts potentially involved in tick transmission. We describe the transcriptome of an obligate intracellular bacterium in depth, while employing massive parallel sequencing to dissect an important trait in bacterial pathogenesis.
KeywordsBacteria Rickettsia SNP RNA-seq Anaplasma
The ongoing revolution in genome sequencing has enabled ever-increasing sequence generation at an ever-decreasing cost. The growing availability of fully sequenced genomes offers new opportunities to identify relationships between genotype and phenotype, one of the major goals of the genomics era. Comparative genomics were first introduced as a tool to predict trait-gene associations in 1998 while trying to define species-specific features of Helicobacter pylori. This approach has been used to predict genomic determinants for well-known phenotypes, including hyperthermophily, flagellar motility and pili assembly[2–4]. These studies share the principle that species with similar phenotypes are likely to utilize orthologous genes in the involved biological process. Thus, the simultaneous presence of genes across species would suggest functional similarity among encoded proteins[5, 6]. While these studies illustrate the advantages and applicability of this principle, they are dependent on previous knowledge of the genetic determinants of a specific trait.
The challenge of associating genes with phenotypes has been highlighted by the development of the pangenome concept and the abundance of intraspecies diversity that has been revealed. The pangenome of a bacterial species encompasses the sum of the genetic repertoire found in all strains. Thus, it consists of the core genome found in all the strains plus the “accessory” genes unique to the different strains. Those bacterial species with a high number of accessory genes are termed “open” pangenomes, whereas those lacking strain specific genes are identified as “closed” pangenomes. While the “openness” of the pangenome is an obvious marker of diversity, sequence heterogeneity within the core gene set has also been shown to be relevant to natural genetic variation[8–10]. When several strains of Streptococcus agalactiae were compared to the 2603 VR strain, 99.2% of the total detected single nucleotide polymorphisms (SNPs) were unique to one strain, while none were common to all strains. A similar scenario was found between three strains of Bacillus anthracis, where all SNPs were unique to one strain. As these two organisms, are classical examples of open and closed pangenomes, respectively, this suggests that the SNP profile of a bacterial species can be open regardless of how “locked” their cores are.
An example of an organism with a closed core genome and a high degree of interstrain diversity is Anaplasma marginale, an obligate intracellular pathogen of both domestic and wild ruminants, with a small genome of 1.2 Mb for which the sequence of multiple strains has been determined[8, 11, 12]. No strain-specific genes and no plasmids were found among sensu stricto strains after sequencing of five strains[8, 11]. In contrast, a high degree of allelic diversity was detected: global comparison of five strains revealed a total of 20,082 sites with SNPs detected in at least one of the analyzed strains and, with approximately 6,000 sites between any given pair. The high degree of gene content conservation suggests that phenotypic differences observed in A. marginale must be due to small polymorphisms between strains rather than whole gene insertions or deletions. Therefore, we exploited the interstrain diversity of A. marginale to map the genetic basis underlying phenotypic differences among strains.
A. marginale genome sequences are available for strains that clearly differ in a measurable phenotype: transmission by the arthropod vector. The Saint Maries, Puerto Rico, Virginia, EMø, 6DE and South Idaho strains are examples of efficiently transmitted strains[13–18]. The Florida strain, has been shown to have a very low transmission efficiency as it was not transmitted using >10 times the number of Dermacentor andersoni ticks routinely used for transmission with the St. Maries strain[17, 19, 20].
Due to the complete gene content conservation, differences in transmission efficiency in A. marginale are likely to be ascribed to sequence variation producing variant proteins or affecting gene transcription. Indeed, precedence is seen in bacterial pathogens, where SNPs have been discovered that provide a selective advantage in host colonization. We combined two genomic sequencing approaches in order to find SNPs and transcriptional changes that segregate with transmission phenotype. We first compared the genome sequences of two strains, St. Maries and Florida, which display contrasting phenotypes with respect to the trait of interest, tick transmissibility. Candidate SNPs included polymorphisms encoding non-synonymous substitutions within genes, as well as SNPs located within putative promoter regions. Each SNP on the resulting list was evaluated through comparative genomics in three efficiently transmissible strains for its consistent segregation with phenotype. The remaining differences were sequenced in two additional efficiently transmissible strains. Only SNPs that were unique to the poorly transmissible Florida strain when compared to six efficiently transmitted strains were retained as candidates. This resulted in a list of candidate genes, consisting of those containing candidate SNPs or located downstream of putative promoter SNPs. Transcriptional analysis of candidate genes by RT-PCR revealed genes that were differentially transcribed in strains with distinctly different transmission efficiencies. To find additional transcriptional changes related to the phenotype of interest, we performed a genome wide transcriptome comparison using RNA-seq technology. Total mRNA populations from two A. marginale strains with different transmission capabilities were sequenced using three different platforms. This study makes use of two sequencing approaches and four different technologies to identify genes involved in a relevant microbial trait. We present, to our knowledge, the deepest analysis of an obligate intracellular bacterial transcriptome during the pathogen’s natural course of infection.
Comparative genomics identifies SNPs that segregate with transmission status
Transcription analysis of candidate genes
Reads mapped to A. marginale from three sequencing platforms
A. marginale reads
% of reads mapped to A. marginale
Average matched read length
Percentage of putative 5 ′ UTRs according to length
5′ UTRs < 40 bp
5′ UTRs ≥ 40 bp
5′ UTR within predicted CDS1
Previously unannotated areas that exhibited high transcriptional activity
Identity through blastX
hypothetical protein AmarV_01231 [Anaplasma marginale str. Virginia]
DNA-binding protein HU [Anaplasma phagocytophilum HZ]
hypothetical protein AmarM_02282 [Anaplasma marginale str. Mississippi]
hypothetical protein PseS9_19739 [Pseudomonas sp. S9]
hypothetical protein AmarM_05569 [Anaplasma marginale str. Mississippi]
Transcriptome comparison identifies transcriptional differences between strains with contrasting transmission phenotypes
Candidate genes involved in transmission phenotype segregated by polymorphisms and differential transcription
Differentially transcribed genes
AP, ER, ECh
AP, ER, ECh, ECa
AP, ER, ECh, ECa
AP, ER, ECh, ECa
AP, ER, ECh
AP, ER, ECh, ECa
Outer membrane protein 4
Differentially transcribed genes w/SNPs genescarrying candidate SNPs
Ribosome-associated inhibitor A
ER, ECh, ECa
NADH Dehydrogenase I chain J
AP, ER, ECh, ECa
AP, ER, ECh, ECa
AP, ER, ECh
Genes carrying candidate SNPs
AP, ER, ECh, ECa
ER, ECh, ECa
DNA gyrase B
AP, ER, ECh, ECa, RB, RC, RR
AP, ER, ECh, ECa
AP, ER, ECh, ECa
AP, ER, ECh, ECa, RB, RC, RR
AP, ER, ECh, ECa, RB, RC, RR
AP, ER, ECh, ECa
AP, ER, ECh, ECa, RB, RC, RR
lipoprotein-releasing transmembrane protein
ER, ECh, ECa
RT-PCR and validation of RNA-seq results
The 18 genes that were the most differentially transcribed across replicates were analyzed by using RT-PCR to confirm the RNA-seq results. Fold change in transcription was evaluated and compared with RNA-seq analysis. As shown in Figure3B, transcriptional changes were confirmed and statistically significant in all but one of the analyzed genes. Gene AMF_209, found to be more highly transcribed in the St. Maries strain by 124 fold was not consistently up-regulated across both replicates by RT-PCR (up-regulated by 18.2 and 1.3 fold in separate replicates) and, therefore, its fold change was not statistically significant.
Table4 shows the 30 genes that were selected as candidates. Genes that were found to be differentially transcribed through RNA-seq and RT-PCR are shown on top of Table4; genes with candidate SNPs and differential transcription are shown in the middle of Table4. The rest of the genes contain non-synonymous SNPs that segregated through comparative genomics. The length of the candidate genes varies, with AMF_530 being the longest at 10,479 bp and AMF_1037 the shortest at 240 bp. Twelve of the candidate genes are annotated as hypothetical proteins (Table4). Genes AMF_474, AMF_553, AMF_480, AMF_762, AMF_764, AMF_824, AMF_893 and AMF_878 are orthologs of genes with known functions. Genes downstream from promoter SNPs included one translation inhibitor (AMF_474) and one gene involved in energy consuming processes, nuoJ (AMF_553). Genes containing non-synonymous substitutions included orthologs for DNA gyrase (AMF_480), a tRNA synthase (AMF_762), an aspartate kinase (AMF_764), a carboxypeptidase involved in cell envelope biogenesis (AMF_824) and a lipoprotein releasing protein (AMF_893). A role in transmission is not immediately apparent for these genes, in fact, it is not surprising that more than half of the candidates were of unknown function due to the lack of information on the determinants of tick transmission. A search for related genes revealed that 18 of the candidate genes had homologs in the tick-transmissible human pathogen Anaplasma phagocytophilum (Table4). Ten genes, AMF_051, AMF_433, AMF_432, AMF_431, AMF_430, AMF_429, AMF_547, AMF_613, AMF_762, AMF_893, AMF_798, and AMF_793 also had homologs in tick transmitted Ehrlichia species. Only genes AMF_197, AMF_264, AMF_269, AMF_480, AMF_703, AMF_824 and AMF_893 had homologs in three tick transmitted Rickettsia species. Additionally, hypothetical candidates AMF_1037, AMF_879, AMF_401 and AMF_530 had no homologs in the GenBank database. These findings provide two mutually exclusive scenarios: if a gene with homologs in the aforementioned tick-transmitted organisms is responsible for the trait of interest, this suggests a common mechanism within a bacterial order or family. Alternatively, a gene unique to A. marginale would favor a species-specific scenario.
Four of the genes encoded transmembrane domains and signal peptides predicted through multiple algorithms: AMF_798, AMF_793, AMF_824 and AMF_878. The results obtained for AMF_878 are not surprising, as it is annotated as outer membrane protein 4 (OMP4). Twenty-three genes had significant scores for transmembrane domain predictions but did not contain signal peptides. Analysis of non-synonymous SNPs using the SIFT algorithm predicted eleven to be deleterious; these substitutions are reported in Additional file8.
Pairing comparative genomics with high throughput RNA-seq analysis allows for identification of sequence and transcriptional differences on a genome wide scale. In the present study, comparative genomics reduced a list of candidate SNPs from 9,609 to 49 SNPs that segregate with transmission status, including 35 that encode non-synonymous substitutions within 18 genes and 14 residing within nine putative promoters that could affect transcription of 11 genes. Of the putative promoter SNPs, we retained only those that affected the transcription of adjacent genes, leaving just 2 SNPs affecting two genes, reducing the overall list to 37 candidate SNPs affecting 20 genes. Deep sequencing and comparative expression analysis found an additional 10 genes whose transcription between strains with distinct transmission efficiencies is significantly different. Transcriptome analysis also revealed two previously un-annotated regions that were differentially transcribed between the strains of interest. This produced a final list of 30 genes and two newly identified transcriptionally active regions that segregate with tick transmission.
Our combined approach allowed us to map SNPs that segregate among A. marginale strains with divergent transmission efficiencies. Such subtle differences have been shown to have dramatic effects on organism biology. A single non-synonymous SNP in the envelope protein gene E1 of the Chikungunya virus is directly responsible for a change in vector specificity that caused an epidemic in the Reunion Island in 2004. One SNP in the FimH adhesion gene from a commensal strain of E. coli modified this strain’s affinity for monomannose receptors, correlating directly with increased uroepithelium affinity and allowing detrimental bladder colonization. Similarly, a SNP within the promoter of the nitrate reductase gene cluster narGHIJ was shown to be responsible for the different nitrate reductase phenotype shown by the almost identical Mycobacterium bovis and Mycobacterium tuberculosis, bacterial species with identical gene content.
Comparative genomics identified 20 genes with at least one SNP that segregated with transmission phenotype. The lack of information on the microbial determinants of tick transmission is consistent with the observation that the majority of the genes containing non-synonymous SNPs are of unknown function. Candidate genes with orthologs in other bacterial species do not appear to have an obvious involvement in the phenotype of interest. Three genes: AMF_798, AMF_793 and AMF_824, were predicted to have both signal peptides and transmembrane domains. The presence of signal peptides and transmembrane domains implies membrane localization of the proteins, and thus, these proteins would be more likely to interact with vector molecules and therefore effect transmission. Out of the 35 non-synonymous candidate SNPs, a little under a third were predicted to be deleterious (Additional file8). Gene AMF_1026 carries the highest number of deleterious substitutions with a total of three non-synonymous SNPs. This gene was also found to be up-regulated in the efficiently transmitted strain through RT-PCR. Interestingly, it also had a SNP in its promoter region. This promoter SNP did not segregate with the rest of the transmissible strains and therefore was not retained as a candidate. Polymorphisms were retained as candidates if six efficiently transmissible strains consistently diverged with the nucleotide found in the poorly transmitted Florida strain. Candidate SNPs included non-synonymous changes in ORFs and SNPs found in putative promoter regions.
Two genes with candidate SNPs in their putative promoter regions were found to be differentially transcribed. AMF_553, more highly transcribed in the St. Maries strain, is annotated as NADH dehydrogenase I chain J (nuoJ). This is part of the membrane arm of respiratory complex I, a conserved proton pumping NADH:ubiquinone oxidoreductase in bacteria. Another closely associated gene from this complex, nuoL, has been found to be up-regulated in the related organism Rickettsia conorii while dealing with osmotic stress, suggesting that enhancement of NADH dehydrogenase expression in a vector-transmitted bacterium could be related to an adaptation strategy necessary to survive in the changing osmolarity of a feeding tick. AMF_474, more highly transcribed in the Florida strain, contains conserved domains for a modulation protein, the ribosome associated inhibitor A (RaiA) also known as protein Y (PY). This protein is a cold-shock induced ribosome binding protein that inhibits translation. PY binds exclusively to the 30S subunit of the 70S ribosome, and preventing the formation of initiation complexes by preventing the binding of mRNA and initiator fMet-tRNA to the ribosome. When temperature levels return to 37°C, initiation of protein synthesis overcomes the PY inhibition as tRNA compete more effectively with PY in elevated temperatures. Related bacterial species which are also transmitted by D. andersoni, such as Rickettsia rickettsii, are known to enter “dormant” stages within ticks. Subsequent reversion of this state, in a process termed “reactivation”, is thought to be due to an increase in temperature when the arthropod feeds on the mammalian host. Therefore these observations suggest an interesting scenario as this gene was up-regulated in the low transmission efficiency Florida strain. The low transmission phenotype could be due to a halt in translation produced by an up regulation of PY during cold shock response.
The three aforementioned differentially transcribed genes were identified through comparative genomics. Although all three carried SNPs in their promoter regions, only two were retained as candidates. This exposes a limitation of the approach that was used in this study: polymorphisms that do not segregate with all the highly transmissible strains may still contribute to the phenotype of interest. In order to confirm the differences in transcription revealed through RT-PCR and to find further changes in transcriptional activity that the strategy might have overlooked; the transcriptome of two strains with contrasting transmission phenotypes were compared. Genome wide comparison of transcriptional activity confirmed our RT-PCR results and found an additional 10 genes that were significantly differentially transcribed. Of these 10 genes only one had a predicted localization: AMF_878 corresponds to OMP4, an outer membrane protein and member of the pfam 01617 superfamily. Among the remaining genes with no functional annotation three genes stood out as they exhibited a complete lack of transcriptional activity in the poorly transmitted Florida strain: AMF_431, AMF_432 and AMF_433. These genes appear to be arranged in an operon along with AMF_429 and AMF_430, according to the tiling of reads mapped in the St. Maries strain. AMF_429 and AMF_430 were also significantly differentially transcribed between the strains. Genes AMF_429, AMF_431 and AMF_433 contain high scoring conserved domains for tail and head/tail connector phage proteins, with the highest similarity found to phage proteins from Wolbachia spp., a related bacterial symbiont of arthropods. Although this could open interesting possibilities, as phages play an important role for Wolbachia spp. within the arthropod host, no mobile elements or intact prophages have been identified in A. marginale.
Typically, pathogenic bacteria that cycle between arthropod and mammalian hosts modify their transcriptional profiles to adapt to these different environments. One of the major difficulties involved in examining gene regulation of obligate intracellular pathogens is the low amounts of bacterial RNA, which is co-isolated with large amounts of host RNA. In order to overcome the limited amount of bacterial RNA, previous transcriptomic studies interrogating genes used for obligate intracellular survival were conducted using mimetic conditions of infection in an in vitro environment[37–39]. While these studies provide insight into a limited number of genes regulated by specific cues, they are not representative of natural infection. Exposing the related pathogen R. rickettsii to different environmental conditions that mimic its transition from arthropod to mammalian host showed a surprisingly minimal transcriptional response, with less than 10 genes changing more than 3-fold in expression level. This could indicate that pathogens in the order Rickettsiales do not regulate genes specifically for growth within mammalian or tick cells but contain a conserved set of genes that are required for growth in both environments. The obligate intracellular habitat of pathogens in this Order may offer such a stable environment that the necessity for gene regulation is much less than that of facultative intracellular pathogens. Our study searched for transcriptional differences between strains with contrasting transmission profiles in the natural host of our model organism.
The use of different sequencing platforms in this study was instrumental in confirming significant and consistent changes in transcriptional activity. It has been shown that different RNA preparation and selection procedures in deep sequencing experiments can lead to measurable over- or under-representation of particular RNAs. This study proved that utilizing different technologies allowed for control of sources of potential bias in RNA sequencing: all three platforms used for our study gave the same results. Making use of various platforms was also instrumental in our goal of describing the A. marginale transcriptome with the highest possible accuracy. In bacteria, the overwhelmingly high numbers of reads in combination with relatively small genome sizes has led to the assumption that complete or nearly complete transcriptomes are being analyzed. However, selecting for prokaryotic sequences in an ocean of eukaryotic RNAs makes accurate representation of RNA populations daunting. Few attempts have been made at describing the transcriptome of obligate intracellular pathogens through RNA-seq; notably, to date, this has been done for Chlamydia species[38, 39] and the tick-transmitted pathogen A. phagocytophilum. The deepest analysis generated 854,242 reads that mapped to the 1.23 Mb Chlamydia pneumonia genome; we mapped up to 2,990,921 reads per replicate to A. marginale’s 1.2 Mb genome. To enrich for prokaryotic sequences, previous attempts at characterizing obligate intracellular microbial transcriptomes used differential centrifugation of in vitro grown bacteria in order to separate the bacteria from host cells. This procedure is likely to stress the bacteria and skew their transcriptional profile. Enrichment for our samples was performed by selective hybridization once RNA populations were collected. Although Mastronunzio et al. used a similar enrichment procedure; they only detected 187,284 reads, representing 11% of the CDSs in the A. phagocytophilum genome. In this study, 99% of the CDSs in the A. marginale were detected through transcriptional analysis.
Analyzing transcriptional profiles with RNA-seq allows us to evaluate “snapshots” in time of bacterial transcriptomes; therefore, it is essential to generate data from more than one replicate to provide a broader more reliable picture of transcriptional changes. The depth and reproducibility of this RNA-seq data set allowed for mapping of the physical structure of the A. marginale transcriptome; including previously unreported transcriptionally active regions and 5′ UTR length. Six regions with no previous annotation were detected in both strains; two of these were differentially transcribed. The role of these transcripts is uncertain as only two of these were predicted to contain ORFs. The majority of the high confidence 5′ UTRs were longer than 40 bp in both strains. Previous studies of TSSs have shown that only a very small portion of 5′ UTRs are longer than 40 bp in bacteria[41, 42]. As 5′ UTRs have been involved in regulation processes in bacteria, further investigation of these elements might reveal translational and transcriptional roles. Additionally, mapping of transcriptional data allowed us to define 70 putative operon structures that involved 292 genes, showing that at least 30% of the genes are polycistronic. Although RNA-seq allows us to study polycistronic messages on a genome wide scale, the depth of this technique coupled with tiling arrays have shown that the concept of simple operons is questionable. Differential expression of consecutive genes within operons and condition dependent modulation highlight the complexity of transcriptional regulation in bacteria.
This study takes advantage of the high interstrain diversity of this intracellular bacterium to significantly reduce the number of candidate differences that could be involved in the tick transmission phenotype. Marrying next generation sequencing approaches allowed us to generate a list of genes differing at the transcriptional and sequence levels in strains with contrasting transmission status. Transformation of the transmission deficient allele into a transmission competent strain will facilitate functional analysis of these genes in order to determine their role in transmission by the arthropod host. Although the successful transformation of A. marginale has been achieved[45, 46], stable targeted gene replacement has not been accomplished and is a necessary next step for determining the role of these genes in tick transmission. Identification of genes involved in tick transmission in our model will provide an important first step toward the development of novel control strategies for tick-borne pathogens, such as transmission-blocking vaccines.
Animal experiments were approved by the Institutional Animal Care and Use Committee at University of Idaho, USA, in accordance with institutional guidelines based on the U.S. National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals.
The Florida, St. Maries, Virginia, Puerto Rico, South Idaho, EMΦ and 6DE strains used in this study have been described in detail elsewhere[47–51]. The St. Maries, Virginia, Puerto Rico, South Idaho, EMΦ and 6DE strains are reproducibly transmitted by the Reynolds Creek stock of D. andersoni[13, 14, 16, 17, 52, 53]. The Florida strain has not been successfully transmitted by any tick species, including the Reynolds Creek stock[15, 18, 19].
The accession numbers for the strains used are: St. Maries: CP000030.1, Florida: CP001079.1, Viriginia: ABOR00000000.1, Puerto Rico: ABOQ00000000.1, South Idaho: AFMY00000000.1. MUMmer v3.1 was used to compare as previously described to compare the Florida and St. Maries strains. SNPs encoding synonymous substitutions were not further analyzed. The runMapping program of the Newbler suite v2.5.3 (454 Life Sciences) was used with default settings to compare all reads from the Virginia, Puerto Rico, and South Idaho strains to the completed Florida and St. Maries genomes. All remaining SNPs from the initial comparison were then checked against the three strains; if the Florida sequence was matched in any of the highly transmissible strains, that SNP was removed as a candidate. Illumina sequencing of the St. Maries strain was used to evaluate the frequencies of the SNPs found between the Florida and St. Maries strain. SNPs that were found at 100% frequencies were highlighted in Additional file1.
The remaining SNPs were examined via targeted sequencing of the South Idaho, EMΦ and 6DE strains[50, 55]. Primers were designed by aligning the SNP-containing region from the Florida and St. Maries strains and selecting primers to flank the polymorphism. The resulting amplicons were generated from genomic DNA, cloned into pCR4-TOPO (Invitrogen) and sequenced in both directions using BigDye v3.1 chemistry on an ABI 3130XL (Applied Biosystems). Sequence analysis eliminated candidates as described above. All candidate SNPs were resequenced in the Florida strain, to verify the original genomic sequence.
Comparative transcriptional analysis
Total RNA was isolated from A. marginale-infected blood using TRIzol (Invitrogen), per manufacturer directions. Expression was measured using quantitative reverse transcription PCR using the SYBR Green ER RT-PCR Kit (Invitrogen). Briefly, 1 ug of RNA was processed with the Superscript III First strand kit (Invitrogen) to obtain cDNA. Copy numbers were corrected to more closely reflect transcript levels based on reverse transcription efficiency. The steady state, single copy gene msp 5 was used to calibrate the RT-PCR. Relative expression ratios were calculated by a mathematical model, which includes efficiency correction of individual transcripts through the REST software. This software uses the Pair Wise Fixed Reallocation Randomization Test to assess the statistical significance of the RT-PCR results when comparing the relative expression of the promoter candidates in both the Florida and St. Maries strains. A differential expression fold cutoff value of 3.2 was established by calculating the mean of the average ratios observed for all genes analyzed in this study plus 2 standard deviations. In order to assess the statistical relevance of the findings across two biological replicates, an adaptation of the standardization method proposed by Willems and coworkers was used; this includes three basic steps: log transformation, mean centering and autoscalling. After standardizing the data, statistical significance of the fold changes observed between the strains across both experiments was determined by calculation of 95% confidence intervals. This procedure was applied to each candidate gene and was also used for verification of transcriptional differences found by RNA-seq.
The accession number for this RNA-seq study is: SRP014580. Two Holstein calves negative for A. marginale by MSP5 cELISA, C1322 and C1323, were inoculated with the Florida and the St. Maries strains, respectively. Infection levels were tracked by analysis of Giemsa-stained blood smears to calculate the percentage of parasitized erythrocytes (PPE). Blood samples were taken at similar levels of parasitemia (3.5 and 4% PPE). Total RNA was isolated from A. marginale-infected blood using TRIzol (Invitrogen) per the manufacturer’s directions. Eukaryotic sequences were negatively selected through hybridization using the MICROBEnrich kit (Ambion). For samples processed for 454 and Ion Torrent technologies, probes for bacterial ribosomal RNAs from the Ribominus kit (Invitrogen) were added during the subtractive hybridization procedure. For samples processed for Illumina, the Duplex‐Specific thermostable nuclease (DSN) normalization protocol was applied. 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 (Reads Per Kilobase per Million mapped reads) values were calculated. Two different tests were applied to evaluate the statistical significance of fold changes: Kal’s and Baggerly’s statistical tests on proportions[58, 59]. Comparisons of replicates were performed in order to account for variation within a strain. These comparisons showed very little variation: a maximum of 2% of genes had fold changes above or below 1. As variation within strains was assessed we proceeded to compare the differentially transmitted strains. In order to establish transcription fold change cutoffs, the relationship between the p-values of the statistical tests applied and the magnitude of the difference in expression values of the samples was plotted and evaluated. This was done in order to arrange genes along dimensions of biological and statistical significance. Genes whose log2 fold change was above and below 2 and -2, respectively, and whose -log10 p-value was above 10 in both replicate comparisons and under both statistical tests were selected for further evaluation (Additional file7).
Areas of the genome that were not previously annotated and showed >0.5 coverage (average sequence data coverage depth) were reported when reads were unambiguously mapped to the A. marginale genome.
The relative performances noted in Table1 for the different sequencing technologies should not be directly compared, as this study was not designed to compare these platforms. As has been noted, different library preparations and sequencing technologies favor recovery of different transcripts. The goal of using multiple technologies was to verify that under- or over-represented transcripts in any strain were not being favored by the technology used.
Putative start site identification
Putative transcript start sites were identified using the rules proposed by Passalacqua et al.: briefly; genes with continuous coverage extending into a codirectional upstream gene were identified as members of an operon. If the signal “dropped off” in the intergenic sequence upstream of the open reading frame, we designated the point at which coverage dropped to 0 as the putative transcriptional start site. Coverage depth was calculated for every position of each genome, and all genes considered had an average coverage score >0.5 above the calculated average coverage signal. Putative TSSs that were found with the highest confidence (i.e. TSSs present in all replicates) were grouped in two different tables according to the length of the 5′ UTRs, less or more than 40 bp.
Bioinformatic analysis of candidates
In order to rank the candidates, two different criteria were established. The first, termed “biological plausibility of association”, examines the annotation of the currently available genomes and the predicted function of the candidate gene, using existing knowledge about biology and the studied phenotype. In other words, is the candidate gene likely to be involved in the examined phenotype according to its known or predicted function? The second criterion involves the use of three in silico analyses. The presence of signal peptides in the candidate genes was assessed by using SignalP 4.0. Transmembrane domains were predicted using two distinct algorithms: TMpred and Dense Alignment Surface (DAS) methods; only genes with transmembrane domains predicted by both algorithms were reported. The “Sorting Tolerant From Intolerant” (SIFT) algorithm uses a sequence homology-based approach to classify amino acid substitutions, and was used to predict if substitutions in the candidate alleles detrimental or tolerated by the protein. The search for ORFs in newly identified transcriptionally active regions was performed using three different tools: CLC Genomics Workbench (CLC Bio), NIH’s ORF finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) and ORF (http://bioinformatics.biol.rug.nl/websoftware/orf/orf_start.php).
Coding DNA Sequence
Open Reading Frame
Reads Per Kilobase per Million mapped reads
The authors would like to acknowledge the expert technical assistance of Ms. Xiaoya Cheng. This work was supported by USDA CREES NRI CGP 2004-35600-14175 and 2005-35604-15440, National Institutes of Health Grant AI44005, and Wellcome Trust GR075800M. SAP was supported in part by fellowships from the Poncin Trust and CONACyT.
- Bork P, Dandekar T, Diaz-Lazcoz Y, Eisenhaber F, Huynen M, Yuan Y: Predicting function: from genes to genomes and back. J Mol Biol. 1998, 283 (4): 707-725. 10.1006/jmbi.1998.2144.View ArticlePubMed
- Jim K, Parmar K, Singh M, Tavazoie S: A cross-genomic approach for systematic mapping of phenotypic traits to genes. Genome Res. 2004, 14 (1): 109-115.PubMed CentralView ArticlePubMed
- Levesque M, Shasha D, Kim W, Surette MG, Benfey PN: Trait-to-gene: a computational method for predicting the function of uncharacterized genes. Curr Biol. 2003, 13 (2): 129-133. 10.1016/S0960-9822(03)00009-5.View ArticlePubMed
- Makarova KS, Wolf YI, Koonin EV: Potential genomic determinants of hyperthermophily. Trends Genet. 2003, 19 (4): 172-176. 10.1016/S0168-9525(03)00047-7.View ArticlePubMed
- Huynen MA, Bork P: Measuring genome evolution. Proc Natl Acad Sci U S A. 1998, 95 (11): 5849-5856. 10.1073/pnas.95.11.5849.PubMed CentralView ArticlePubMed
- Pellegrini M, Marcotte EM, Thompson MJ, Eisenberg D, Yeates TO: Assigning protein functions by comparative genome analysis: protein phylogenetic profiles. Proc Natl Acad Sci U S A. 1999, 96 (8): 4285-4288. 10.1073/pnas.96.8.4285.PubMed CentralView ArticlePubMed
- Tettelin H, Masignani V, Cieslewicz MJ, Donati C, Medini D, Ward NL, Angiuoli SV, Crabtree J, Jones AL, Durkin AS, et al: Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial “pan-genome”. Proc Natl Acad Sci U S A. 2005, 102 (39): 13950-13955. 10.1073/pnas.0506758102.PubMed CentralView ArticlePubMed
- Dark MJ, Herndon DR, Kappmeyer LS, Gonzales MP, Nordeen E, Palmer GH, Knowles DP, Brayton KA: Conservation in the face of diversity: multistrain analysis of an intracellular bacterium. BMC Genomics. 2009, 10: 16-10.1186/1471-2164-10-16.PubMed CentralView ArticlePubMed
- Pearson T, Busch JD, Ravel J, Read TD, Rhoton SD, U’Ren JM, Simonson TS, Kachur SM, Leadem RR, Cardon ML, et al: Phylogenetic discovery bias in Bacillus anthracis using single-nucleotide polymorphisms from whole-genome sequencing. Proc Natl Acad Sci U S A. 2004, 101 (37): 13536-13541. 10.1073/pnas.0403844101.PubMed CentralView ArticlePubMed
- Stratilo CW, Lewis CT, Bryden L, Mulvey MR, Bader D: Single-nucleotide repeat analysis for subtyping Bacillus anthracis isolates. J Clin Microbiol. 2006, 44 (3): 777-782. 10.1128/JCM.44.3.777-782.2006.PubMed CentralView ArticlePubMed
- Brayton KA, Kappmeyer LS, Herndon DR, Dark MJ, Tibbals DL, Palmer GH, McGuire TC, Knowles DP: Complete genome sequencing of Anaplasma marginale reveals that the surface is skewed to two superfamilies of outer membrane proteins. Proc Natl Acad Sci U S A. 2005, 102 (3): 844-849. 10.1073/pnas.0406656102.PubMed CentralView ArticlePubMed
- Dark MJ, Al-Khedery B, Barbet AF: Multistrain genome analysis identifies candidate vaccine antigens of Anaplasma marginale. Vaccine. 2011, 29 (31): 4923-4932. 10.1016/j.vaccine.2011.04.131.PubMed CentralView ArticlePubMed
- Eriks IS, Stiller D, Palmer GH: Impact of persistent Anaplasma marginale rickettsemia on tick infection and transmission. J Clin Microbiol. 1993, 31 (8): 2091-2096.PubMed CentralPubMed
- Futse JE, Brayton KA, Dark MJ, Knowles DP, Palmer GH: Superinfection as a driver of genomic diversification in antigenically variant pathogens. Proc Natl Acad Sci U S A. 2008, 105 (6): 2123-2127. 10.1073/pnas.0710333105.PubMed CentralView ArticlePubMed
- Futse JE, Ueti MW, Knowles DP, Palmer GH: Transmission of Anaplasma marginale by Boophilus microplus: retention of vector competence in the absence of vector-pathogen interaction. J Clin Microbiol. 2003, 41 (8): 3829-3834. 10.1128/JCM.41.8.3829-3834.2003.PubMed CentralView ArticlePubMed
- Rurangirwa FR, Stiller D, Palmer GH: Strain diversity in major surface protein 2 expression during tick transmission of Anaplasma marginale. Infect Immun. 2000, 68 (5): 3023-3027. 10.1128/IAI.68.5.3023-3027.2000.PubMed CentralView ArticlePubMed
- Scoles GA, Broce AB, Lysyk TJ, Palmer GH: Relative efficiency of biological transmission of Anaplasma marginale (Rickettsiales: Anaplasmataceae) by Dermacentor andersoni (Acari: Ixodidae) compared with mechanical transmission by Stomoxys calcitrans (Diptera: Muscidae). J Med Entomol. 2005, 42 (4): 668-675. 10.1603/0022-2585(2005)042[0668:REOBTO]2.0.CO;2.View ArticlePubMed
- Scoles GA, Ueti MW, Noh SM, Knowles DP, Palmer GH: Conservation of transmission phenotype of Anaplasma marginale (Rickettsiales: Anaplasmataceae) strains among Dermacentor and Rhipicephalus ticks (Acari: Ixodidae). J Med Entomol. 2007, 44 (3): 484-491. 10.1603/0022-2585(2007)44[484:COTPOA]2.0.CO;2.View ArticlePubMed
- Ueti MW, Reagan JO, Knowles DP, Scoles GA, Shkap V, Palmer GH: Identification of midgut and salivary glands as specific and distinct barriers to efficient tick-borne transmission of Anaplasma marginale. Infect Immun. 2007, 75 (6): 2959-2964. 10.1128/IAI.00284-07.PubMed CentralView ArticlePubMed
- Wickwire KB, Kocan KM, Barron SJ, Ewing SA, Smith RD, Hair JA: Infectivity of three Anaplasma marginale isolates for Dermacentor andersoni. Am J Vet Res. 1987, 48 (1): 96-99.PubMed
- Sokurenko EV, Hasty DL, Dykhuizen DE: Pathoadaptive mutations: gene loss and variation in bacterial pathogens. Trends Microbiol. 1999, 7 (5): 191-195. 10.1016/S0966-842X(99)01493-6.View ArticlePubMed
- Ramabu SS, Ueti MW, Brayton KA, Baszler TV, Palmer GH: Identification of Anaplasma marginale proteins specifically up-regulated during colonization of the tick vector. Infect Immun. 2010, 78 (7): 3047-3052. 10.1128/IAI.00300-10.PubMed CentralView ArticlePubMed
- Visser ES, McGuire TC, Palmer GH, Davis WC, Shkap V, Pipano E, Knowles DP: The Anaplasma marginale msp5 gene encodes a 19-kilodalton protein conserved in all recognized Anaplasma species. Infect Immun. 1992, 60 (12): 5139-5144.PubMed CentralPubMed
- Willems E, Leyns L, Vandesompele J: Standardization of real-time PCR gene expression data from independent biological replicates. Anal Biochem. 2008, 379 (1): 127-129. 10.1016/j.ab.2008.04.036.View ArticlePubMed
- Raabe CA, Hoe CH, Randau G, Brosius J, Tang TH, Rozhdestvensky TS: The rocks and shallows of deep RNA sequencing: Examples in the Vibrio cholerae RNome. RNA (New York, NY. 2011, 17 (7): 1357-1366. 10.1261/rna.2682311.View Article
- Kumar P, Henikoff S, Ng PC: Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat Protoc. 2009, 4 (7): 1073-1081.View ArticlePubMed
- Tsetsarkin KA, Vanlandingham DL, McGee CE, Higgs S: A single mutation in chikungunya virus affects vector specificity and epidemic potential. PLoS pathogens. 2007, 3 (12): e201-10.1371/journal.ppat.0030201.PubMed CentralView ArticlePubMed
- Stermann M, Sedlacek L, Maass S, Bange FC: A promoter mutation causes differential nitrate reductase activity of Mycobacterium tuberculosis and Mycobacterium bovis. J Bacteriol. 2004, 186 (9): 2856-2861. 10.1128/JB.186.9.2856-2861.2004.PubMed CentralView ArticlePubMed
- Kao MC, Di Bernardo S, Nakamaru-Ogiso E, Miyoshi H, Matsuno-Yagi A, Yagi T: Characterization of the membrane domain subunit NuoJ (ND6) of the NADH-quinone oxidoreductase from Escherichia coli by chromosomal DNA manipulation. Biochemistry. 2005, 44 (9): 3562-3571. 10.1021/bi0476477.View ArticlePubMed
- Renesto P, Rovery C, Schrenzel J, Leroy Q, Huyghe A, Li W, Lepidi H, Francois P, Raoult D: Rickettsia conorii transcriptional response within inoculation eschar. PLoS One. 2008, 3 (11): e3681-10.1371/journal.pone.0003681.PubMed CentralView ArticlePubMed
- Munderloh UG, Kurtti TJ: Cellular and molecular interrelationships between ticks and prokaryotic tick-borne pathogens. Annu Rev Entomol. 1995, 40: 221-243. 10.1146/annurev.en.40.010195.001253.View ArticlePubMed
- Agafonov DE, Kolb VA, Spirin AS: Ribosome-associated protein that inhibits translation at the aminoacyl-tRNA binding stage. EMBO Reports. 2001, 2 (5): 399-402.PubMed CentralView ArticlePubMed
- Vila-Sanjurjo A, Schuwirth BS, Hau CW, Cate JH: Structural basis for the control of translation initiation during stress. Nat Struct Mol Biol. 2004, 11 (11): 1054-1059. 10.1038/nsmb850.View ArticlePubMed
- McDade JE, Newhouse VF: Natural history of Rickettsia rickettsii. Annu Rev Microbiol. 1986, 40: 287-309. 10.1146/annurev.mi.40.100186.001443.View ArticlePubMed
- Bordenstein SR, Marshall ML, Fry AJ, Kim U, Wernegreen JJ: The tripartite associations between bacteriophage, Wolbachia, and arthropods. PLoS pathogens. 2006, 2 (5): e43-10.1371/journal.ppat.0020043.PubMed CentralView ArticlePubMed
- Obonyo M, Munderloh UG, Fingerle V, Wilske B, Kurtti TJ: Borrelia burgdorferi in tick cell culture modulates expression of outer surface proteins A and C in response to temperature. J Clin Microbiol. 1999, 37 (7): 2137-2141.PubMed CentralPubMed
- Ellison DW, Clark TR, Sturdevant DE, Virtaneva K, Hackstadt T: Limited transcriptional responses of Rickettsia rickettsii exposed to environmental stimuli. PLoS One. 2009, 4 (5): e5612-10.1371/journal.pone.0005612.PubMed CentralView ArticlePubMed
- Albrecht M, Sharma CM, Reinhardt R, Vogel J, Rudel T: Deep sequencing-based discovery of the Chlamydia trachomatis transcriptome. Nucleic Acids Res. 2010, 38 (3): 868-877. 10.1093/nar/gkp1032.PubMed CentralView ArticlePubMed
- Albrecht M, Sharma CM, Dittrich MT, Muller T, Reinhardt R, Vogel J, Rudel T: The transcriptional landscape of Chlamydia pneumoniae. Genome Biol. 2011, 12 (10): R98-10.1186/gb-2011-12-10-r98.PubMed CentralView ArticlePubMed
- Mastronunzio JE, Kurscheid S, Fikrig E: Post-genomic analyses reveal development of infectious Anaplasma phagocytophilum during transmission from ticks to mice. J Bacteriol. 2012, 194 (9): 2238-2247. 10.1128/JB.06791-11.PubMed CentralView ArticlePubMed
- McGrath PT, Lee H, Zhang L, Iniesta AA, Hottes AK, Tan MH, Hillson NJ, Hu P, Shapiro L, McAdams HH: High-throughput identification of transcription start sites, conserved promoter motifs and predicted regulons. Nat Biotechnol. 2007, 25 (5): 584-592. 10.1038/nbt1294.View ArticlePubMed
- Passalacqua KD, Varadarajan A, Ondov BD, Okou DT, Zwick ME, Bergman NH: Structure and complexity of a bacterial transcriptome. J Bacteriol. 2009, 191 (10): 3203-3211. 10.1128/JB.00122-09.PubMed CentralView ArticlePubMed
- Naville M, Gautheret D: Transcription attenuation in bacteria: theme and variations. Brief Funct Genomic Proteomic. 2009, 8 (6): 482-492. 10.1093/bfgp/elp025.View ArticlePubMed
- Guell M, Yus E, Lluch-Senar M, Serrano L: Bacterial transcriptomics: what is beyond the RNA horiz-ome?. Nature Reviews. 2011, 9 (9): 658-669. 10.1038/nrmicro2620.PubMed
- Felsheim RF, Chavez AS, Palmer GH, Crosby L, Barbet AF, Kurtti TJ, Munderloh UG: Transformation of Anaplasma marginale. Vet Parasitol. 2010, 167 (2–4): 167-174.PubMed CentralView ArticlePubMed
- Noh SM, Ueti MW, Palmer GH, Munderloh UG, Felsheim RF, Brayton KA: Stability and tick transmission phenotype of gfp-transformed Anaplasma marginale through a complete in vivo infection cycle. Appl Environ Microbiol. 2011, 77 (1): 330-334. 10.1128/AEM.02096-10.PubMed CentralView ArticlePubMed
- Eriks IS, Stiller D, Goff WL, Panton M, Parish SM, McElwain TF, Palmer GH: Molecular and biological characterization of a newly isolated Anaplasma marginale strain. J Vet Diagn Invest. 1994, 6 (4): 435-441. 10.1177/104063879400600406.View ArticlePubMed
- Kuttler KL, Winward LD: Serologic comparisons of 4 Anaplasma isolates as measured by the complement-fixation test. Vet Microbiol. 1984, 9 (2): 181-186. 10.1016/0378-1135(84)90033-6.View ArticlePubMed
- McGuire TC, Palmer GH, Goff WL, Johnson MI, Davis WC: Common and isolate-restricted antigens of Anaplasma marginale detected with monoclonal antibodies. Infect Immun. 1984, 45 (3): 697-700.PubMed CentralPubMed
- Palmer GH, Knowles DP, Rodriguez JL, Gnad DP, Hollis LC, Marston T, Brayton KA: Stochastic transmission of multiple genotypically distinct Anaplasma marginale strains in a herd with high prevalence of Anaplasma infection. J Clin Microbiol. 2004, 42 (11): 5381-5384. 10.1128/JCM.42.11.5381-5384.2004.PubMed CentralView ArticlePubMed
- Ristic CCA: Methods of immunoprophylaxis against bovine anaplasmosis with emphasis on use of attenuated Anaplasma marginale vaccine. Adv Exp Med Biol. 1977, 93: 151-188. 10.1007/978-1-4615-8855-9_10.View ArticlePubMed
- Lohr CV, Rurangirwa FR, McElwain TF, Stiller D, Palmer GH: Specific expression of Anaplasma marginale major surface protein 2 salivary gland variants occurs in the midgut and is an early event during tick transmission. Infect Immun. 2002, 70 (1): 114-120. 10.1128/IAI.70.1.114-120.2002.PubMed CentralView ArticlePubMed
- Ueti MW, Knowles DP, Davitt CM, Scoles GA, Baszler TV, Palmer GH: Quantitative differences in salivary pathogen load during tick transmission underlie strain-specific variation in transmission efficiency of Anaplasma marginale. Infect Immun. 2009, 77 (1): 70-75. 10.1128/IAI.01164-08.PubMed CentralView ArticlePubMed
- Kurtz S, Phillippy A, Delcher AL, Smoot M, Shumway M, Antonescu C, Salzberg SL: Versatile and open software for comparing large genomes. Genome Biol. 2004, 5 (2): R12-10.1186/gb-2004-5-2-r12.PubMed CentralView ArticlePubMed
- Rodriguez JL, Palmer GH, Knowles DP, Brayton KA: Distinctly different msp2 pseudogene repertoires in Anaplasma marginale strains that are capable of superinfection. Gene. 2005, 361: 127-132.View ArticlePubMed
- Pfaffl MW, Horgan GW, Dempfle L: Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002, 30 (9): e36-10.1093/nar/30.9.e36.PubMed CentralView ArticlePubMed
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods. 2008, 5 (7): 621-628. 10.1038/nmeth.1226.View ArticlePubMed
- Baggerly KA, Deng L, Morris JS, Aldaz CM: Differential expression in SAGE: accounting for normal between-library variation. Bioinformatics (Oxford, England). 2003, 19 (12): 1477-1483. 10.1093/bioinformatics/btg173.View Article
- Kal AJ, van Zonneveld AJ, Benes V, van den Berg M, Koerkamp MG, Albermann K, Strack N, Ruijter JM, Richter A, Dujon B, et al: Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources. Mol Biol Cell. 1999, 10 (6): 1859-1872.PubMed CentralView ArticlePubMed
- Cui X, Churchill GA: Statistical tests for differential expression in cDNA microarray experiments. Genome Biol. 2003, 4 (4): 210-10.1186/gb-2003-4-4-210.PubMed CentralView ArticlePubMed
- Tabor HK, Risch NJ, Myers RM: Candidate-gene approaches for studying complex genetic traits: practical considerations. Nat Rev Genet. 2002, 3 (5): 391-397. 10.1038/nrg796.View ArticlePubMed
- Emanuelsson O, Brunak S, von Heijne G, Nielsen H: Locating proteins in the cell using TargetP, SignalP and related tools. Nat Protoc. 2007, 2 (4): 953-971. 10.1038/nprot.2007.131.View ArticlePubMed
- Cserzo M, Wallin E, Simon I, von Heijne G, Elofsson A: Prediction of transmembrane alpha-helices in prokaryotic membrane proteins: the dense alignment surface method. Protein Eng. 1997, 10 (6): 673-676. 10.1093/protein/10.6.673.View ArticlePubMed
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.