Comparative genomics and transcriptomics of trait-gene association

Background 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. Results 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. Conclusions 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.


Background
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 [1]. This approach has been used to predict genomic determinants for well-known phenotypes, including hyperthermophily, flagellar motility and pili assembly [2][3][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 [7]. 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][9][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][14][15][16][17][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 [21]. 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
Comparison of the poorly transmissible Florida strain with the efficiently transmitted St. Maries strain produced a total of 9,609 SNPs evenly distributed throughout the genome ( Figure 1, Figure 2, and Additional file 1). Two types of SNPs were further characterized: those that resulted in non-synonymous amino acid changes within genes and SNPs located in putative promoter regions. For the purposes of this study, putative promoters were defined as intergenic regions immediately 5 0 to translation start sites. Global comparison of these SNPs with genome sequences of three efficiently transmitted strains, Puerto Rico, Virginia and South Idaho yielded 241 NS changes within genes, and 62 SNPs distributed in 27 putative promoters. These genes and promoters were then further analyzed in two additional efficiently transmitted strains, 6DE and EMø, by performing targeted sequencing of the regions of interest. The final candidate list included 18 genes that contained at least one SNP encoding a non-synonymous substitution that segregated with transmission status, and 14 SNPs within nine intergenic regions that could potentially affect the transcription of 11 genes (Figure 1, Additional file 1).

Transcription analysis of candidate genes
These 29 genes with SNPs in their coding regions or in their putative promoter regions were analyzed for transcriptional activity by using RT-PCR, which revealed that the 29 candidate genes were transcribed in both the efficiently transmitted St. Maries and the poorly transmissible Florida strain (Additional file 2). For the 11 genes flanking candidate promoter regions, the relative expression ratio was analyzed from two separate infections using msp5 as a steady state calibrator [22,23]. The fold changes were tested for statistical significance by the pairwise randomization test in two separate infections. Statistical significance of the average fold changes across both biological replicates was tested using an adaptation of the method proposed by Willems et al. [24] ( Figure 3A). Four genes were differentially expressed in two biological replicates: AMF_553 showed 4.3 times increased expression in the efficiently transmitted strain (P < 0.05). AMF_474, AMF_505 and AMF_142 showed decreased expression in the highly transmissible strain by ratios of 0.2, 0.6 and 0.7 respectively (P < 0.05). We calculated an expression cutoff by adding 2 standard deviations to the average fold change seen in all the studied genes. Of these differentially expressed candidates, only genes AMF_474 and AMF_553 were below and above the calculated cutoff, respectively.

RNA-seq
The transcriptomes of the Florida and St. Maries strains of A. marginale were sequenced using three different technologies: 454, Illumina, and Ion Torrent. Roche's 454 technology provided the longest reads, as expected (Table 1). Interestingly, this technology also yielded the highest percentage of A. marginale reads at 37.1%. Although the Illumina platform had the lowest percentage of A. marginale reads (4.7% for the Florida strain), this was compensated by depth and was sufficient for quantitative analysis. The use of different platforms allowed us to address some of the challenges of working with obligate intracellular pathogens; as these microbes are dependent on their eukaryotic host cells, RNA samples are significantly contaminated by host transcripts, and RNA preps have been shown to be biased [25]. Our results were corroborated with the different platforms.
Transcriptome analysis allowed us to identify putative transcription start sites (TSS) for both strains as a byproduct of our study. Seventy putative TSSs were found in the Florida strain and 109 were found for the St. Maries strain (Additional files 3 and 4). The majority of these TSSs are present in both lists, the larger number of high confidence TSSs found in the St. Maries strain can be attributed to the deeper coverage obtained for this strain. Most of the Genome wide comparison of the non-transmissible Florida strain (red) with the efficiently transmitted St. Maries (green) strain produced 9609 SNPs. From this list we subtracted SNPs that encode for synonymous changes, leaving two types of SNPs that were further characterized: those that resulted in non-synonymous (NS) amino acid changes within ORFs and SNPs located in putative promoter regions. Comparison of these SNPs with genome sequences of three tick transmissible strains was then performed. SNPs that consistently segregated with phenotype were retained. The remaining differences were then targeted sequenced in two additional efficiently transmissible strains. B. A total of 9609 SNPs were found between the transmissible St. Maries and the nontransmissible Florida strain (SNPs). This comparison found 4498 nonsynonymous SNPs (represented in black), 1630 SNPs found within putative promoter regions (shown in dark grey) and synonymous SNPs (shown in light gray). Whole genome comparison with three transmissible strains allowed removal of 4127 non-synonymous SNPs and 1568 promoter SNPs from further consideration. Finally, Targeted sequencing in additional transmissible strains of 241 nonsynonymous and 62 promoter SNPs allowed retention of 35 NS and 14 promoter SNPs as candidate SNPs involved in tick transmission. putative 5 0 untranslated regions (UTR) found were longer than 40 bp in both strains (63.3% in St. Maries and 48.6% in Florida). Fewer putative 5 0 UTRs were found to be smaller than 40 bp, and the minority were found within the predicted open reading frame (ORF) ( Table 2). The 5 0 UTRs found within annotated CDS suggest that the predictions for these genes were inaccurate and an adjustment in annotation is required. We also identified 70 high confidence operon structures that involved 292 different genes (Additional file 5). Finally, six regions with no previous annotation were found to display high transcriptional activity ( Table 3). The six regions showed transcriptional activity in both strains, with two of these newly identified transcripts showing significant differential transcription between the strains. These regions are shown in Table 3 and Figure 4, and are further discussed in the next section.

Transcriptome comparison identifies transcriptional differences between strains with contrasting transmission phenotypes
After normalization, the distributions of the expression values across replicates were compared before evaluating changes in transcription (Additional file 6). Comparing the transcriptomes of the highly transmissible St. Maries strain with the poorly transmissible Florida strain produced a list of 14 genes that are significantly differentially transcribed using our criteria (see Methods) and across replicates ( Figure 3B, Figure 5, and Additional file 7). Genes that were found to have a lower transcription level in the poorly transmitted Florida strain are of particular interest considering the examined trait. Significant fold change differences that were constant across replicates ranged from 3.5 to 413.0 ( Figure 5, Figure 3B). Of the 10 genes that had significantly low (or absent) transcriptional activity in the Florida strain (Table 4), only one is annotated with a predicted localization: gene AMF_878, coding for outer membrane protein 4 (OMP4). The other nine genes are annotated as hypothetical proteins. Three of these had no mapped reads in any of the different sequencing technologies: AMF_431, AMF_432 and AMF_433. An additional two genes, AMF_429 and AMF_430, which appear to be arranged in an operon with AMF_431-3 (based on reads mapped to the St. Maries genome), are also significantly differentially transcribed ( Figure 3B). RNA-seq analysis of the promoter candidates identified by comparative genomics confirmed the RT-PCR results ( Figure 3A). Examination of candidates carrying non-synonymous SNPs found significant differential transcription of two genes; AMF_793 and AMF_1026 (shown in Table 4, along with differentially transcribed promoter candidates AMF_474   and AMF_553). These genes have a lower transcription level in the Florida strain by fold changes of 3.5 and 1.5 respectively (p < 1E-10). Finally, two newly identified regions were differentially transcribed. The regions between bp 336042 and 336685 and bp 1084944 and 1085520 in the St. Maries genome (Table 3) were up-regulated by fold changes of 23.6 and 6.9, respectively (p < 1E-10). These regions are shown in detail in Figure 4. In order to determine if these newly identified transcripts could indicate the presence of genes, we searched for ORFs that would overlap these regions. Only two regions: from bp 393765 to 394740 and 1084944 to 1085520 contained ORFs that would span the uninterrupted transcript. Comparisons of replicates were performed in order to account for variation of transcription values within a strain; importantly, the genes that were consistently differentially transcribed across replicates were found to be homogeneously transcribed when strain replicates were compared to each other.

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 Figure 3B, 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 upregulated 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. Table 4 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 Table 4; genes with candidate SNPs and differential transcription are shown in the middle of Table 4. 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 (Table 4). 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 (Table 4). 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 Gen-Bank 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.

Gene characteristics/bioinformatics
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). Twentythree genes had significant scores for transmembrane domain predictions but did not contain signal peptides. Analysis of non-synonymous SNPs using the SIFT algorithm [26] predicted eleven to be deleterious; these substitutions are reported in Additional file 8.

Discussion
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 [27]. 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 [21]. 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 [28].
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 nonsynonymous 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 file 8). 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 [29]. 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 [30], 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 [31]. 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 [32]. 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 [33]. 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 [34]. 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 [11]. 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 [35], no mobile elements or intact prophages have been identified in A. marginale [11].
Typically, pathogenic bacteria that cycle between arthropod and mammalian hosts modify their transcriptional profiles to adapt to these different environments [36]. 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][38][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 [37]. 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 [25]. 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 [40]. The deepest analysis generated 854,242 reads that mapped to the 1.23 Mb Chlamydia pneumonia genome [39]; 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 [40]. 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 0 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 0 UTRs were longer than 40 bp in both strains. Previous studies of TSSs have shown that only a very small portion of 5 0 UTRs are longer than 40 bp in bacteria [41,42]. As 5 0 UTRs have been involved in regulation processes in bacteria, further investigation of these elements might reveal translational and transcriptional roles [43]. 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 RNAseq 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 [44].

Conclusions
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 transmissionblocking vaccines.

Ethics statement
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.

Comparative genomics
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 [54] was used to compare as previously described [8] 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 file 1.

Targeted sequencing
The remaining SNPs were examined via targeted sequencing of the South Idaho, EMΦ and 6DE strains [50,55]. Primers were designed by aligning the SNPcontaining 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 [52]. The steady state, single copy gene msp5 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 [56]. 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 [24]; 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.

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 [57]. 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 [60]. 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 file 7).
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 [42].
The relative performances noted in Table 1 for the different sequencing technologies should not be directly compared, as this study was not designed to compare these platforms. As has been noted [25], 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. [42]: 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 0 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 [61]. 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 [62]. Transmembrane domains were predicted using two distinct algorithms: TMpred and Dense Alignment Surface (DAS) methods [63]; only genes with transmembrane domains predicted by both algorithms were reported. The "Sorting Tolerant From Intolerant" (SIFT) algorithm [26] 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).