Comparative genomics and transcriptomics of trait-gene association

Pierl6 et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13


BM ics
Genomics


Comparative genomics and transcriptomics of

trait-gene association

Sebastian Aguilar Pierl6*, Michael J Dark 23, Dani Dahmen Guy H Palmer' and Kelly A Brayton"


Abstract
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.
Keywords: Bacteria, Rickettsia, SNP, RNA-seq, Anaplasma


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

* Correspondence saguilar@vetmed wsu edu; kbrayton@vetmed wsu edu
Program in Genomics, Department of Veterinary Microbiology and
Pathology, Paul G Allen School for Global Animal Health, Washington State
University, Pullman, WA 99164-7040, USA
Full list of author information is available at the end of the article


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 pres-
ence of genes across species would suggest functional simi-
larity 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


S 2012 Pierle et al., licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative
Biol led Central 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.






Pierle et aL BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


different strains. Those bacterial species with a high num-
ber of accessory genes are termed "open" pangenomes,
whereas those lacking strain specific genes are identified as
"closed" pangenomes. While the "openness" of the pangen-
ome is an obvious marker of diversity, sequence heterogen-
eity 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 nucleo-
tide 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 pan-
genomes, 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 do-
mestic and wild ruminants, with a small genome of 1.2 Mb
for which the sequence of multiple strains has been deter-
mined [8,11,12]. No strain-specific genes and no plasmids
were found among sensu strict strains after sequencing of
five strains [8,11]. In contrast, a high degree of allelic diver-
sity 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 be-
tween any given pair. The high degree of gene content con-
servation suggests that phenotypic differences observed in
A. marginale must be due to small polymorphisms be-
tween strains rather than whole gene insertions or dele-
tions. 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, EMo, 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 transmis-
sion with the St. Maries strain [17,19,20].
Due to the complete gene content conservation, differ-
ences in transmission efficiency in A. marginale are likely
to be ascribed to sequence variation producing variant pro-
teins or affecting gene transcription. Indeed, precedence is
seen in bacterial pathogens, where SNPs have been discov-
ered that provide a selective advantage in host colonization
[21]. We combined two genomic sequencing approaches
in order to find SNPs and transcriptional changes that seg-
regate 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 sub-
stitutions within genes, as well as SNPs located within pu-
tative promoter regions. Each SNP on the resulting list was
evaluated through comparative genomics in three effi-
ciently 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 candi-
date genes, consisting of those containing candidate SNPs
or located downstream of putative promoter SNPs. Tran-
scriptional analysis of candidate genes by RT-PCR revealed
genes that were differentially transcribed in strains with
distinctly different transmission efficiencies. To find add-
itional transcriptional changes related to the phenotype
of interest, we performed a genome wide transcriptome
comparison using RNA-seq technology. Total mRNA po-
pulations from two A. marginale strains with different
transmission capabilities were sequenced using three differ-
ent 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 intracel-
lular bacterial transcriptome during the pathogen's natural
course of infection.


Results
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 gen-
ome (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 pur-
poses of this study, putative promoters were defined as
intergenic regions immediately 5' 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 add-
itional efficiently transmitted strains, 6DE and EM0, 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).


Page 2 of 15






Pierl6 et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


A
Florida vs St. Maries


9609 SNPs



Syn NS + Promoter
SNPs SNPs SNPs

+
Whole genome comparisons
in 3 strains

+
Targeted sequencing
in 2 transmissible strains


-I *4 I.


StM-FL Whole
SNPs genome
comparisons


Targeted Candidates
sequencing


Figure 1 SNPs segregated with transmission status through
whole genome comparison and targeted sequencing. A.
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 non-
transmissible Florida strain (SNPs). This comparison found 4498 non-
synonymous 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 non-
synonymous and 62 promoter SNPs allowed retention of 35 NS and
14 promoter SNPs as candidate SNPs involved in tick transmission.


Altogether, comparative genomics identified 29 candidate
genes.

Transcription analysis of candidate genes
These 29 genes with SNPs in their coding regions or in
their putative promoter regions were analyzed for transcrip-
tional 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 candi-
date promoter regions, the relative expression ratio was
analyzed from two separate infections using mspS 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 differen-
tially 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 trans-
missible 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 candi-
dates, 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%. Al-
though 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 quanti-
tative analysis. The use of different platforms allowed us
to address some of the challenges of working with obli-
gate 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 by-
product 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 confi-
dence TSSs found in the St. Maries strain can be attributed
to the deeper coverage obtained for this strain. Most of the


Page 3 of 15


B
* 12,000



"8,000



4,000


7-






Pierle et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


putative 5' 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' UTRs were found to be smal-
ler than 40 bp, and the minority were found within the
predicted open reading frame (ORF) (Table 2). The 5'
UTRs found within annotated CDS suggest that the pre-
dictions for these genes were inaccurate and an adjustment
in annotation is required. We also identified 70 high confi-
dence operon structures that involved 292 different genes
(Additional file 5). Finally, six regions with no previous an-
notation 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 differen-
tially 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 particu-
lar interest considering the examined trait. Significant fold
change differences that were constant across replicates ran-
ged 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 tech-
nologies: AMF_431, AMF_432 and AMF_433. An add-
itional 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 signifi-
cantly differentially transcribed (Figure 3B). RNA-seq ana-
lysis 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


Page 4 of 15


600000
Figure 2 Location of candidate SNPs on the Florida strain genome. This circular representation of the Florida genome shows in light blue
annotated CDSs; outer circle represetsCDSs on the wastrand, inner circle epesents the reverse strand, in grey the 9609 SNPs found
between the St. Maries and the Florida strain genomes. The elements in light green are miscellaneous features annotated in the genome. In the
inner most circle 49 candidate SNPs found through comparative genomics are shown. Red bars show the position of candidate non-synonymous
SNPs within CDSs. Dark blue bars show candidate SNPs found within putative promoter regions.







Pierl6 et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


Table 1 Reads mapped to A. marginale from three sequencing platforms
Platform Strain Total reads A. marginale reads % of reads mapped to A. marginale Average matched read length
Roche 454 STM' 726,051 269,730 37.1 381.7


Ion Torrent STM


1,018,447
1,004,747
2,043,607


STM 88,650,713
FL 81,507,967


1STM: St. Maries strain.
2FL: Florida strain.


Page 5 of 15


A
10





CD







0.1






B
1000


)1001


S10





0.1 -




Figure 3 RNA-seq and qPCR confirm trends in transcriptional changes between strains that differ in their tick transmission status. A.
Fold change in the transmissible St. Maries strain relative to the non-transmissible Florida strain for all promoter candidates expressed in log scale
10. Locus tags for all genes are given on the X axis. Blue bars show the results obtained after evaluating two biological replicates with RT-PCR.
Red bars show the fold change obtained using RNA-seq analysis for the promoter candidates across two biological replicates. The asterisk
indicates statistical at p < 0.05. B. Fold change in the transmissible St. Maries strain relative to the non-transmissible florida strain
expressed in logi. The top 18 differentially transcribed genes identified through RNA-seq across two replicates and two statistical tests and their
fold changes are shown. Red bars show results obtained with RNA-seq, blue bars show validation through qPCR. The asterisk indicates statistical
at p <0.01.


Illumina


326,440
295,629
577,284
4,604,993
3,845,853







Pierle et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


Table 2 Percentage of putative 5' UTRs according to
length
Strain 5' UTRs <40 bp 5' UTRs 40 bp 5' UTR within
predicted CDS'


St. Maries
Florida


63.30


24.28


'5' UTR within predicted CDS: in this column we list cases where transcript
mapping shows that the 5' UTR and TSS are found within the previously
predicted and annotated CDS indicating that the previous annotation was
incorrect.


and AMF_553). These genes have a lower transcription
level in the Florida strain by fold changes of 3.5 and 1.5 re-
spectively (p < IE-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 con-
sistently differentially transcribed across replicates were
found to be homogeneously transcribed when strain repli-
cates 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 con-
firm 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 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.


Gene characteristics/bioinformatics
Table 4 shows the 30 genes that were selected as candi-
dates. Genes that were found to be differentially tran-
scribed through RNA-seq and RT-PCR are shown on top
of Table 4; genes with candidate SNPs and differential tran-
scription 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 candi-
date 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 pro-
teins (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 down-
stream from promoter SNPs included one translation in-
hibitor (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 sur-
prising that more than half of the candidates were of un-
known 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 homo-
logs 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 Ehrli-
chia 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. Add-
itionally, hypothetical candidates AMF_1037, AMF_879,
AMF_401 and AMF_530 had no homologs in the Gen-
Bank database. These findings provide two mutually ex-
clusive scenarios: if a gene with homologs in the
aforementioned tick-transmitted organisms is responsible
for the trait of interest, this suggests a common mechan-
ism within a bacterial order or family. Alternatively, a gene


Table 3 Previously unannotated areas that exhibited high transcriptional activity


Region' Len
251313..251855 543
336042..336685 644
393765..394740 976
459343..459783 441
887245..887579 335
1084944..1085520 577


gth Identity through blastX
hypothetical protein AmarV 01231 [Anaplasma marginal str
n/a
DNA-binding protein HU [Anaplasma phagocytophilum HZ]


Gene before Gene after Fold change
ia] AM294 pepl AM259 thiD 0.7
AM380 AM382 23.6*
AM434 pdxJ AM435 1.1


hypothetical protein AmarM_02282 [Anaplasma marginal str. Mississippi] AM504


tRNA-Asr


hypothetical protein PseS9_19739 [Pseudomonas sp. S9] AM969 bioB AM973 pur
hypothetical protein AmarM_05569 [Anaplasma marginal str. Mississippi] AM1214 polA AM1216


base pair positions spanned by the newly identified regions in the St. Maries genome.
*Statistically significant fold changes are indicated with an asterisk.


Page 6 of 15







Pierl6 et al. BMC Genomics 2012, 13:669 Page 7 of 15
http://www.biomedcentral.com/1471-2164/13/669





A
I I I I I I I II I I
I I I I I I II I I
I I I II I I I

:.... : I: I..... ... .i i
AM382
I I I I <
I I I I I II I II I I
I I I I I I II I I I
















B
I I I I I
II I I IIII I II I

PUTATIVE CDS AM1216
AMI216
385000 L085100 L085200 L085300 L1085400 L085500 L1085600 10857

I I I I I I
I I I I
I I I I I I I I I I I I I I













Figure 4 Newly identified transcriptionally active regions of the genome. Mapping of cDNA reads to the A. morginale genome allowed us
to detect regions without previous annotation that exhibited transcriptional activity. A shows the region of the St. Maries genome that spans
from bp 336042 to 336685. Three different gene identification algorithms did not detect a CDS that would span the length of the transcript. The
top panel shows the six reading frames containing forty-five stop codons, shown as black bars. The bottom panel shows some of the mapped
cDNA reads in green and red (indicating direction of the read). The grey histogram under the reads represents depth (read height). This transcript
was up-regulated in the St. Maries strain by a fold change of 23.7 at p < 1E-10. B shows the region of the St. Maries genome that spans from bp
1084944 to 1085520. The newly identified gene is found between genes polA (not shown) and AMI216. One ORF on the leading strand seems
to span the length of this transcript and is shown as PUTATIVECDS in this figure. This new gene was found to be up-regulated in the
transmissible St. Maries strain by a fold change of 6.9 at p < 1E-10.







Pierl6 et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


Page 8 of 15


Genes
Figure 5 Whole genome comparison of transcriptional activity in the St. Maries and Florida strains. The RPKM values for 955 genes found


he Y axis. Features are a arranged from left to right as they appear
strain. RPKM values for the transmissible St. Maries strain are sh
ida strain are shown in light blue in the lower part of the graph
for ease of comparison; they do not represent negative values.


in the Florida strain genome of A. marginal. RPKM values are shown on tl
genome on the X axis. The normalized RPKM values were plotted for each
red in the upper part of the graph; numbers for the non-transmissible Flor
values for the Florida strain are plotted on the opposite side of the x axis
Ribosomal RNA (rRNA) genes were subtracted from this comparison.


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 do-
main predictions but did not contain signal peptides. Ana-
lysis 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 segre-
gate 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 transcrip-
tion 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 ana-
lysis 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 biol-
ogy. A single non-synonymous SNP in the envelope pro-
tein gene El 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 mono-
mannose 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 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


in the
)wn in
RPKM







Pierle et al. BMC Genomics 2012, 13:669 Page 9 of 15
http://www.biomedcentral.com/1471-2164/13/669




Table 4 Candidate genes involved in transmission phenotype segregated by polymorphisms and differential
transcription


Category

Differentially transcribed
genes















Differentially transcribed
genes w/SNPs
genescarrying candidate
SNPs





Genes carrying candidate
SNPs


FL' St.M1


AM579
AM579
AM580
AM576
AM574
AM1055
AM1165
AM1164
AM540
AM347
AM632


AMF 553 AM748


793 AM 1048
1026 AM1352
051 AM071
197 AM265
264 AM354
265 AM356
269 AM368
480 AM644


Product

Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Outer membrane protein
Hypothetical protein
Hypothetical protein
Ribosome-associated
inhibitor A
NADH Dehydrogenase
chain J
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
Hypothetical protein
DNA gyrase B


SNPs candidate
SNPs
1 0
1 0
3 0
22 0
15 0


SNP location

Gene
Gene
Gene
Gene
Gene
Gene
Gene
Gene
Gene/Prom
Gene
Promoter


Homologs2 B13


AP,
AP, ER
AP, ER
AP, ER
AP,
AP, ER


ER, ECh TM
, ECh, ECa
, ECh, ECa TM
, ECh, ECa TM/DS
ER, ECh TM
, ECh, ECa TM/SP
TM
AP TM/SP
TM
AP TM


R, ECh, ECa


Promoter AP, ER, ECh, ECa


77 2 Gene/Prom
18 5 Gene/Prom
5 1 Gene
19 1 Gene
21 2 Gene
65 1 Gene
43 1 Gene
14 1 Gene


AP, ER, ECh, ECa
AP, ER, ECh
AP, ER, ECh, ECa
AP, RB
ER, ECh, ECa
AP
RB
AP, ER, ECh, ECa, RE
RC, RR


TM


TM/SP
TM/DS
TM
TM
TM/DS


TM/DS
TM


ethical pr<
ethical pr<
ethical pr<
etical nr2


Hypothetical protel


AMF 762 AM1001 methionyl-tRNA synthetase


AMF 764 AM1005
AMF 824 AM 091


AMF 893 AM 1183


AMF 1037 AM345


aspartokinase
D-Ala-D-Ala
carboxypeptidase
lipoprotein-releasing
transmembrane protel
Hypothetical protein


Gene AP, ER, ECh, ECa, RB,
RC, RR
Gene AP, ER, ECh, ECa
Gene AP, ER, ECh, ECa, RB,
RC, RR
Gene ER, ECh, ECa


1FL: Locus tag in the Florida strain St.M: locus tag in the St. Maries strain.
2AP: Anaplasma phagocytophilum, ER: Ehrlichia ruminantium, ECh: Ehrlichia chaffeensis, ECa: Ehrlichia canis, RB: Rickettsia bellii, RC: Rickettsia conorii, RR: Rickettsia
rickettsia.
3BI: Bioinformatics TM: Transmembrane domain, SP: Signal peptide, Prom: promoter, DS: Deleterious substitution.


and transmembrane domains. The presence of signal
peptides and transmembrane domains implies mem-
brane localization of the proteins, and thus, these pro-
teins 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 deleteri-
ous substitutions with a total of three non-synonymous
SNPs. This gene was also found to be up-regulated in


AM712
AM689
AM742
AM823
AM919


Gene
Gene
Gene
Gene
Gene


AP, E
AP, EF
AP, ER,


AP
, ECh, ECa
R, ECh,ECa
ECh, ECa, RE
KC, RR


TM/DS


TM/DS
TM/SP


TM/DS






Pierle et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


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 con-
sistently 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 pro-
moter regions were found to be differentially tran-
scribed. 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 re-
spiratory complex I, a conserved proton pumping
NADH:ubiquinone oxidoreductase in bacteria [29]. An-
other closely associated gene from this complex, nuoL,
has been found to be up-regulated in the related organ-
ism Rickettsia conorii while dealing with osmotic stress
[30], suggesting that enhancement of NADH dehydro-
genase 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 ribo-
some, and preventing the formation of initiation com-
plexes by preventing the binding of mRNA and initiator
fMet-tRNA to the ribosome [33]. When temperature
levels return to 37oC, initiation of protein synthesis over-
comes the PY inhibition as tRNA compete more effect-
ively 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 observa-
tions 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 acti-
vity that the strategy might have overlooked; the tran-
scriptome 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 sig-
nificantly differentially transcribed. Of these 10 genes
only one had a predicted localization: AMF_878 corre-
sponds 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 con-
served domains for tail and head/tail connector phage
proteins, with the highest similarity found to phage
proteins from Wolbachia spp., a related bacterial sym-
biont of arthropods. Although this could open interest-
ing possibilities, as phages play an important role for
Wolbachia spp. within the arthropod host [35], no mo-
bile elements or intact prophages have been identified
in A. marginale [11].
Typically, pathogenic bacteria that cycle between
arthropod and mammalian hosts modify their transcrip-
tional 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 transcrip-
tomic 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 lim-
ited 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 tran-
scriptional response, with less than 10 genes changing
more than 3-fold in expression level [37]. This could in-
dicate 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 obli-
gate intracellular habitat of pathogens in this Order may


Page 10 of 15






Pierle et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


offer such a stable environment that the necessity for
gene regulation is much less than that of facultative
intracellular pathogens. Our study searched for tran-
scriptional 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 consist-
ent 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. mar-
ginale 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 popula-
tions daunting. Few attempts have been made at describ-
ing 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 patho-
gen A. phagocytophilum [40]. The deepest analysis gener-
ated 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. marginal'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 sam-
ples 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 tran-
scriptional analysis.
Analyzing transcriptional profiles with RNA-seq allows
us to evaluate "snapshots" in time of bacterial transcrip-
tomes; therefore, it is essential to generate data from more
than one replicate to provide a broader more reliable pic-
ture of transcriptional changes. The depth and reproduci-
bility of this RNA-seq data set allowed for mapping of the
physical structure of the A. marginale transcriptome; in-
cluding 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 [43]. Additionally, mapping of
transcriptional data allowed us to define 70 putative op-
eron 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 high-
light the complexity of transcriptional regulation in bac-
teria [44].

Conclusions
This study takes advantage of the high interstrain diver-
sity 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 transmis-
sion 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. Al-
though the successful transformation of A. marginale
has been achieved [45,46], stable targeted gene replace-
ment 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 strat-
egies for tick-borne pathogens, such as transmission-
blocking vaccines.

Methods
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.

Strains
The Florida, St. Maries, Virginia, Puerto Rico, South
Idaho, EM() and 6DE strains used in this study have been


Page 11 of 15






Pierle et al. BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


described in detail elsewhere [47-51]. The St. Maries,
Virginia, Puerto Rico, South Idaho, EMO 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, in-
cluding the Reynolds Creek stock [15,18,19].

Comparative genomics
The accession numbers for the strains used are: St.
Maries: CP000030.1, Florida: CP001079.1, Viriginia:
ABOROOOOOOOO.1, Puerto Rico: ABOQ00000000.1, South
Idaho: AFMYOOOO00000000.1. MUMmer v3.1 [54] was used to
compare as previously described [8] to compare the Flor-
ida 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 se-
quencing of the South Idaho, EMO 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 elimi-
nated candidates as described above. All candidate SNPs
were resequenced in the Florida strain, to verify the ori-
ginal 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 tran-
scription 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 re-
flect transcript levels based on reverse transcription effi-
ciency [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 stat-
istical 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 bio-
logical 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 deter-
mined 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. Infec-
tion 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 pro-
cessed for 454 and Ion Torrent technologies, probes for
bacterial ribosomal RNAs from the Ribominus kit (Invi-
trogen) were added during the subtractive hybridization
procedure. For samples processed for Illumina, the Du-
plex--Specific thermostable nuclease (DSN) normalization
protocol was applied. Data was processed using CLC Gen-
omics Workbench (CLC Bio). Mapping parameters were
adjusted to map a maximum number of reads to the refer-
ence bacterial genomes. The distribution of the expression
values for all samples was analyzed and compared.
Normalization by quantiles was applied to adjust the dis-
tributions for further comparison. Fold changes with re-
spect 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 propor-
tions [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


Page 12 of 15







Pierle et aL BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


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 statis-
tical significance [60]. Genes whose log2 fold change was
above and below 2 and -2, respectively, and whose -logl0
p-value was above 10 in both replicate comparisons and
under both statistical tests were selected for further evalu-
ation (Additional file 7).
Areas of the genome that were not previously annotated
and showed >0.5 coverage (average sequence data cover-
age depth) were reported when reads were unambiguously
mapped to the A. marginale genome [42].
The relative performances noted in Table 1 for the differ-
ent sequencing technologies should not be directly com-
pared, as this study was not designed to compare these
platforms. As has been noted [25], different library prepara-
tions and sequencing technologies favor recovery of differ-
ent 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 up-
stream of the open reading frame, we designated the point
at which coverage dropped to 0 as the putative transcrip-
tional 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 as-
sociation", examines the annotation of the currently avail-
able 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 substitu-
tions, and was used to predict if substitutions in the candi-
date 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).



Additional files

Additional file 1: Nucleotide polymorphisms between the St.
Maries and Florida strains. SNPs between the St Maries and Florida
strains are listed here together with the nucleotides reported for all the
additional reported nucleotides
Additional file 2: Absolute expression values of candidate genes.
Gene identify cations are provided on the x axis and the copy number per
ml of blood on the y axis The black bars represent the numbers
obtained for the Florida strain and the white bars the numbers for the St
Maries strain Transcription of all candidate genes is shown together with
the calibrator MSP5
Additional file 3: Mapping of putative TSS and 5' UTRs length in
the St. Maries strain. The location and length of 5' UTRs in the St
Maries strain are reported
Additional file 4: Mapping of putative TSS and 5' UTRs length in
the Florida strain. The location and length of 5' UTRs in the Florida
strain are reported
Additional file 5: Operon strucutres found through transcriptome
sequencing in the St. Maries strain. Genes involved in the different
operon structures are reported
Additional file 6: Distribution of the normalized expression values
of all replicates analyzed in this study with RNA-seq. The distribution
of the normalized RPKM values for all replicates is plotted in a box plot
RNA-SeqFL1 and RNA-SeqFL2 designate distributions for Florida strain
replicates 1 and 2 respectively RNA-Seq STM 1 and RNA-Seq STM 2
designate RPKM distributions for St Maries replicates 1 and 2 respectively
The distributions allow for comparisons
Additional file 7: A. marginale genes arranged along dimensions of
biological and statistical significance. A volcano plot shows the
relationship between the p-values of a statistical test and the magnitude
of the difference in expression values On the y axis the negative loglO
p-values are plotted On the x-axis the log 2 values of the fold changes
seen in whole transcriptome comparison The red lines highlight the
cutoffs for genes that were analyzed further Only genes populating the
upper right and left quadrants of the plot under two different statistical
tests (Kal's and Baggerly's) were chosen This plot shows results obtained
for Kal's test
Additional file 8: Candidate non-synonymous changes predicted to
be deleterious by the SIFT algorithm. Non-synonymous changes
predicted to be deleterious by the SIFT algorithm found between the St
Maries and Florida strains are reported


Abbreviations
CDS Coding DNA Sequence; NS Non-Synonymous; ORF Open Reading
Frame; RPKM Reads Per Kilobase per Million mapped reads; SNP Single-
Nucleotide Polymorphism; UTR UnTrasnIated Region


Page 13 of 15








Pierle et aL BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


Competing interests
The authors declare that they have no competing interests


Authors' contributions
SAP, MJD, GHP, KAB conceived the experiments; SAP, MJD, DD performed
the experiments; SAP, MJD, DD, GHP, KAB analyzed the data; SAP, MJD, GHP
KAB wrote and edited the manuscript All authors read and approved the
final manuscript


Acknowledgments
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

Author details
Program in Genomics, Department of Veterinary Microbiology and
Pathology, Paul G Allen School for Global Animal Health, Washington State
University, Pullman, WA 99164-7040, USA 2Department of Infectious Diseases
and Pathobiology, College of Veterinary Medicine, University of Florida,
Gainesville, FL 32611-0880, USA Emerging Pathogens Institute, University of
Florida, Gainesville, FL 32611-0880, USA

Received: 7 August 2012 Accepted: 16 November 2012
Published: 26 November 2012


References
1 Bork P, Dandekar T, Diaz-Lazcoz Y, Eisenhaber F, Huynen M, Yuan Y'
Predicting function: from genes to genomes and back. Mol Biol 1998,
283(4) 707-725
2 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
3 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
4 Makarova KS, Wolf YI, Koonin EV Potential genomic determinants of
hyperthermophily. Trends Genet 2003, 19(4)'172-176
5 Huynen MA, Bork P Measuring genome evolution. Proc Nat/ Acod Sc USA
1998, 95(11)'5849-5856
6 Pellegrini M, Marcotte EM, Thompson MJ, Eisenberg D, Yeates TO'
Assigning protein functions by comparative genome analysis: protein
phylogenetic profiles. Proc Nat/ Acod Sc U SA 1999, 96(8)4285-4288
7 Tettelin H, Masignani V, Cieslewicz MJ, Donati C, Medini D, Ward NL,
Angiuoli SV, Crabtree J, Jones AL, Durkin AS, et o/ Genome analysis of
multiple pathogenic isolates of Streptococcus agalactiae: implications
for the microbial "pan-genome". Proc Nat/ Acod Sci USA 2005,
102(39)13950-13955
8 Dark MJ, Herndon DR, Kappmeyer LS, Gonzales MP, Nordeen E, Palmer GH,
Knowles DP Jr, Brayton KA Conservation in the face of diversity:
multistrain analysis of an intracellular bacterium. BMC Genomics 2009,
10:16
9 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 Nat Acod Sci US A 2004,101 (37) 13536-13541
10 Stratilo CW, Lewis CT, Bryden L, Mulvey MR, Bader D Single-nucleotide
repeat analysis for subtyping Bacillus anthracis isolates. J Clin Microbio/
2006, 44(3)777-782
11 Brayton KA, Kappmeyer LS, Herndon DR, Dark MJ, Tibbals DL, Palmer GH,
McGuire TC, Knowles DP Jr Complete genome sequencing of Anaplasma
marginale reveals that the surface is skewed to two superfamilies of
outer membrane proteins. Proc Nat/ Acod Sc USA 2005, 102(3)'844-849
12 Dark MJ, Al-Khedery B, Barbet AF Multistrain genome analysis identifies
candidate vaccine antigens of Anaplasma marginale. Vaccine 2011,
29(31) 4923-4932
13 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


14 Futse JE, Brayton KA, Dark MJ, Knowles DP Jr, Palmer GH Superinfection as
a driver of genomic diversification in antigenically variant pathogens.
Proc Nat/ Acod Sc U S A 2008, 105(6)'2123-2127
15 Futse JE, Ueti MW, Knowles DP Jr, Palmer GH Transmission of Anaplasma
marginal by Boophilus microplus: retention of vector competence in
the absence of vector-pathogen interaction. J Clin Microbiol 2003,
41(8)3829-3834
16 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
17 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
18 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
19 Ueti MW, Reagan JO Jr, Knowles DP Jr, 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
20 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
21 Sokurenko EV, Hasty DL, Dykhuizen DE Pathoadaptive mutations: gene
loss and variation in bacterial pathogens. Trends Microbiol 1999,
7(5)'191-195
22 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
23 Visser ES, McGuire TC, Palmer GH, Davis WC, Shkap V, Pipano E, Knowles DP
Jr The Anaplasma marginale msp5 gene encodes a 19-kilodalton
protein conserved in all recognized Anaplasma species. Infect Immun
1992, 60(12)'5139-5144
24 Willems E, Leyns L, Vandesompele J Standardization of real-time PCR
gene expression data from independent biological replicates. Anoal
Biochem 2008, 379(1 )127-129
25 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
26 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
27 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
28 Stermann M, Sedlacek L, Maass S, Bange FC A promoter mutation
causes differential nitrate reductase activity of Mycobacterium
tuberculosis and Mycobacterium bovis. J Bocterio 2004,
186(9) 2856-2861
29 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
30 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
31 Munderloh UG, Kurtti TJ Cellular and molecular interrelationships
between ticks and prokaryotic tick-borne pathogens. Annu Rev Entomol
1995, 40:221-243
32 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
33 Vila-Sanjurjo A, Schuwirth BS, Hau CW, Cate JH Structural basis for the
control of translation initiation during stress. Nat Struct Mol Bio/ 2004,
11(11)'1054-1059
34 McDade JE, Newhouse VF Natural history of Rickettsia rickettsii. Annu Rev
Microbio 1986, 40:287-309


Page 14 of 15








Pierle et aL BMC Genomics 2012, 13:669
http://www.biomedcentral.com/1471-2164/13/669


35 Bordenstein SR, Marshall ML, Fry AJ, Kim U, Wernegreen JJ The tripartite
associations between bacteriophage, Wolbachia, and arthropods. PLoS
pathogens 2006, 2(5)'e43
36 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
37 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
38 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
39 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
40 Mastronunzio JE, Kurscheid S, Fikrig E Post-genomic analyses reveal
development of infectious Anaplasma phagocytophilum during
transmission from ticks to mice. J Bacterio/ 2012, 194(9) 2238-2247
41 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
42 Passalacqua KD, Varadarajan A, Ondov BD, Okou DT, Zwick ME, Bergman
NH Structure and complexity of a bacterial transcriptome. J Bacterio/
2009, 191(10)'3203-3211
43 Naville M, Gautheret D Transcription attenuation in bacteria: theme and
variations. Brief Funct Genomic Proteomic 2009, 8(6)482-492
44 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
45 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
46 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. App Environ
Microbio 2011, 77(1) 330-334
47 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
48 Kuttler KL, Winward LD Serologic comparisons of 4 Anaplasma isolates as
measured by the complement-fixation test. Vet Microbiol 1984,
9(2) 181-186
49 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
50 Palmer GH, Knowles DP Jr, Rodriguez JIL, 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 Microbio 2004, 42(11 )5381-5384
51 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
52 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
53 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
54 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
55 Rodriguez JIL, Palmer GH, Knowles DP Jr, Brayton KA Distinctly different
msp2 pseudogene repertoires in Anaplasma marginale strains that are
capable of superinfection. Gene 2005, 361:127-132
56 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


57 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
58 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
59 Kal AJ, van Zonneveld AJ, Benes V, van den Berg M, Koerkamp MG,
Albermann K, Strack N, Ruijter JM, Richter A, Dujon B, et o/ 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
60 Cui X, Churchill GA Statistical tests for differential expression in cDNA
microarray experiments. Genome Bio/ 2003, 4(4) 210
61 Tabor HK, Risch NJ, Myers RM Candidate-gene approaches for studying
complex genetic traits: practical considerations. Not Rev Gener 2002,
3(5)'391-397
62 Emanuelsson 0, 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
63 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

doi:10.1186/1471-2164-13-669
Cite this article as: Pierle et al Comparative genomics and
transcriptomics of trait-gene association. BMC Genomics 2012 13'669


Page 15 of 15


Submit your next manuscript to BioMed Central
and take full advantage of:

* Convenient online submission
* Thorough peer review
* No space constraints or color figure charges
* Immediate publication on acceptance
* Inclusion in PubMed, CAS, Scopus and Google Scholar
* Research which is freely available for redistribution

Submit your manuscript at nt
www.biomedcentral.com/submit 0 I IlIaId Central

PAGE 1

9 8 7 6 5 4 3 2 1 0 Copy #(log 10 )/ml of blood


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E87BPAHMH_A4OI4T INGEST_TIME 2013-03-05T20:09:41Z PACKAGE AA00013515_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

Supplementary Table S 5 : Candidate non synonymous changes predicted to be deleterious by the SIFT algorithm Gene Description Position in the Florida genome Amino Acid in Florida Amino acid in St. Maries AMF_264 hypothetical protein 299040 Glutamine Arginine AMF_264 hypothetical protein 299131 Serine Glycine AMF_269 hypothetical protein 321307 Aspartic acid Cysteine AMF_430 hypothetical protein 528313 Glycine Aspartic acid AMF_547 hypothetical protein 680650 Cysteine Tyrosine AMF_762 methionyl tRNA synthetase 911553 Serine Arginine AMF_764 aspartate kinase 913078 Methionine Leucine AMF_893 lipoprotein releasing system transmembrane protein 1061636 Histidine Tyrosine AMF_1026 hypothetical protein 1195613 Alanine Valine AMF_1026 hypothetical protein 1195626 Alanine Threonine AMF_1026 hypothetical protein 1195634 Valine Alanine



PAGE 1

RESEARCHARTICLEOpenAccessComparativegenomicsandtranscriptomicsof trait-geneassociationSebastinAguilarPierl1*,MichaelJDark2,3,DaniDahmen1,GuyHPalmer1andKellyABrayton1*AbstractBackground: TheOrder Rickettsiales includesimportanttick-bornepathogens,from Rickettsiarickettsii ,whichcauses RockyMountainspottedfever,to Anaplasmamarginale ,themostprevalentvector-bornepathogenofcattle. AlthoughmostpathogensinthisOrderaretransmittedbyarthropodvectors,littleisknownaboutthemicrobial determinantsoftransmission. A.marginale providesuniquetoolsforstudyingthedeterminantsoftransmission, withmultiplestrainsequencesavailablethatdisplaydistinctandreproducibletransmissionphenotypes.Theclosed core A.marginale genomesuggeststhatanyphenotypicdifferencesareduetosinglenucleotidepolymorphisms (SNPs).WecombinedDNA/RNAcomparativegenomicapproachesusingstrainswithdifferentticktransmission phenotypesandidentifiedgenesthatsegregatewithtransmissibility. Results: ComparisonofsevenstrainswithdifferenttransmissionphenotypesgeneratedalistofSNPsaffecting18 genesandninepromoters.TranscriptionalanalysisfoundtwocandidategenesdownstreamfrompromoterSNPs thatweredifferentiallytranscribed.TocorroboratethecomparativegenomicsapproachweusedthreeRNA-seq platformstoanalyzethetranscriptomesfromtwo A.marginale strainswithdifferenttransmissionphenotypes. RNA-seqanalysisconfirmedthecomparativegenomicsdataandfound10additionalgeneswhosetranscription betweenstrainswithdistincttransmissionefficiencieswassignificantlydifferent.Sixregionsofthegenomethat containednoannotationwerefoundtobetranscriptionallyactive,andtwoofthesenewlyidentifiedtranscripts weredifferentiallytranscribed. Conclusions: Thisapproachidentified30genesandtwonoveltranscriptspotentiallyinvolvedinticktransmission. Wedescribethetranscriptomeofanobligateintracellularbacteriumindepth,whileemployingmassiveparallel sequencingtodissectanimportanttraitinbacterialpathogenesis. Keywords: Bacteria,Rickettsia,SNP,RNA-seq,AnaplasmaBackgroundTheongoingrevolutioningenomesequencinghasenabled ever-increasingsequencegenerationatanever-decreasing cost.Thegrowingavailabilityoffullysequencedgenomes offersnewopportunitiestoidentifyrelationshipsbetween genotypeandphenotype,oneofthemajorgoalsofthe genomicsera.Comparativegenomicswerefirstintroduced asatooltopredicttrait-geneassociationsin1998while tryingtodefinespecies-specificfeaturesof Helicobacter pylori [1].Thisapproachhasbeenusedtopredictgenomic determinantsforwell-knownphenotypes,including hyperthermophily,flagellarmotilityandpiliassembly[2-4]. Thesestudiessharetheprinciplethatspecieswithsimilar phenotypesarelikelytoutilizeorthologousgenesinthe involvedbiologicalprocess.Thus,thesimultaneouspresenceofgenesacrossspecieswouldsuggestfunctionalsimilarityamongencodedproteins[5,6].Whilethesestudies illustratetheadvantagesandapplicabilityofthisprinciple, theyaredependentonpreviousknowledgeofthegenetic determinantsofaspecifictrait. Thechallengeofassociatinggeneswithphenotypeshas beenhighlightedbythedevelopmentofthepangenome conceptandtheabundanceofintraspeciesdiversitythat hasbeenrevealed.Thepangenomeofabacterialspecies encompassesthesumofthegeneticrepertoirefoundinall strains[7].Thus,itconsistsofthecoregenomefoundin allthestrainsplusthe “ accessory ” genesuniquetothe *Correspondence: saguilar@vetmed.wsu.edu ; kbrayton@vetmed.wsu.edu1PrograminGenomics,DepartmentofVeterinaryMicrobiologyand Pathology,PaulG.AllenSchoolforGlobalAnimalHealth,WashingtonState University,Pullman,WA99164-7040,USA Fulllistofauthorinformationisavailableattheendofthearticle 2012Pierletal.;licenseeBioMedCentralLtd.ThisisanOpenAccessarticledistributedunderthetermsoftheCreative CommonsAttributionLicense(http://creativecommons.org/licenses/by/2.0),whichpermitsunrestricteduse,distribution,and reproductioninanymedium,providedtheoriginalworkisproperlycited.Pierl etal.BMCGenomics 2012, 13 :669 http://www.biomedcentral.com/1471-2164/13

PAGE 2

differentstrains.Thosebacterialspecieswithahighnumberofaccessorygenesaretermed “ open ” pangenomes, whereasthoselackingstrainspecificgenesareidentifiedas “ closed ” pangenomes.Whilethe “ openness ” ofthepangenomeisanobviousmarkerofdiversity,sequenceheterogeneitywithinthecoregenesethasalsobeenshowntobe relevanttonaturalgeneticvariation[8-10].Whenseveral strainsof Streptococcusagalactiae werecomparedtothe 2603VRstrain,99.2%ofthetotaldetectedsinglenucleotidepolymorphisms(SNPs)wereuniquetoonestrain, whilenonewerecommontoallstrains.Asimilarscenario wasfoundbetweenthreestrainsof Bacillusanthracis whereallSNPswereuniquetoonestrain.Asthesetwo organisms,areclassicalexamplesofopenandclosedpangenomes,respectively,thissuggeststhattheSNPprofileof abacterialspeciescanbeopenregardlessofhow “ locked ” theircoresare. Anexampleofanorganismwithaclosedcoregenome andahighdegreeofinterstraindiversityis Anaplasma marginale ,anobligateintracellularpathogenofbothdomesticandwildruminants,withasmallgenomeof1.2Mb forwhichthesequenceofmultiplestrainshasbeendetermined[8,11,12].Nostrain-specificgenesandnoplasmids werefoundamong sensustricto strainsaftersequencingof fivestrains[8,11].Incontrast,ahighdegreeofallelicdiversitywasdetected:globalcomparisonoffivestrainsrevealed atotalof20,082siteswithSNPsdetectedinatleastoneof theanalyzedstrainsand,withapproximately6,000sitesbetweenanygivenpair.Thehighdegreeofgenecontentconservationsuggeststhatphenotypicdifferencesobservedin A.marginale mustbeduetosmallpolymorphismsbetweenstrainsratherthanwholegeneinsertionsordeletions.Therefore,weexploitedtheinterstraindiversityof A. marginale tomapthegeneticbasisunderlyingphenotypic differencesamongstrains. A.marginale genomesequencesareavailableforstrains thatclearlydifferinameasurablephenotype:transmission bythearthropodvector.TheSaintMaries,PuertoRico, Virginia,EM,6DEandSouthIdahostrainsareexamples ofefficientlytransmittedstrains[13-18].TheFloridastrain, hasbeenshowntohaveaverylowtransmissionefficiency asitwasnottransmittedusing>10timesthenumberof Dermacentorandersoni ticksroutinelyusedfortransmissionwiththeSt.Mariesstrain[17,19,20]. Duetothecompletegenecontentconservation,differencesintransmissionefficiencyin A.marginale arelikely tobeascribedtosequencevariationproducingvariantproteinsoraffectinggenetranscription.Indeed,precedenceis seeninbacterialpathogens,whereSNPshavebeendiscoveredthatprovideaselectiveadvantageinhostcolonization [21].Wecombinedtwogenomicsequencingapproaches inordertofindSNPsandtranscriptionalchangesthatsegregatewithtransmissionphenotype.Wefirstcompared thegenomesequencesoftwostrains,St.Mariesand Florida,whichdisplaycontrastingphenotypeswithrespect tothetraitofinterest,ticktransmissibility.CandidateSNPs includedpolymorphismsencodingnon-synonymoussubstitutionswithingenes,aswellasSNPslocatedwithinputativepromoterregions.EachSNPontheresultinglistwas evaluatedthroughcomparativegenomicsinthreeefficientlytransmissiblestrainsforitsconsistentsegregation withphenotype.Theremainingdifferencesweresequenced intwoadditionalefficientlytransmissiblestrains.Only SNPsthatwereuniquetothepoorlytransmissibleFlorida strainwhencomparedtosixefficientlytransmittedstrains wereretainedascandidates.Thisresultedinalistofcandidategenes,consistingofthosecontainingcandidateSNPs orlocateddownstreamofputativepromoterSNPs.TranscriptionalanalysisofcandidategenesbyRT-PCRrevealed genesthatweredifferentiallytranscribedinstrainswith distinctlydifferenttransmissionefficiencies.Tofindadditionaltranscriptionalchangesrelatedtothephenotype ofinterest,weperformedagenomewidetranscriptome comparisonusingRNA-seqtechnology.TotalmRNApopulationsfromtwo A.marginale strainswithdifferent transmissioncapabilitiesweresequencedusingthreedifferentplatforms.Thisstudymakesuseoftwosequencing approachesandfourdifferenttechnologiestoidentify genesinvolvedinarelevantmicrobialtrait.Wepresent,to ourknowledge,thedeepestana lysisofanobligateintracellularbacterialtranscriptomeduringthepathogen ’ snatural courseofinfection.ResultsComparativegenomicsidentifiesSNPsthatsegregate withtransmissionstatusComparisonofthepoorlytransmissibleFloridastrainwith theefficientlytransmittedSt.Mariesstrainproduceda totalof9,609SNPsevenlydistributedthroughoutthegenome(Figure1,Figure2,andAdditionalfile1).Twotypes ofSNPswerefurthercharacterized:thosethatresultedin non-synonymousaminoacidchangeswithingenesand SNPslocatedinputativepromoterregions.Forthepurposesofthisstudy,putativepromotersweredefinedas intergenicregionsimmediately50totranslationstartsites. GlobalcomparisonoftheseSNPswithgenomesequences ofthreeefficientlytransmittedstrains,PuertoRico,Virginia andSouthIdahoyielded241NSchangeswithingenes,and 62SNPsdistributedin27putativepromoters.Thesegenes andpromoterswerethenfurtheranalyzedintwoadditionalefficientlytransmittedstrains,6DEandEM,by performingtargetedsequencingoftheregionsofinterest. Thefinalcandidatelistincluded18genesthatcontainedat leastoneSNPencodinganon-synonymoussubstitution thatsegregatedwithtransmissionstatus,and14SNPs withinnineintergenicregionsthatcouldpotentiallyaffect thetranscriptionof11genes(Figure1,Additionalfile1).Pierl etal.BMCGenomics 2012, 13 :669 Page2of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 3

Altogether,comparativegenomicsidentified29candidate genes.TranscriptionanalysisofcandidategenesThese29geneswithSNPsintheircodingregionsorin theirputativepromoterregionswereanalyzedfortranscriptionalactivitybyusingRT-PCR,whichrevealedthatthe29 candidategenesweretranscribedinboththeefficiently transmittedSt.MariesandthepoorlytransmissibleFlorida strain(Additionalfile2).For the11genesflankingcandidatepromoterregions,therelativeexpressionratiowas analyzedfromtwosepar ateinfectionsusing msp5 asa steadystatecalibrator[22,23] .Thefoldchangesweretested forstatisticalsignificance bythepairwiserandomization testintwoseparateinfections. Statisticalsignificanceofthe averagefoldchangesacrossbothbiologicalreplicateswas testedusinganadaptationofthemethodproposedby Willemsetal.[24](Figure3A).Fourgenesweredifferentiallyexpressedintwobiologicalreplicates:AMF_553 showed4.3timesincreasedexpressionintheefficiently transmittedstrain(P<0.05).AMF_474,AMF_505and AMF_142showeddecreasedexpressioninthehighlytransmissiblestrainbyratiosof0.2,0.6and0.7respectively(P< 0.05).Wecalculatedanexpressioncutoffbyadding2 standarddeviationstotheaveragefoldchangeseeninall thestudiedgenes.Ofthesedifferentiallyexpressedcandidates,onlygenesAMF_474andAMF_553werebelowand abovethecalculatedcutoff,respectively.RNA-seqThetranscriptomesoftheFloridaandSt.Mariesstrains of A.marginale weresequencedusingthreedifferent technologies:454,Illumina,andIonTorrent.Roche ’ s 454technologyprovidedthelongestreads,asexpected (Table1).Interestingly,thistechnologyalsoyieldedthe highestpercentageof A.marginale readsat37.1%.AlthoughtheIlluminaplatformhadthelowestpercentage of A.marginale reads(4.7%fortheFloridastrain),this wascompensatedbydepthandwassufficientforquantitativeanalysis.Theuseofdifferentplatformsallowedus toaddresssomeofthechallengesofworkingwithobligateintracellularpathogens;asthesemicrobesare dependentontheireukaryotichostcells,RNAsamples aresignificantlycontaminatedbyhosttranscripts,and RNAprepshavebeenshowntobebiased[25].Our resultswerecorroboratedwiththedifferentplatforms. Transcriptomeanalysisallowedustoidentifyputative transcriptionstartsites( TSS)forbothstrainsasabyproductofourstudy.SeventyputativeTSSswerefoundin theFloridastrainand109werefoundfortheSt.Maries strain(Additionalfiles3and4).ThemajorityoftheseTSSs arepresentinbothlists,thelargernumberofhighconfidenceTSSsfoundintheSt.Mariesstraincanbeattributed tothedeepercoverageobtainedforthisstrain.Mostofthe Figure1 SNPssegregatedwithtransmissionstatusthrough wholegenomecomparisonandtargetedsequencing.A Genomewidecomparisonofthenon-transmissibleFloridastrain (red)withtheefficientlytransmittedSt.Maries(green)strain produced9609SNPs.FromthislistwesubtractedSNPsthatencode forsynonymouschanges,leavingtwotypesofSNPsthatwere furthercharacterized:thosethatresultedinnon-synonymous(NS) aminoacidchangeswithinORFsandSNPslocatedinputative promoterregions.ComparisonoftheseSNPswithgenome sequencesofthreeticktransmissiblestrainswasthenperformed. SNPsthatconsistentlysegregatedwithphenotypewereretained. Theremainingdifferenceswerethentargetedsequencedintwo additionalefficientlytransmissiblestrains. B .Atotalof9609SNPs werefoundbetweenthetransmissibleSt.MariesandthenontransmissibleFloridastrain(SNPs).Thiscomparisonfound4498nonsynonymousSNPs(representedinblack),1630SNPsfoundwithin putativepromoterregions(shownindarkgrey)andsynonymous SNPs(showninlightgray).Wholegenomecomparisonwiththree transmissiblestrainsallowedremovalof4127non-synonymousSNPs and1568promoterSNPsfromfurtherconsideration.Finally, Targetedsequencinginadditionaltransmissiblestrainsof241nonsynonymousand62promoterSNPsallowedretentionof35NSand 14promoterSNPsascandidateSNPsinvolvedinticktransmission. Pierl etal.BMCGenomics 2012, 13 :669 Page3of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 4

putative50untranslatedregions(UTR)foundwerelonger than40bpinbothstrains(63.3%inSt.Mariesand48.6% inFlorida).Fewerputative50UTRswerefoundtobesmallerthan40bp,andtheminoritywerefoundwithinthe predictedopenreadingframe(ORF)(Table2).The50UTRsfoundwithinannotatedCDSsuggestthatthepredictionsforthesegeneswereinaccurateandanadjustment inannotationisrequired.Wealsoidentified70highconfidenceoperonstructuresthatinvolved292differentgenes (Additionalfile5).Finally,sixregionswithnopreviousannotationwerefoundtodisplayhightranscriptionalactivity (Table3).Thesixregionsshowedtranscriptionalactivityin bothstrains,withtwoofthesenewlyidentifiedtranscripts showingsignificantdifferentialtranscriptionbetweenthe strains.TheseregionsareshowninTable3andFigure4, andarefurtherdiscussedinthenextsection.Transcriptomecomparisonidentifiestranscriptional differencesbetweenstrainswithcontrastingtransmission phenotypesAfternormalization,thedistributionsoftheexpression valuesacrossreplicateswerecomparedbeforeevaluating changesintranscription(Additionalfile6).Comparing thetranscriptomesofthehighlytransmissibleSt.Maries strainwiththepoorlytransmissibleFloridastrain producedalistof14genesthataresignificantlydifferentiallytranscribedusingourcriteria(seeMethods)and acrossreplicates(Figure3B,Figure5,andAdditional file7).Genesthatwerefoundtohavealowertranscription levelinthepoorlytransmittedFloridastrainareofparticularinterestconsideringtheexaminedtrait.Significantfold changedifferencesthatwereconstantacrossreplicatesrangedfrom3.5to413.0(Figure5,Figure3B).Ofthe10 genesthathadsignificantlylow(orabsent)transcriptional activityintheFloridastrain(Table4),onlyoneisannotated withapredictedlocalization:geneAMF_878,codingfor outermembraneprotein4(OMP4).Theotherninegenes areannotatedashypotheticalproteins.Threeofthesehad nomappedreadsinanyofthedifferentsequencingtechnologies:AMF_431,AMF_432andAMF_433.Anadditionaltwogenes,AMF_429andAMF_430,whichappear tobearrangedinanoperonwithAMF_431-3(basedon readsmappedtotheSt.Mariesgenome),arealsosignificantlydifferentiallytranscribed(Figure3B).RNA-seqanalysisofthepromotercandidat esidentifiedbycomparative genomicsconfirmedtheRT-PCRresults(Figure3A). Examinationofcandidatescar ryingnon-synonymousSNPs foundsignificantdifferentia ltranscriptionoftwogenes; AMF_793andAMF_1026(showninTable4,alongwith differentiallytranscribedpromotercandidatesAMF_474 Figure2 LocationofcandidateSNPsontheFloridastraingenome. ThiscircularrepresentationoftheFloridagenomeshowsinlightblue annotatedCDSs;outercirclerepresentsCDSsontheforwardstrand,innercirclerepresentsthereversestrand,ingreythe9609SNPsfound betweentheSt.MariesandtheFloridastraingenomes.Theelementsinlightgreenaremiscellaneousfeaturesannotatedinthegenome.Inthe innermostcircle49candidateSNPsfoundthroughcomparativegenomicsareshown.Redbarsshowthepositionofcandidatenon-synonymous SNPswithinCDSs.DarkbluebarsshowcandidateSNPsfoundwithinputativepromoterregions. Pierl etal.BMCGenomics 2012, 13 :669 Page4of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 5

***** *** * *** 0.1 1 10 100 1000 0.1 1 10 A BFoldchangein theSt. MariesstrainFoldchangein theSt. Mariesstrain Figure3 RNA-seqandqPCRconfirmtrendsintranscriptionalchangesbetweenstrainsthatdifferintheirticktransmissionstatus.A FoldchangeinthetransmissibleSt.Mariesstrainrelativetothenon-transmissibleFloridastrainforallpromotercandidatesexpressedinlogsca le 10.LocustagsforallgenesaregivenontheXaxis.BluebarsshowtheresultsobtainedafterevaluatingtwobiologicalreplicateswithRT-PCR. RedbarsshowthefoldchangeobtainedusingRNA-seqanalysisforthepromotercandidatesacrosstwobiologicalreplicates.Theasterisk indicatesstatisticalsignificanceatp<0.05. B .FoldchangeinthetransmissibleSt.Mariesstrainrelativetothenon-transmissiblefloridastrain expressedinlog10.Thetop18differentiallytranscribedgenesidentifiedthroughRNA-seqacrosstworeplicatesandtwostatisticaltestsandtheir foldchangesareshown.RedbarsshowresultsobtainedwithRNA-seq,bluebarsshowvalidationthroughqPCR.Theasteriskindicatesstatistical significanceatp<0.01. Table1Readsmappedto A.marginale fromthreesequencingplatformsPlatformStrainTotalreads A.marginale reads%ofreadsmappedto A.marginale Averagematchedreadlength Roche454 STM1726,051269,73037.1381.7 FL21,018,447326,44032.0396.5 IonTorrent STM1,004,747295,62929.4111.0 FL2,043,607577,28428.2111.2 Illumina STM88,650,7134,604,9935.2100 FL81,507,9673,845,8534.71001STM :St.Mariesstrain.2FL :Floridastrain.Pierl etal.BMCGenomics 2012, 13 :669 Page5of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 6

andAMF_553).Thesegeneshavealowertranscription levelintheFloridastrainbyfoldchangesof3.5and1.5respectively(p<1E-10).Finally,twonewlyidentifiedregions weredifferentiallytranscribed.Theregionsbetweenbp 336042and336685andbp1084944and1085520intheSt. Mariesgenome(Table3)wereup -regulatedbyfoldchanges of23.6and6.9,respectively(p<1E-10).Theseregionsare shownindetailinFigure4.Inordertodetermineifthese newlyidentifiedtranscriptscouldindicatethepresenceof genes,wesearchedforORFsthatwouldoverlapthese regions.Onlytworegions:frombp393765to394740and 1084944to1085520containedORFsthatwouldspanthe uninterruptedtranscript.Com parisonsofreplicateswere performedinordertoaccountfor variationoftranscription valueswithinastrain;import antly,thegenesthatwereconsistentlydifferentiallytran scribedacrossreplicateswere foundtobehomogeneouslytranscribedwhenstrainreplicateswerecomparedtoeachother.RT-PCRandvalidationofRNA-seqresultsThe18genesthatwerethemostdifferentiallytranscribed acrossreplicateswereanalyzedbyusingRT-PCRtoconfirmtheRNA-seqresults.Foldchangeintranscriptionwas evaluatedandcomparedwithRNA-seqanalysis.Asshown inFigure3B,transcriptionalchangeswereconfirmedand statisticallysignificantinallbutoneoftheanalyzedgenes. GeneAMF_209,foundtobemorehighlytranscribedin theSt.Mariesstrainby124foldwasnotconsistentlyupregulatedacrossbothreplicatesbyRT-PCR(up-regulated by18.2and1.3foldinseparatereplicates)and,therefore, itsfoldchangewasnotstatisticallysignificant.Genecharacteristics/bioinformaticsTable4showsthe30genesthatwereselectedascandidates.GenesthatwerefoundtobedifferentiallytranscribedthroughRNA-seqandRT-PCRareshownontop ofTable4;geneswithcandidateSNPsanddifferentialtranscriptionareshowninthemiddleofTable4.Therestof thegenescontainnon-synonymousSNPsthatsegregated throughcomparativegenomi cs.Thelengthofthecandidategenesvaries,withAMF_530beingthelongestat 10,479bpandAMF_1037theshortestat240bp.Twelve ofthecandidategenesareannotatedashypotheticalproteins(Table4).GenesAMF_474,AMF_553,AMF_480, AMF_762,AMF_764,AMF_824,AMF_893andAMF_878 areorthologsofgeneswithknownfunctions.GenesdownstreamfrompromoterSNPsincludedonetranslationinhibitor(AMF_474)andonegeneinvolvedinenergy consumingprocesses,nuoJ(AMF_553).Genescontaining non-synonymoussubstituti onsincludedorthologsfor DNAgyrase(AMF_480),atRNAsynthase(AMF_762),an aspartatekinase(AMF_764),acarboxypeptidaseinvolved incellenvelopebiogenesis(AMF_824)andalipoprotein releasingprotein(AMF_893).Aroleintransmissionisnot immediatelyapparentforthesegenes,infact,itisnotsurprisingthatmorethanhalfofthecandidateswereofunknownfunctionduetothelackofinformationonthe determinantsofticktransmission.Asearchforrelated genesrevealedthat18ofthecandidategeneshadhomologsinthetick-transmissiblehumanpathogen Anaplasma phagocytophilum (Table4).Tengenes,AMF_051, AMF_433,AMF_432,AMF_431,AMF_430,AMF_429, AMF_547,AMF_613,AMF_762,AMF_893,AMF_798, andAMF_793alsohadhomologsinticktransmitted Ehrlichia species.OnlygenesAMF_197,AMF_264,AMF_269, AMF_480,AMF_703,AMF_824andAMF_893had homologsinthreeticktransmitted Rickettsia species.Additionally,hypotheticalcandidatesAMF_1037,AMF_879, AMF_401andAMF_530hadnohomologsintheGenBankdatabase.Thesefindingsprovidetwomutuallyexclusivescenarios:ifagenewithhomologsinthe aforementionedtick-transmittedorganismsisresponsible forthetraitofinterest,thissuggestsacommonmechanismwithinabacterialorderorfamily.Alternatively,agene Table2Percentageofputative50UTRsaccordingto lengthStrain50UTRs<40bp50UTRs 40bp50UTRwithin predictedCDS1St.Maries 25.6863.3011.02 Florida 27.1548.5724.28150UTRwithinpredictedCDS:inthiscolumnwelistcaseswheretranscript mappingshowsthatthe50UTRandTSSarefoundwithinthepreviously predictedandannotatedCDSindicatingthatthepreviousannotationwas incorrect. Table3PreviouslyunannotatedareasthatexhibitedhightranscriptionalactivityRegion1LengthIdentitythroughblastXGenebeforeGeneafterFoldchange 251313..251855 543hypotheticalproteinAmarV_01231[Anaplasmamarginalestr.Virginia]AM294pep1AM259thiD0.7 336042..336685 644n/aAM380AM38223.6* 393765..394740 976DNA-bindingproteinHU[AnaplasmaphagocytophilumHZ]AM434pdxJAM4351.1 459343..459783 441hypotheticalproteinAmarM_02282[Anaplasmamarginalestr.Mississippi]AM504tRNA-Asn-11.3 887245..887579 335hypotheticalproteinPseS9_19739[Pseudomonassp.S9]AM969bioBAM973purL1.5 1084944..1085520 577hypotheticalproteinAmarM_05569[Anaplasmamarginalestr.Mississippi]AM1214polAAM12166.9*1basepairpositionsspannedbythenewlyidentifiedregionsintheSt.Mariesgenome. *Statisticallysignificantfoldchangesareindicatedwithanasterisk.Pierl etal.BMCGenomics 2012, 13 :669 Page6of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 7

Figure4 Newlyidentifiedtranscriptionallyactiveregionsofthegenome. MappingofcDNAreadstothe A.marginale genomeallowedus todetectregionswithoutpreviousannotationthatexhibitedtranscriptionalactivity. A showstheregionoftheSt.Mariesgenomethatspans frombp336042to336685.ThreedifferentgeneidentificationalgorithmsdidnotdetectaCDSthatwouldspanthelengthofthetranscript.The toppanelshowsthesixreadingframescontainingforty-fivestopcodons,shownasblackbars.Thebottompanelshowssomeofthemapped cDNAreadsingreenandred(indicatingdirectionoftheread).Thegreyhistogramunderthereadsrepresentsdepth(readheight).Thistranscript wasup-regulatedintheSt.Mariesstrainbyafoldchangeof23.7atp<1E-10. B showstheregionoftheSt.Mariesgenomethatspansfrombp 1084944to1085520.ThenewlyidentifiedgeneisfoundbetweengenespolA(notshown)andAM1216.OneORFontheleadingstrandseems tospanthelengthofthistranscriptandisshownasPUTATIVE_CDSinthisfigure.Thisnewgenewasfoundtobeup-regulatedinthe transmissibleSt.Mariesstrainbyafoldchangeof6.9atp<1E-10. Pierl etal.BMCGenomics 2012, 13 :669 Page7of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 8

uniqueto A.marginale wouldfavoraspecies-specific scenario. Fourofthegenesencodedtransmembranedomainsand signalpeptidespredictedthr oughmultiplealgorithms: AMF_798,AMF_793,AMF_824andAMF_878.The resultsobtainedforAMF_878arenotsurprising,asitis annotatedasoutermembraneprotein4(OMP4).Twentythreegeneshadsignificantscoresfortransmembranedomainpredictionsbutdidnotco ntainsignalpeptides.Analysisofnon-synonymousSNPsusingtheSIFTalgorithm [26]predictedeleventobedeleterious;thesesubstitutions arereportedinAdditionalfile8.DiscussionPairingcomparativegenomicswithhighthroughput RNA-seqanalysisallowsforidentificationofsequence andtranscriptionaldifferencesonagenomewidescale. Inthepresentstudy,comparativegenomicsreduceda listofcandidateSNPsfrom9,609to49SNPsthatsegregatewithtransmissionstatus,including35thatencode non-synonymoussubstitutionswithin18genesand14 residingwithinnineputativepromotersthatcouldaffect transcriptionof11genes.Oftheputativepromoter SNPs,weretainedonlythosethataffectedthetranscriptionofadjacentgenes,leavingjust2SNPsaffectingtwo genes,reducingtheoveralllistto37candidateSNPs affecting20genes.Deepsequencingandcomparative expressionanalysisfoundanadditional10geneswhose transcriptionbetweenstrainswithdistincttransmission efficienciesissignificantlydifferent.Transcriptomeanalysisalsorevealedtwopreviouslyun-annotatedregions thatweredifferentiallytranscribedbetweenthestrains ofinterest.Thisproducedafinallistof30genesand twonewlyidentifiedtranscriptionallyactiveregionsthat segregatewithticktransmission. OurcombinedapproachallowedustomapSNPsthat segregateamong A.marginale strainswithdivergent transmissionefficiencies.Suchsubtledifferenceshave beenshowntohavedramaticeffectsonorganismbiology.Asinglenon-synonymousSNPintheenvelopeproteingeneE1oftheChikungunyavirusisdirectly responsibleforachangeinvectorspecificitythatcaused anepidemicintheReunionIslandin2004[27].One SNPintheFimHadhesiongenefromacommensal strainof E.coli modifiedthisstrain ’ saffinityformonomannosereceptors,correlatingdirectlywithincreased uroepitheliumaffinityandallowingdetrimentalbladder colonization[21].Similarly,aSNPwithinthepromoter ofthenitratereductasegeneclusternarGHIJwasshown toberesponsibleforthedifferentnitratereductase phenotypeshownbythealmostidentical Mycobacterium bovis and Mycobacteriumtuberculosis ,bacterialspecies withidenticalgenecontent[28]. Comparativegenomicsidentified20geneswithatleast oneSNPthatsegregatedwithtransmissionphenotype. Thelackofinformationonthemicrobialdeterminants ofticktransmissionisconsistentwiththeobservation thatthemajorityofthegenescontainingnonsynonymousSNPsareofunknownfunction.Candidate geneswithorthologsinotherbacterialspeciesdonot appeartohaveanobviousinvolvementinthephenotype ofinterest.Threegenes:AMF_798,AMF_793and AMF_824,werepredictedtohavebothsignalpeptides RPKM Genes Figure5 WholegenomecomparisonoftranscriptionalactivityintheSt.MariesandFloridastrains. TheRPKMvaluesfor955genesfound intheFloridastraingenomeof A.marginale. RPKMvaluesareshownontheYaxis.Featuresarearrangedfromlefttorightastheyappearinthe genomeontheXaxis.ThenormalizedRPKMvalueswereplottedforeachstrain.RPKMvaluesforthetransmissibleSt.Mariesstrainareshownin redintheupperpartofthegraph;numbersforthenon-transmissibleFloridastrainareshowninlightblueinthelowerpartofthegraph.RPKM valuesfortheFloridastrainareplottedontheoppositesideofthexaxisforeaseofcomparison;theydonotrepresentnegativevalues. RibosomalRNA(rRNA)genesweresubtractedfromthiscomparison. Pierl etal.BMCGenomics 2012, 13 :669 Page8of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 9

andtransmembranedomains.Thepresenceofsignal peptidesandtransmembranedomainsimpliesmembranelocalizationoftheproteins,andthus,theseproteinswouldbemorelikelytointeractwithvector moleculesandthereforeeffecttransmission.Outofthe 35non-synonymouscandidateSNPs,alittleundera thirdwerepredictedtobedeleterious(Additionalfile8). GeneAMF_1026carriesthehighestnumberofdeleterioussubstitutionswithatotalofthreenon-synonymous SNPs.Thisgenewasalsofoundtobeup-regulatedin Table4Candidategenesinvolvedintransmissionphenotypesegregatedbypolymorphismsanddifferential transcriptionCategoryFL1St.M1ProductSNPscandidate SNPs SNPlocationHomologs2BI3Differentiallytranscribed genes AMF_433AM579Hypotheticalprotein10GeneAP,ER,EChTM AMF_432AM579Hypotheticalprotein10GeneAP,ER,ECh,ECaAMF_431AM580Hypotheticalprotein30GeneAP,ER,ECh,ECaTM AMF_430AM576Hypotheticalprotein220GeneAP,ER,ECh,ECaTM/DS AMF_429AM574Hypotheticalprotein150GeneAP,ER,EChTM AMF_798AM1055Hypotheticalprotein90GeneAP,ER,ECh,ECaTM/SP AMF_879AM1165Hypotheticalprotein140Gene-TM AMF_878AM1164Outermembraneprotein420GeneAPTM/SP AMF_401AM540Hypotheticalprotein120Gene/Prom-TM AMF_258AM347Hypotheticalprotein190GeneAPTM Differentiallytranscribed genesw/SNPs genescarryingcandidate SNPs AMF_474AM632Ribosome-associated inhibitorA 11PromoterER,ECh,ECaAMF_553AM748NADHDehydrogenaseI chainJ 11PromoterAP,ER,ECh,ECaTM AMF_793AM1048Hypotheticalprotein772Gene/PromAP,ER,ECh,ECaTM/SP AMF_1026AM1352Hypotheticalprotein185Gene/PromAP,ER,EChTM/DS Genescarryingcandidate SNPs AMF_051AM071Hypotheticalprotein51GeneAP,ER,ECh,ECaTM AMF_197AM265Hypotheticalprotein191GeneAP,RBTM AMF_264AM354Hypotheticalprotein212GeneER,ECh,ECaTM/DS AMF_265AM356Hypotheticalprotein651GeneAPAMF_269AM368Hypotheticalprotein431GeneRBTM/DS AMF_480AM644DNAgyraseB141GeneAP,ER,ECh,ECa,RB, RC,RR TM AMF_530AM712Hypotheticalprotein1381Gene-TM AMF_518AM689Hypotheticalprotein201GeneAPAMF_547AM742Hypotheticalprotein11GeneAP,ER,ECh,ECaDS AMF_613AM823Hypotheticalprotein174GeneAP,ER,ECh,ECaAMF_703AM919Hypotheticalprotein11GeneAP,ER,ECh,ECa,RB, RC,RR TM AMF_762AM1001methionyl-tRNAsynthetase231GeneAP,ER,ECh,ECa,RB, RC,RR TM/DS AMF_764AM1005aspartokinase131GeneAP,ER,ECh,ECaTM/DS AMF_824AM1091D-Ala-D-Ala carboxypeptidase 336GeneAP,ER,ECh,ECa,RB, RC,RR TM/SP AMF_893AM1183lipoprotein-releasing transmembraneprotein 83GeneER,ECh,ECaTM/DS AMF_1037AM345Hypotheticalprotein41Gene--1FL :LocustagintheFloridastrain St.M :locustagintheSt.Mariesstrain.2AP : Anaplasmaphagocytophilum ER : Ehrlichiaruminantium ECh : Ehrlichiachaffeensis ECa : Ehrlichiacanis RB : Rickettsiabellii RC : Rickettsiaconorii RR : Rickettsia rickettsia.3BI :Bioinformatics TM :Transmembranedomain, SP :Signalpeptide, Prom :promoter, DS :Deleterioussubstitution.Pierl etal.BMCGenomics 2012, 13 :669 Page9of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 10

theefficientlytransmittedstrainthroughRT-PCR. Interestingly,italsohadaSNPinitspromoterregion. ThispromoterSNPdidnots egregatewiththerestof thetransmissiblestrainsandthereforewasnot retainedasacandidate.Polymorphismswereretained ascandidatesifsixefficientlytransmissiblestrainsconsistentlydivergedwiththenucleotidefoundinthe poorlytransmittedFloridastrain.CandidateSNPs includednon-synonymouschangesinORFsandSNPs foundinputativepromoterregions. TwogeneswithcandidateSNPsintheirputativepromoterregionswerefoundtobedifferentiallytranscribed.AMF_553,morehighlytranscribedintheSt. Mariesstrain,isannotatedasNADHdehydrogenaseI chainJ(nuoJ).ThisispartofthemembranearmofrespiratorycomplexI,aconservedprotonpumping NADH:ubiquinoneoxidoreductaseinbacteria[29].Anothercloselyassociatedgenefromthiscomplex,nuoL, hasbeenfoundtobeup-regulatedintherelatedorganism Rickettsiaconorii whiledealingwithosmoticstress [30],suggestingthatenhancementofNADHdehydrogenaseexpressioninavector-transmittedbacterium couldberelatedtoanadaptationstrategynecessaryto surviveinthechangingosmolarityofafeedingtick[31]. AMF_474,morehighlytranscribedintheFloridastrain, containsconserveddomainsforamodulationprotein, theribosomeassociatedinhibitorA(RaiA)alsoknown asproteinY(PY).Thisproteinisacold-shockinduced ribosomebindingproteinthatinhibitstranslation[32]. PYbindsexclusivelytothe30Ssubunitofthe70Sribosome,andpreventingtheformationofinitiationcomplexesbypreventingthebindingofmRNAandinitiator fMet-tRNAtotheribosome[33].Whentemperature levelsreturnto37C,initiationofproteinsynthesisovercomesthePYinhibitionastRNAcompetemoreeffectivelywithPYinelevatedtemperatures.Relatedbacterial specieswhicharealsotransmittedby D.andersoni ,such as Rickettsiarickettsii ,areknowntoenter “ dormant ” stageswithinticks[34].Subsequentreversionofthis state,inaprocesstermed “ reactivation ” ,isthoughttobe duetoanincreaseintemperaturewhenthearthropod feedsonthemammalianhost.Thereforetheseobservationssuggestaninterestingscenarioasthisgenewas up-regulatedinthelowtransmissionefficiencyFlorida strain.Thelowtransmissionphenotypecouldbedueto ahaltintranslationproducedbyanupregulationofPY duringcoldshockresponse. Thethreeaforementioneddi fferentiallytranscribed geneswereidentifiedthroughcomparativegenomics. AlthoughallthreecarriedSNPsintheirpromoter regions,onlytwowereretainedascandidates.This exposesalimitationoftheapproachthatwasusedin thisstudy:polymorphismst hatdonotsegregatewith allthehighlytransmissibles trainsmaystillcontribute tothephenotypeofinterest.Inordertoconfirmthe differencesintranscriptionrevealedthroughRT-PCR andtofindfurtherchangesintranscriptionalactivitythatthestrategymighthaveoverlooked;thetranscriptomeoftwostrainswithcontrastingtransmission phenotypeswerecompared .Genomewidecomparison oftranscriptionalactivityconfirmedourRT-PCR resultsandfoundanadditional10genesthatweresignificantlydifferentiallytranscribed.Ofthese10genes onlyonehadapredictedlocalization:AMF_878correspondstoOMP4,anoutermembraneproteinand memberofthepfam01617superfamily[11].Among theremaininggeneswithn ofunctionalannotation threegenesstoodoutastheyexhibitedacompletelack oftranscriptionalactivityinthepoorlytransmitted Floridastrain:AMF_431,AMF_432andAMF_433. Thesegenesappeartobearrangedinanoperonalong withAMF_429andAMF_430,accordingtothetiling ofreadsmappedintheSt.Mariesstrain.AMF_429 andAMF_430werealsosignificantlydifferentially transcribedbetweenthestrains.GenesAMF_429, AMF_431andAMF_433containhighscoringconserveddomainsfortailandhead/tailconnectorphage proteins,withthehighestsimilarityfoundtophage proteinsfrom Wolbachiaspp .,arelatedbacterialsymbiontofarthropods.Althoughthiscouldopeninterestingpossibilities,asphagesplayanimportantrolefor Wolbachiaspp .withinthearthropodhost[35],nomobileelementsorintactprophageshavebeenidentified in A.marginale [11]. Typically,pathogenicbacteriathatcyclebetween arthropodandmammalianhostsmodifytheirtranscriptionalprofilestoadapttothesedifferentenvironments [36].Oneofthemajordifficultiesinvolvedinexamining generegulationofobligateintracellularpathogensisthe lowamountsofbacterialRNA,whichisco-isolatedwith largeamountsofhostRNA.Inordertoovercomethe limitedamountofbacterialRNA,previoustranscriptomicstudiesinterrogatinggenesusedforobligate intracellularsurvivalwereconductedusingmimetic conditionsofinfectioninan invitro environment [37-39].Whilethesestudiesprovideinsightintoalimitednumberofgenesregulatedbyspecificcues,theyare notrepresentativeofnaturalinfection.Exposingthe relatedpathogen R.rickettsii todifferentenvironmentalconditionsthatmimicitstransitionfromarthropodto mammalianhostshowedasurprisinglyminimaltranscriptionalresponse,withlessthan10geneschanging morethan3-foldinexpressionlevel[37].Thiscouldindicatethatpathogensintheorder Rickettsiales donot regulategenesspecificallyforgrowthwithinmammalian ortickcellsbutcontainaconservedsetofgenesthat arerequiredforgrowthinbothenvironments.TheobligateintracellularhabitatofpathogensinthisOrdermayPierl etal.BMCGenomics 2012, 13 :669 Page10of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 11

offersuchastableenvironmentthatthenecessityfor generegulationismuchlessthanthatoffacultative intracellularpathogens.Ourstudysearchedfortranscriptionaldifferencesbetweenstrainswithcontrasting transmissionprofilesinthenaturalhostofourmodel organism. Theuseofdifferentsequencingplatformsinthisstudy wasinstrumentalinconfirmingsignificantandconsistentchangesintranscriptionalactivity.Ithasbeenshown thatdifferentRNApreparationandselectionprocedures indeepsequencingexperimentscanleadtomeasurable over-orunder-representationofparticularRNAs[25]. Thisstudyprovedthatutilizingdifferenttechnologies allowedforcontrolofsourcesofpotentialbiasinRNA sequencing:allthreeplatformsusedforourstudygave thesameresults.Makinguseofvariousplatformswas alsoinstrumentalinourgoalofdescribingthe A.marginale transcriptomewiththehighestpossibleaccuracy. Inbacteria,theoverwhelminglyhighnumbersofreads incombinationwithrelativelysmallgenomesizeshas ledtotheassumptionthatcompleteornearlycomplete transcriptomesarebeinganalyzed.However,selecting forprokaryoticsequencesinanoceanofeukaryotic RNAsmakesaccuraterepresentationofRNApopulationsdaunting.Fewattemptshavebeenmadeatdescribingthetranscriptomeofobligateintracellularpathogens throughRNA-seq;notably,todate,thishasbeendonefor Chlamydia species[38,39]andthetick-transmittedpathogen A.phagocytophilum [40].Thedeepestanalysisgenerated854,242readsthatmappedtothe1.23Mb Chlamydiapneumonia genome[39];wemappedup to2,990,921readsperreplicateto A.marginale ’ s1.2Mb genome.Toenrichforprokaryoticsequences,previous attemptsatcharacterizingobligateintracellularmicrobial transcriptomesuseddifferentialcentrifugationof invitro grownbacteriainordertoseparatethebacteriafromhost cells.Thisprocedureislikelytostressthebacteriaand skewtheirtranscriptionalprofile.EnrichmentforoursampleswasperformedbyselectivehybridizationonceRNA populationswerecollected.AlthoughMastronunzioetal. usedasimilarenrichmentprocedure;theyonlydetected 187,284reads,representing11%oftheCDSsinthe A. phagocytophilum genome[40].Inthisstudy,99%ofthe CDSsinthe A.marginale weredetectedthroughtranscriptionalanalysis. AnalyzingtranscriptionalprofileswithRNA-seqallows ustoevaluate “ snapshots ” intimeofbacterialtranscriptomes;therefore,itisessentialtogeneratedatafrommore thanonereplicatetoprovideabroadermorereliablepictureoftranscriptionalchanges.ThedepthandreproducibilityofthisRNA-seqdatasetallowedformappingofthe physicalstructureofthe A.marginale transcriptome;includingpreviouslyunreportedtranscriptionallyactive regionsand50UTRlength.Sixregionswithnoprevious annotationweredetectedinbothstrains;twoofthesewere differentiallytranscribed.Theroleofthesetranscriptsis uncertainasonlytwoofthesewerepredictedtocontain ORFs.Themajorityofthehighconfidence50UTRswere longerthan40bpinbothstrains.PreviousstudiesofTSSs haveshownthatonlyaverysmallportionof50UTRsare longerthan40bpinbacteria[41,42].As50UTRshave beeninvolvedinregulationprocessesinbacteria,further investigationoftheseelementsmightrevealtranslational andtranscriptionalroles[43].Additionally,mappingof transcriptionaldataallowedustodefine70putativeoperonstructuresthatinvolved292genes,showingthatat least30%ofthegenesarepolycistronic.AlthoughRNAseqallowsustostudypolycistronicmessagesonagenome widescale,thedepthofthistechniquecoupledwithtiling arrayshaveshownthattheconceptofsimpleoperonsis questionable.Differentialexpressionofconsecutivegenes withinoperonsandconditiondependentmodulationhighlightthecomplexityoftranscriptionalregulationinbacteria[44].ConclusionsThisstudytakesadvantageofthehighinterstraindiversityofthisintracellularbacteriumtosignificantlyreduce thenumberofcandidatedifferencesthatcouldbe involvedintheticktransmissionphenotype.Marrying nextgenerationsequencingapproachesallowedusto generatealistofgenesdifferingatthetranscriptional andsequencelevelsinstrainswithcontrastingtransmissionstatus.Transformationofthetransmissiondeficient alleleintoatransmissioncompetentstrainwillfacilitate functionalanalysisofthesegenesinordertodetermine theirroleintransmissionbythearthropodhost.Althoughthesuccessfultransformationof A.marginale hasbeenachieved[45,46],stabletargetedgenereplacementhasnotbeenaccomplishedandisanecessarynext stepfordeterminingtheroleofthesegenesintick transmission.Identificationofgenesinvolvedintick transmissioninourmodelwillprovideanimportant firststeptowardthedevelopmentofnovelcontrolstrategiesfortick-bornepathogens,suchastransmissionblockingvaccines.MethodsEthicsstatementAnimalexperimentswereapprovedbytheInstitutional AnimalCareandUseCommitteeatUniversityofIdaho, USA,inaccordancewithinstitutionalguidelinesbased ontheU.S.NationalInstitutesofHealth(NIH)Guide fortheCareandUseofLaboratoryAnimals.StrainsTheFlorida,St.Maries,Virginia,PuertoRico,South Idaho,EM and6DEstrainsusedinthisstudyhavebeenPierl etal.BMCGenomics 2012, 13 :669 Page11of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 12

describedindetailelsewhere[47-51].TheSt.Maries, Virginia,PuertoRico,SouthIdaho,EM and6DEstrains arereproduciblytransmittedbytheReynoldsCreekstock of D.andersoni [13,14,16,17,52,53].TheFloridastrainhas notbeensuccessfullytransmittedbyanytickspecies,includingtheReynoldsCreekstock[15,18,19].ComparativegenomicsTheaccessionnumbersforthestrainsusedare:St. Maries:CP000030.1,Florida:CP001079.1,Viriginia: ABOR00000000.1,PuertoRico:ABOQ00000000.1,South Idaho:AFMY00000000.1.MUMmerv3.1[54]wasusedto compareaspreviouslydescribed[8]tocomparetheFloridaandSt.Mariesstrains.SNPsencodingsynonymous substitutionswerenotfurtheranalyzed.TherunMapping programoftheNewblersuitev2.5.3(454LifeSciences) wasusedwithdefaultsettingstocompareallreadsfrom theVirginia,PuertoRico,andSouthIdahostrainstothe completedFloridaandSt.Mariesgenomes.Allremaining SNPsfromtheinitialcomparisonwerethenchecked againstthethreestrains;iftheFloridasequencewas matchedinanyofthehighlytransmissiblestrains,that SNPwasremovedasacandidate.Illuminasequencingof theSt.Mariesstrainwasusedtoevaluatethefrequencies oftheSNPsfoundbetweentheFloridaandSt.Maries strain.SNPsthatwerefoundat100%frequencieswere highlightedinAdditionalfile1.TargetedsequencingTheremainingSNPswereexaminedviatargetedsequencingoftheSouthIdaho,EM and6DEstrains [50,55].PrimersweredesignedbyaligningtheSNPcontainingregionfromtheFloridaandSt.Mariesstrains andselectingprimerstoflankthepolymorphism.The resultingampliconsweregeneratedfromgenomicDNA, clonedintopCR4-TOPO(Invitrogen)andsequencedin bothdirectionsusingBigDyev3.1chemistryonanABI 3130XL(AppliedBiosystems).Sequenceanalysiseliminatedcandidatesasdescribedabove.AllcandidateSNPs wereresequencedintheFloridastrain,toverifytheoriginalgenomicsequence.ComparativetranscriptionalanalysisTotalRNAwasisolatedfrom A.marginale -infectedblood usingTRIzol(Invitrogen),permanufacturerdirections. ExpressionwasmeasuredusingquantitativereversetranscriptionPCRusingtheSYBRGreenERRT-PCRKit (Invitrogen).Briefly,1ugofRNAwasprocessedwiththe SuperscriptIIIFirststrandkit(Invitrogen)toobtain cDNA.Copynumberswerecorrectedtomorecloselyreflecttranscriptlevelsbasedonreversetranscriptionefficiency[52].Thesteadystate,singlecopygene msp 5was usedtocalibratetheRT-PCR.Relativeexpressionratios werecalculatedbyamathematicalmodel,whichincludes efficiencycorrectionofindividualtranscriptsthroughthe RESTsoftware[56].ThissoftwareusesthePairWise FixedReallocationRandomizationTesttoassessthestatisticalsignificanceoftheRT-PCRresultswhencomparing therelativeexpressionofthepromotercandidatesinboth theFloridaandSt.Mariesstrains.Adifferentialexpression foldcutoffvalueof3.2wasestablishedbycalculatingthe meanoftheaverageratiosobservedforallgenesanalyzed inthisstudyplus2standarddeviations.Inordertoassess thestatisticalrelevanceofthefindingsacrosstwobiologicalreplicates,anadaptationofthestandardization methodproposedbyWillemsandcoworkerswasused [24];thisincludesthreebasicsteps:logtransformation, meancenteringandautoscalling.Afterstandardizingthe data,statisticalsignificanceofthefoldchangesobserved betweenthestrainsacrossbothexperimentswasdeterminedbycalculationof95%confidenceintervals.This procedurewasappliedtoeachcandidategeneandwas alsousedforverificationoftranscriptionaldifferences foundbyRNA-seq.RNA-seqTheaccessionnumberforthisRNA-seqstudyis: SRP014580.TwoHolsteincalvesnegativefor A.marginale byMSP5cELISA,C1322andC1323,wereinoculatedwith theFloridaandtheSt.Mariesstrains,respectively.InfectionlevelsweretrackedbyanalysisofGiemsa-stained bloodsmearstocalculatethepercentageofparasitized erythrocytes(PPE).Bloodsamplesweretakenatsimilar levelsofparasitemia(3.5and4%PPE).TotalRNAwas isolatedfrom A.marginale -infectedbloodusingTRIzol (Invitrogen)perthemanufacturer ’ sdirections.Eukaryotic sequenceswerenegativelyselectedthroughhybridization usingtheMICROBEnrichkit(Ambion).Forsamplesprocessedfor454andIonTorrenttechnologies,probesfor bacterialribosomalRNAsfromtheRibominuskit(Invitrogen)wereaddedduringthesubtractivehybridization procedure.ForsamplesprocessedforIllumina,theDuplexSpecificthermostablenuclease(DSN)normalization protocolwasapplied.DatawasprocessedusingCLCGenomicsWorkbench(CLCBio).Mappingparameterswere adjustedtomapamaximumnumberofreadstothereferencebacterialgenomes.Thedistributionoftheexpression valuesforallsampleswasanalyzedandcompared. Normalizationbyquantileswasappliedtoadjustthedistributionsforfurthercomparison.FoldchangeswithrespecttoRPKM(ReadsPerKilobaseperMillionmapped reads)valueswerecalculated[57].Twodifferenttests wereappliedtoevaluatethestatisticalsignificanceoffold changes:Kal ’ sandBaggerly ’ sstatisticaltestsonproportions[58,59].Comparisonsofreplicateswereperformed inordertoaccountforvariationwithinastrain.These comparisonsshowedverylittlevariation:amaximumof 2%ofgeneshadfoldchangesaboveorbelow1.AsPierl etal.BMCGenomics 2012, 13 :669 Page12of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 13

variationwithinstrainswasassessedweproceededto comparethedifferentiallytransmittedstrains.Inorderto establishtranscriptionfoldchangecutoffs,therelationship betweenthep-valuesofthestatisticaltestsappliedandthe magnitudeofthedifferenceinexpressionvaluesofthe sampleswasplottedandevaluated.Thiswasdoneinorder toarrangegenesalongdimensionsofbiologicalandstatisticalsignificance[60].Geneswhoselog2foldchangewas aboveandbelow2and-2,respectively,andwhose-log10 p-valuewasabove10inbothreplicatecomparisonsand underbothstatisticaltestswereselectedforfurtherevaluation(Additionalfile7). Areasofthegenomethatwerenotpreviouslyannotated andshowed>0.5coverage(averagesequencedatacoveragedepth)werereportedwhenreadswereunambiguously mappedtothe A.marginale genome[42]. TherelativeperformancesnotedinTable1forthedifferentsequencingtechnologiesshouldnotbedirectlycompared,asthisstudywasnotdesignedtocomparethese platforms.Ashasbeennoted[25],differentlibrarypreparationsandsequencingtechnologiesfavorrecoveryofdifferenttranscripts.Thegoalofusingmultipletechnologieswas toverifythatunder-orover-representedtranscriptsinany strainwerenotbeingfavoredbythetechnologyused.PutativestartsiteidentificationPutativetranscriptstartsiteswereidentifiedusingthe rulesproposedbyPassalacquaetal.[42]:briefly;genes withcontinuouscoverageextendingintoacodirectional upstreamgenewereidentifiedasmembersofanoperon. Ifthesignal “ droppedoff ” intheintergenicsequenceupstreamoftheopenreadingframe,wedesignatedthepoint atwhichcoveragedroppedto0astheputativetranscriptionalstartsite.Coveragedepthwascalculatedforevery positionofeachgenome,andallgenesconsideredhadan averagecoveragescore>0.5abovethecalculatedaverage coveragesignal.PutativeTSSsthatwerefoundwiththe highestconfidence(i.e.TSSspresentinallreplicates)were groupedintwodifferenttablesaccordingtothelengthof the50UTRs,lessormorethan40bp.BioinformaticanalysisofcandidatesInordertorankthecandidates,twodifferentcriteriawere established.Thefirst,termed “ biologicalplausibilityofassociation ” ,examinestheannotationofthecurrentlyavailablegenomesandthepredictedfunctionofthecandidate gene,usingexistingknowledgeaboutbiologyandthe studiedphenotype[61].Inotherwords,isthecandidate genelikelytobeinvolvedintheexaminedphenotype accordingtoitsknownorpredictedfunction?Thesecond criterioninvolvestheuseofthree insilico analyses.The presenceofsignalpeptidesinthecandidategeneswas assessedbyusingSignalP4.0[62].Transmembrane domainswerepredictedusingtwodistinctalgorithms: TMpredandDenseAlignmentSurface(DAS)methods [63];onlygeneswithtransmembranedomainspredicted bybothalgorithmswerereported.The “ SortingTolerant FromIntolerant ” (SIFT)algorithm[26]usesasequence homology-basedapproachtoclassifyaminoacidsubstitutions,andwasusedtopredictifsubstitutionsinthecandidateallelesdetrimentalortoleratedbytheprotein.The searchforORFsinnewlyidentifiedtranscriptionallyactive regionswasperformedusingthreedifferenttools:CLC GenomicsWorkbench(CLCBio),NIH ’ sORFfinder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html)andORF (http://bioinformatics.biol.rug.nl/websoftware/orf/orf_start.php).AdditionalfilesAdditionalfile1: NucleotidepolymorphismsbetweentheSt. MariesandFloridastrains. SNPsbetweentheSt.MariesandFlorida strainsarelistedheretogetherwiththenucleotidesreportedforallthe additionalreportednucleotides. Additionalfile2: Absoluteexpressionvaluesofcandidategenes. Geneidentificationsareprovidedonthexaxisandthecopynumberper mlofbloodontheyaxis.Theblackbarsrepresentthenumbers obtainedfortheFloridastrainandthewhitebarsthenumbersfortheSt. Mariesstrain.Transcriptionofallcandidategenesisshowntogetherwith thecalibratorMSP5. Additionalfile3: MappingofputativeTSSand50UTRslengthin theSt.Mariesstrain. Thelocationandlengthof50UTRsintheSt. Mariesstrainarereported. Additionalfile4: MappingofputativeTSSand50UTRslengthin theFloridastrain. Thelocationandlengthof50UTRsintheFlorida strainarereported. Additionalfile5: Operonstrucutresfoundthroughtranscriptome sequencingintheSt.Mariesstrain. Genesinvolvedinthedifferent operonstructuresarereported. Additionalfile6: Distributionofthenormalizedexpressionvalues ofallreplicatesanalyzedinthisstudywithRNA-seq. Thedistribution ofthenormalizedRPKMvaluesforallreplicatesisplottedinaboxplot. RNA-SeqFL1andRNA-SeqFL2designatedistributionsforFloridastrain replicates1and2respectively.RNA-SeqSTM1andRNA-SeqSTM2 designateRPKMdistributionsforSt.Mariesreplicates1and2respectively. Thedistributionsallowforcomparisons. Additionalfile7: A.marginale genesarrangedalongdimensionsof biologicalandstatisticalsignificance. Avolcanoplotshowsthe relationshipbetweenthep-valuesofastatisticaltestandthemagnitude ofthedifferenceinexpressionvalues.Ontheyaxisthenegativelog10 p-valuesareplotted.Onthex-axisthelog2valuesofthefoldchanges seeninwholetranscriptomecomparison.Theredlineshighlightthe cutoffsforgenesthatwereanalyzedfurther.Onlygenespopulatingthe upperrightandleftquadrantsoftheplotundertwodifferentstatistical tests(Kal ’ sandBaggerly ’ s)werechosen.Thisplotshowsresultsobtained forKal ’ stest. Additionalfile8: Candidatenon-synonymouschangespredictedto bedeleteriousbytheSIFTalgorithm. Non-synonymouschanges predictedtobedeleteriousbytheSIFTalgorithmfoundbetweentheSt. MariesandFloridastrainsarereported. Abbreviations CDS:CodingDNASequence;NS:Non-Synonymous;ORF:OpenReading Frame;RPKM:ReadsPerKilobaseperMillionmappedreads;SNP:SingleNucleotidePolymorphism;UTR:UnTrasnlatedRegion.Pierl etal.BMCGenomics 2012, 13 :669 Page13of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 14

Competinginterests Theauthorsdeclarethattheyhavenocompetinginterests. Authors ’ contributions SAP,MJD,GHP,KABconceivedtheexperiments;SAP,MJD,DDperformed theexperiments;SAP,MJD,DD,GHP,KABanalyzedthedata;SAP,MJD,GHP KABwroteandeditedthemanuscript.Allauthorsreadandapprovedthe finalmanuscript. Acknowledgments Theauthorswouldliketoacknowledgetheexperttechnicalassistanceof Ms.XiaoyaCheng.ThisworkwassupportedbyUSDACREESNRICGP200435600-14175and2005-35604-15440,NationalInstitutesofHealthGrant AI44005,andWellcomeTrustGR075800M.SAPwassupportedinpartby fellowshipsfromthePoncinTrustandCONACyT. Authordetails1PrograminGenomics,DepartmentofVeterinaryMicrobiologyand Pathology,PaulG.AllenSchoolforGlobalAnimalHealth,WashingtonState University,Pullman,WA99164-7040,USA.2DepartmentofInfectiousDiseases andPathobiology,CollegeofVeterinaryMedicine,UniversityofFlorida, Gainesville,FL32611-0880,USA.3EmergingPathogensInstitute,Universityof Florida,Gainesville,FL32611-0880,USA. Received:7August2012Accepted:16November2012 Published:26November2012 References1.BorkP,DandekarT,Diaz-LazcozY,EisenhaberF,HuynenM,YuanY: Predictingfunction:fromgenestogenomesandback. JMolBiol 1998, 283 (4):707 – 725. 2.JimK,ParmarK,SinghM,TavazoieS: Across-genomicapproachfor systematicmappingofphenotypictraitstogenes. GenomeRes 2004, 14 (1):109 – 115. 3.LevesqueM,ShashaD,KimW,SuretteMG,BenfeyPN: Trait-to-gene:a computationalmethodforpredictingthefunctionofuncharacterized genes. CurrBiol 2003, 13 (2):129 – 133. 4.MakarovaKS,WolfYI,KooninEV: Potentialgenomicdeterminantsof hyperthermophily. TrendsGenet 2003, 19 (4):172 – 176. 5.HuynenMA,BorkP: Measuringgenomeevolution. ProcNatlAcadSciUSA 1998, 95 (11):5849 – 5856. 6.PellegriniM,MarcotteEM,ThompsonMJ,EisenbergD,YeatesTO: Assigningproteinfunctionsbycomparativegenomeanalysis:protein phylogeneticprofiles. ProcNatlAcadSciUSA 1999, 96 (8):4285 – 4288. 7.TettelinH,MasignaniV,CieslewiczMJ,DonatiC,MediniD,WardNL, AngiuoliSV,CrabtreeJ,JonesAL,DurkinAS, etal : Genomeanalysisof multiplepathogenicisolatesofStreptococcusagalactiae:implications forthemicrobial “ pan-genome ” ProcNatlAcadSciUSA 2005, 102 (39):13950 – 13955. 8.DarkMJ,HerndonDR,KappmeyerLS,GonzalesMP,NordeenE,PalmerGH, KnowlesDPJr,BraytonKA: Conservationinthefaceofdiversity: multistrainanalysisofanintracellularbacterium. BMCGenomics 2009, 10: 16. 9.PearsonT,BuschJD,RavelJ,ReadTD,RhotonSD,U ’ RenJM,SimonsonTS, KachurSM,LeademRR,CardonML, etal : Phylogeneticdiscoverybiasin Bacillusanthracisusingsingle-nucleotidepolymorphismsfromwholegenomesequencing. ProcNatlAcadSciUSA 2004, 101 (37):13536 – 13541. 10.StratiloCW,LewisCT,BrydenL,MulveyMR,BaderD: Single-nucleotide repeatanalysisforsubtypingBacillusanthracisisolates. JClinMicrobiol 2006, 44 (3):777 –782. 11.BraytonKA,KappmeyerLS,HerndonDR,DarkMJ,TibbalsDL,PalmerGH, McGuireTC,KnowlesDPJr: CompletegenomesequencingofAnaplasma marginalerevealsthatthesurfaceisskewedtotwosuperfamiliesof outermembraneproteins. ProcNatlAcadSciUSA 2005, 102 (3):844 – 849. 12.DarkMJ,Al-KhederyB,BarbetAF: Multistraingenomeanalysisidentifies candidatevaccineantigensofAnaplasmamarginale. Vaccine 2011, 29 (31):4923 – 4932. 13.EriksIS,StillerD,PalmerGH: ImpactofpersistentAnaplasmamarginale rickettsemiaontickinfectionandtransmission. JClinMicrobiol 1993, 31 (8):2091 – 2096. 14.FutseJE,BraytonKA,DarkMJ,KnowlesDPJr,PalmerGH: Superinfectionas adriverofgenomicdiversificationinantigenicallyvariantpathogens. ProcNatlAcadSciUSA 2008, 105 (6):2123 – 2127. 15.FutseJE,UetiMW,KnowlesDPJr,PalmerGH: TransmissionofAnaplasma marginalebyBoophilusmicroplus:retentionofvectorcompetencein theabsenceofvector-pathogeninteraction. JClinMicrobiol 2003, 41 (8):3829 – 3834. 16.RurangirwaFR,StillerD,PalmerGH: Straindiversityinmajorsurface protein2expressionduringticktransmissionofAnaplasmamarginale. InfectImmun 2000, 68 (5):3023 – 3027. 17.ScolesGA,BroceAB,LysykTJ,PalmerGH: Relativeefficiencyofbiological transmissionofAnaplasmamarginale(Rickettsiales:Anaplasmataceae) byDermacentorandersoni(Acari:Ixodidae)comparedwithmechanical transmissionbyStomoxyscalcitrans(Diptera:Muscidae). JMedEntomol 2005, 42 (4):668 – 675. 18.ScolesGA,UetiMW,NohSM,KnowlesDP,PalmerGH: Conservationof transmissionphenotypeofAnaplasmamarginale(Rickettsiales: Anaplasmataceae)strainsamongDermacentorandRhipicephalusticks (Acari:Ixodidae). JMedEntomol 2007, 44 (3):484 – 491. 19.UetiMW,ReaganJOJr,KnowlesDPJr,ScolesGA,ShkapV,PalmerGH: Identificationofmidgutandsalivaryglandsasspecificanddistinct barrierstoefficienttick-bornetransmissionofAnaplasmamarginale. InfectImmun 2007, 75 (6):2959 – 2964. 20.WickwireKB,KocanKM,BarronSJ,EwingSA,SmithRD,HairJA: Infectivity ofthreeAnaplasmamarginaleisolatesforDermacentorandersoni. AmJVetRes 1987, 48 (1):96 – 99. 21.SokurenkoEV,HastyDL,DykhuizenDE: Pathoadaptivemutations:gene lossandvariationinbacterialpathogens. TrendsMicrobiol 1999, 7 (5):191– 195. 22.RamabuSS,UetiMW,BraytonKA,BaszlerTV,PalmerGH: Identification ofAnaplasmamarginaleproteinsspecificallyup-regulated duringcolonizationofthetickvector. InfectImmun 2010, 78 (7):3047 – 3052. 23.VisserES,McGuireTC,PalmerGH,DavisWC,ShkapV,PipanoE,KnowlesDP Jr: TheAnaplasmamarginalemsp5geneencodesa19-kilodalton proteinconservedinallrecognizedAnaplasmaspecies. InfectImmun 1992, 60 (12):5139 – 5144. 24.WillemsE,LeynsL,VandesompeleJ: Standardizationofreal-timePCR geneexpressiondatafromindependentbiologicalreplicates. Anal Biochem 2008, 379 (1):127 – 129. 25.RaabeCA,HoeCH,RandauG,BrosiusJ,TangTH,RozhdestvenskyTS: TherocksandshallowsofdeepRNAsequencing:Examples intheVibriocholeraeRNome. RNA(NewYork,NY) 2011, 17 (7):1357 – 1366. 26.KumarP,HenikoffS,NgPC: PredictingtheeffectsofcodingnonsynonymousvariantsonproteinfunctionusingtheSIFTalgorithm. Nat Protoc 2009, 4 (7):1073 – 1081. 27.TsetsarkinKA,VanlandinghamDL,McGeeCE,HiggsS: Asinglemutationin chikungunyavirusaffectsvectorspecificityandepidemicpotential. PLoS pathogens 2007, 3 (12):e201. 28.StermannM,SedlacekL,MaassS,BangeFC: Apromotermutation causesdifferentialnitratereductaseactivityofMycobacterium tuberculosisandMycobacteriumbovis. JBacteriol 2004, 186 (9):2856 – 2861. 29.KaoMC,DiBernardoS,Nakamaru-OgisoE,MiyoshiH,Matsuno-YagiA,Yagi T: CharacterizationofthemembranedomainsubunitNuoJ(ND6)ofthe NADH-quinoneoxidoreductasefromEscherichiacolibychromosomal DNAmanipulation. Biochemistry 2005, 44 (9):3562 – 3571. 30.RenestoP,RoveryC,SchrenzelJ,LeroyQ,HuygheA,LiW,LepidiH, FrancoisP,RaoultD: Rickettsiaconoriitranscriptionalresponsewithin inoculationeschar. PLoSOne 2008, 3 (11):e3681. 31.MunderlohUG,KurttiTJ: Cellularandmolecularinterrelationships betweenticksandprokaryotictick-bornepathogens. AnnuRevEntomol 1995, 40: 221 – 243. 32.AgafonovDE,KolbVA,SpirinAS: Ribosome-associatedproteinthat inhibitstranslationattheaminoacyl-tRNAbindingstage. EMBOReports 2001, 2(5):399 – 402. 33.Vila-SanjurjoA,SchuwirthBS,HauCW,CateJH: Structuralbasisforthe controloftranslationinitiationduringstress. NatStructMolBiol 2004, 11 (11):1054 – 1059. 34.McDadeJE,NewhouseVF: NaturalhistoryofRickettsiarickettsii. AnnuRev Microbiol 1986, 40: 287 – 309.Pierl etal.BMCGenomics 2012, 13 :669 Page14of15 http://www.biomedcentral.com/1471-2164/13/669

PAGE 15

35.BordensteinSR,MarshallML,FryAJ,KimU,WernegreenJJ: Thetripartite associationsbetweenbacteriophage,Wolbachia,andarthropods. PLoS pathogens 2006, 2 (5):e43. 36.ObonyoM,MunderlohUG,FingerleV,WilskeB,KurttiTJ: Borrelia burgdorferiintickcellculturemodulatesexpressionofoutersurface proteinsAandCinresponsetotemperature. JClinMicrobiol 1999, 37 (7):2137 – 2141. 37.EllisonDW,ClarkTR,SturdevantDE,VirtanevaK,HackstadtT: Limited transcriptionalresponsesofRickettsiarickettsiiexposedto environmentalstimuli. PLoSOne 2009, 4 (5):e5612. 38.AlbrechtM,SharmaCM,ReinhardtR,VogelJ,RudelT: DeepsequencingbaseddiscoveryoftheChlamydiatrachomatistranscriptome. Nucleic AcidsRes 2010, 38 (3):868 – 877. 39.AlbrechtM,SharmaCM,DittrichMT,MullerT,ReinhardtR,VogelJ,RudelT: ThetranscriptionallandscapeofChlamydiapneumoniae. GenomeBiol 2011, 12 (10):R98. 40.MastronunzioJE,KurscheidS,FikrigE: Post-genomicanalysesreveal developmentofinfectiousAnaplasmaphagocytophilumduring transmissionfromtickstomice. JBacteriol 2012, 194 (9):2238 – 2247. 41.McGrathPT,LeeH,ZhangL,IniestaAA,HottesAK,TanMH,HillsonNJ,Hu P,ShapiroL,McAdamsHH: High-throughputidentificationof transcriptionstartsites,conservedpromotermotifsandpredicted regulons. NatBiotechnol 2007, 25 (5):584 – 592. 42.PassalacquaKD,VaradarajanA,OndovBD,OkouDT,ZwickME,Bergman NH: Structureandcomplexityofabacterialtranscriptome. JBacteriol 2009, 191 (10):3203 – 3211. 43.NavilleM,GautheretD: Transcriptionattenuationinbacteria:themeand variations. BriefFunctGenomicProteomic 2009, 8 (6):482 – 492. 44.GuellM,YusE,Lluch-SenarM,SerranoL: Bacterialtranscriptomics:whatis beyondtheRNAhoriz-ome? NatureReviews 2011, 9 (9):658 – 669. 45.FelsheimRF,ChavezAS,PalmerGH,CrosbyL,BarbetAF,KurttiTJ, MunderlohUG: TransformationofAnaplasmamarginale. VetParasitol 2010, 167 (2 – 4):167 – 174. 46.NohSM,UetiMW,PalmerGH,MunderlohUG,FelsheimRF,BraytonKA: Stabilityandticktransmissionphenotypeofgfp-transformedAnaplasma marginalethroughacompleteinvivoinfectioncycle. ApplEnvironMicrobiol 2011, 77 (1):330 – 334. 47.EriksIS,StillerD,GoffWL,PantonM,ParishSM,McElwainTF,PalmerGH: MolecularandbiologicalcharacterizationofanewlyisolatedAnaplasma marginalestrain. JVetDiagnInvest 1994, 6 (4):435 – 441. 48.KuttlerKL,WinwardLD: Serologiccomparisonsof4Anaplasmaisolatesas measuredbythecomplement-fixationtest. VetMicrobiol 1984, 9 (2):181 – 186. 49.McGuireTC,PalmerGH,GoffWL,JohnsonMI,DavisWC: Commonand isolate-restrictedantigensofAnaplasmamarginaledetectedwith monoclonalantibodies. InfectImmun 1984, 45 (3):697 – 700. 50.PalmerGH,KnowlesDPJr,RodriguezJL,GnadDP,HollisLC,MarstonT, BraytonKA: Stochastictransmissionofmultiplegenotypicallydistinct Anaplasmamarginalestrainsinaherdwithhighprevalenceof Anaplasmainfection. JClinMicrobiol 2004, 42 (11):5381 – 5384. 51.RisticCCA: Methodsofimmunoprophylaxisagainstbovineanaplasmosis withemphasisonuseofattenuatedAnaplasmamarginalevaccine. AdvExpMedBiol 1977, 93: 151 – 188. 52.LohrCV,RurangirwaFR,McElwainTF,StillerD,PalmerGH: Specific expressionofAnaplasmamarginalemajorsurfaceprotein2salivary glandvariantsoccursinthemidgutandisanearlyeventduringtick transmission. InfectImmun 2002, 70 (1):114 – 120. 53.UetiMW,KnowlesDP,DavittCM,ScolesGA,BaszlerTV,PalmerGH: Quantitativedifferencesinsalivarypathogenloadduringtick transmissionunderliestrain-specificvariationintransmissionefficiency ofAnaplasmamarginale. InfectImmun 2009, 77 (1):70 – 75. 54.KurtzS,PhillippyA,DelcherAL,SmootM,ShumwayM,AntonescuC, SalzbergSL: Versatileandopensoftwareforcomparinglargegenomes. GenomeBiol 2004, 5 (2):R12. 55.RodriguezJL,PalmerGH,KnowlesDPJr,BraytonKA: Distinctlydifferent msp2pseudogenerepertoiresinAnaplasmamarginalestrainsthatare capableofsuperinfection. Gene 2005, 361: 127 – 132. 56.PfafflMW,HorganGW,DempfleL: Relativeexpressionsoftwaretool (REST)forgroup-wisecomparisonandstatisticalanalysisofrelative expressionresultsinreal-timePCR. NucleicAcidsRes 2002, 30 (9):e36. 57.MortazaviA,WilliamsBA,McCueK,SchaefferL,WoldB: Mappingand quantifyingmammaliantranscriptomesbyRNA-Seq. NatureMethods 2008,5 (7):621 – 628. 58.BaggerlyKA,DengL,MorrisJS,AldazCM: DifferentialexpressioninSAGE: accountingfornormalbetween-libraryvariation. Bioinformatics(Oxford, England) 2003, 19 (12):1477 – 1483. 59.KalAJ,vanZonneveldAJ,BenesV,vandenBergM,KoerkampMG, AlbermannK,StrackN,RuijterJM,RichterA,DujonB, etal : Dynamicsof geneexpressionrevealedbycomparisonofserialanalysisofgene expressiontranscriptprofilesfromyeastgrownontwodifferentcarbon sources. MolBiolCell 1999, 10 (6):1859 – 1872. 60.CuiX,ChurchillGA: StatisticaltestsfordifferentialexpressionincDNA microarrayexperiments. GenomeBiol 2003, 4 (4):210. 61.TaborHK,RischNJ,MyersRM: Candidate-geneapproachesforstudying complexgenetictraits:practicalconsiderations. NatRevGenet 2002, 3 (5):391 – 397. 62.EmanuelssonO,BrunakS,vonHeijneG,NielsenH: Locatingproteinsin thecellusingTargetP,SignalPandrelatedtools. NatProtoc 2007, 2 (4):953 – 971. 63.CserzoM,WallinE,SimonI,vonHeijneG,ElofssonA: Predictionof transmembranealpha-helicesinprokaryoticmembraneproteins:the densealignmentsurfacemethod. ProteinEng 1997, 10 (6):673 – 676.doi:10.1186/1471-2164-13-669 Citethisarticleas: Pierl etal. : Comparativegenomicsand transcriptomicsoftrait-geneassociation. BMCGenomics 2012 13 :669. Submit your next manuscript to BioMed Central and take full advantage of: € Convenient online submission € Thorough peer review € No space constraints or color “gure charges € Immediate publication on acceptance € Inclusion in PubMed, CAS, Scopus and Google Scholar € Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit Pierl etal.BMCGenomics 2012, 13 :669 Page15of15 http://www.biomedcentral.com/1471-2164/13/669



PAGE 1

Normalized expression values RNA SeqFL1 RNA SeqFL2 RNA Seq STM2 RNA Seq STM1



PAGE 1

log 10 (p values ) Log 2 fold change ( proportions )


!DOCTYPE art SYSTEM 'http:www.biomedcentral.comxmlarticle.dtd'
ui 1471-2164-13-669
ji 1471-2164
fm
dochead Research article
bibl
title
p Comparative genomics and transcriptomics of trait-gene association
aug
au id A1 ca yes snm Pierlémnm Aguilarfnm Sebastiáninsr iid I1 email saguilar@vetmed.wsu.edu
A2 Darkmi JMichaelI2 I3 darkmich@ufl.edu
A3 DahmenDaniddahmen@vetmed.wsu.edu
A4 PalmerHGuygpalmer@vetmed.wsu.edu
A5 BraytonAKellykbrayton@vetmed.wsu.edu
insg
ins Program in Genomics, Department of Veterinary Microbiology and Pathology, Paul G. Allen School for Global Animal Health, Washington State University, Pullman, WA, 99164-7040, USA
Department of Infectious Diseases and Pathobiology, College of Veterinary Medicine, University of Florida, Gainesville, FL, 32611-0880, USA
Emerging Pathogens Institute, University of Florida, Gainesville, FL, 32611-0880, USA
source BMC Genomics
section Comparative and evolutionary genomicsissn 1471-2164
pubdate 2012
volume 13
issue 1
fpage 669
url http://www.biomedcentral.com/1471-2164/13/669
xrefbib pubidlist pubid idtype doi 10.1186/1471-2164-13-669pmpid 23181781
history rec date day 7month 8year 2012acc 16112012pub 26112012
cpyrt 2012collab Pierlé et al.; licensee BioMed Central Ltd.note 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.
kwdg
kwd Bacteria
Rickettsia
SNP
RNA-seq
Anaplasma
abs
sec
st
Abstract
Background
The Order it 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.
bdy
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
abbrgrp
abbr bid B1 1
. This approach has been used to predict genomic determinants for well-known phenotypes, including hyperthermophily, flagellar motility and pili assembly
B2 2
B3 3
B4 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
B5 5
B6 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
B7 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
B8 8
B9 9
B10 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
B11 11
B12 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
B13 13
B14 14
B15 15
B16 16
B17 17
B18 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
B19 19
B20 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
B21 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.
Results
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
figr fid F1 1, Figure
F2 2, and Additional file
supplr sid S1 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 5sup ′ 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). Altogether, comparative genomics identified 29 candidate genes.
suppl
Additional file 1
text
b Nucleotide polymorphisms between the St. Maries and Florida strains. SNPs between the St. Maries and Florida strains are listed here together with the nucleotides reported for all the additional reported nucleotides.
file name 1471-2164-13-669-S1.xlsx
Click here for file
fig Figure 1caption SNPs segregated with transmission status through whole genome comparison and targeted sequencing
SNPs segregated with transmission status through whole genome comparison and targeted sequencing. A. 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 non-transmissible Florida strain (SNPs). This comparison found 4498 non-synonymous 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 non-synonymous and 62 promoter SNPs allowed retention of 35 NS and 14 promoter SNPs as candidate SNPs involved in tick transmission.
graphic 1471-2164-13-669-1
Figure 2Location of candidate SNPs on the Florida strain genome
Location of candidate SNPs on the Florida strain genome. This circular representation of the Florida genome shows in light blue annotated CDSs; outer circle represents CDSs on the forward strand, inner circle represents the reverse strand, in grey the 9609 SNPs found between the St. Maries and the Florida strain genomes. The elements in light green are miscellaneous features annotated in the genome. In the inner most circle 49 candidate SNPs found through comparative genomics are shown. Red bars show the position of candidate non-synonymous SNPs within CDSs. Dark blue bars show candidate SNPs found within putative promoter regions.
1471-2164-13-669-2
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
S2 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
B22 22
B23 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.
B24 24
(Figure
F3 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.
Additional file 2
Absolute expression values of candidate genes. Gene identifications are provided on the x axis and the copy number per ml of blood on the y axis. The black bars represent the numbers obtained for the Florida strain and the white bars the numbers for the St. Maries strain. Transcription of all candidate genes is shown together with the calibrator MSP5.
1471-2164-13-669-S2.pptx
Click here for file
Figure 3RNA-seq and qPCR confirm trends in transcriptional changes between strains that differ in their tick transmission status
RNA-seq and qPCR confirm trends in transcriptional changes between strains that differ in their tick transmission status. A. Fold change in the transmissible St. Maries strain relative to the non-transmissible Florida strain for all promoter candidates expressed in log scale 10. Locus tags for all genes are given on the X axis. Blue bars show the results obtained after evaluating two biological replicates with RT-PCR. Red bars show the fold change obtained using RNA-seq analysis for the promoter candidates across two biological replicates. The asterisk indicates statistical significance at p < 0.05. B. Fold change in the transmissible St. Maries strain relative to the non-transmissible florida strain expressed in logsub 10. The top 18 differentially transcribed genes identified through RNA-seq across two replicates and two statistical tests and their fold changes are shown. Red bars show results obtained with RNA-seq, blue bars show validation through qPCR. The asterisk indicates statistical significance at p < 0.01.
1471-2164-13-669-3
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
tblr tid T1 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
B25 25
. Our results were corroborated with the different platforms.
table
Table 1
Reads mapped to
A. marginale
from three sequencing platforms
tgroup align left cols 6
colspec colname c1 colnum 1 colwidth 1*
center c2 2
c3 3
c4 4
c5 5
c6
thead valign top
row rowsep
entry
Platform
Strain
Total reads
A. marginale
reads
% of reads mapped to
A. marginale
Average matched read length
tfoot
1
STM: St. Maries strain.
2
FL: Florida strain.
tbody
morerows
Roche 454
STM1
726,051
269,730
37.1
381.7
FL2
1,018,447
326,440
32.0
396.5
Ion Torrent
STM
1,004,747
295,629
29.4
111.0
FL
2,043,607
577,284
28.2
111.2
Illumina
STM
88,650,713
4,604,993
5.2
100
FL
81,507,967
3,845,853
4.7
100
Transcriptome analysis allowed us to identify putative transcription start sites (TSS) for both strains as a by-product of our study. Seventy putative TSSs were found in the Florida strain and 109 were found for the St. Maries strain (Additional files
S3 3 and
S4 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 putative 5′ 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′ UTRs were found to be smaller than 40 bp, and the minority were found within the predicted open reading frame (ORF) (Table
T2 2). The 5′ 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
S5 5). Finally, six regions with no previous annotation were found to display high transcriptional activity (Table
T3 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
F4 4, and are further discussed in the next section.
Additional file 3
Mapping of putative TSS and 5
′ UTRs length in the St. Maries strain. The location and length of 5′ UTRs in the St. Maries strain are reported.
1471-2164-13-669-S3.xlsx
Click here for file
Additional file 4
Mapping of putative TSS and 5
′ UTRs length in the Florida strain. The location and length of 5′ UTRs in the Florida strain are reported.
1471-2164-13-669-S4.xlsx
Click here for file
Additional file 5
Operon strucutres found through transcriptome sequencing in the St. Maries strain. Genes involved in the different operon structures are reported.
1471-2164-13-669-S5.xlsx
Click here for file
Table 2
Percentage of putative 5
′ UTRs according to length
Strain
5
′ UTRs < 40 bp
5
′ UTRs ≥ 40 bp
5
′ UTR within predicted CDS
1
15′ UTR within predicted CDS: in this column we list cases where transcript mapping shows that the 5′ UTR and TSS are found within the previously predicted and annotated CDS indicating that the previous annotation was incorrect.
St. Maries
25.68
63.30
11.02
Florida
27.15
48.57
24.28
Table 3
Previously unannotated areas that exhibited high transcriptional activity
Region
1
Length
Identity through blastX
Gene before
Gene after
Fold change
1base pair positions spanned by the newly identified regions in the St. Maries genome.
*Statistically significant fold changes are indicated with an asterisk.
251313.251855
543
hypothetical protein AmarV_01231 [Anaplasma marginale str. Virginia]
AM294 pep1
AM259 thiD
0.7
336042.336685
644
n/a
AM380
AM382
23.6*
393765.394740
976
DNA-binding protein HU [Anaplasma phagocytophilum HZ]
AM434 pdxJ
AM435
1.1
459343.459783
441
hypothetical protein AmarM_02282 [Anaplasma marginale str. Mississippi]
AM504
tRNA-Asn-1
1.3
887245.887579
335
hypothetical protein PseS9_19739 [Pseudomonas sp. S9]
AM969 bioB
AM973 purL
1.5
1084944.1085520
577
hypothetical protein AmarM_05569 [Anaplasma marginale str. Mississippi]
AM1214 polA
AM1216
6.9*
Figure 4Newly identified transcriptionally active regions of the genome
Newly identified transcriptionally active regions of the genome. Mapping of cDNA reads to the A. marginale genome allowed us to detect regions without previous annotation that exhibited transcriptional activity. A shows the region of the St. Maries genome that spans from bp 336042 to 336685. Three different gene identification algorithms did not detect a CDS that would span the length of the transcript. The top panel shows the six reading frames containing forty-five stop codons, shown as black bars. The bottom panel shows some of the mapped cDNA reads in green and red (indicating direction of the read). The grey histogram under the reads represents depth (read height). This transcript was up-regulated in the St. Maries strain by a fold change of 23.7 at p < 1E-10. B shows the region of the St. Maries genome that spans from bp 1084944 to 1085520. The newly identified gene is found between genes polA (not shown) and AM1216. One ORF on the leading strand seems to span the length of this transcript and is shown as PUTATIVE_CDS in this figure. This new gene was found to be up-regulated in the transmissible St. Maries strain by a fold change of 6.9 at p < 1E-10.
1471-2164-13-669-4
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
S6 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
F5 5, and Additional file
S7 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
T4 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.
Additional file 6
Distribution of the normalized expression values of all replicates analyzed in this study with RNA-seq. The distribution of the normalized RPKM values for all replicates is plotted in a box plot. RNA-SeqFL1 and RNA-SeqFL2 designate distributions for Florida strain replicates 1 and 2 respectively. RNA-Seq STM 1 and RNA-Seq STM 2 designate RPKM distributions for St. Maries replicates 1 and 2 respectively. The distributions allow for comparisons.
1471-2164-13-669-S6.pptx
Click here for file
Additional file 7
A. marginale
genes arranged along dimensions of biological and statistical significance. A volcano plot shows the relationship between the p-values of a statistical test and the magnitude of the difference in expression values. On the y axis the negative log10 p-values are plotted. On the x-axis the log 2 values of the fold changes seen in whole transcriptome comparison. The red lines highlight the cutoffs for genes that were analyzed further. Only genes populating the upper right and left quadrants of the plot under two different statistical tests (Kal’s and Baggerly’s) were chosen. This plot shows results obtained for Kal’s test.
1471-2164-13-669-S7.pptx
Click here for file
Figure 5Whole genome comparison of transcriptional activity in the St. Maries and Florida strains
Whole genome comparison of transcriptional activity in the St. Maries and Florida strains. The RPKM values for 955 genes found in the Florida strain genome of A. marginale. RPKM values are shown on the Y axis. Features are arranged from left to right as they appear in the genome on the X axis. The normalized RPKM values were plotted for each strain. RPKM values for the transmissible St. Maries strain are shown in red in the upper part of the graph; numbers for the non-transmissible Florida strain are shown in light blue in the lower part of the graph. RPKM values for the Florida strain are plotted on the opposite side of the x axis for ease of comparison; they do not represent negative values. Ribosomal RNA (rRNA) genes were subtracted from this comparison.
1471-2164-13-669-5
Table 4
Candidate genes involved in transmission phenotype segregated by polymorphisms and differential transcription
9
c7 7
c8 8
c9
Category
FL
1
St.M
1
Product
SNPs
candidate SNPs
SNP location
Homologs
2
BI
3
1
FL: Locus tag in the Florida strain St.M: locus tag in the St. Maries strain.
2
AP: Anaplasma phagocytophilum, ER: Ehrlichia ruminantium, ECh: Ehrlichia chaffeensis, ECa: Ehrlichia canis, RB: Rickettsia bellii, RC: Rickettsia conorii, RR: Rickettsia rickettsia.
3
BI: Bioinformatics TM: Transmembrane domain, SP: Signal peptide, Prom: promoter, DS: Deleterious substitution.
Differentially transcribed genes
AMF_433
AM579
Hypothetical protein
1
0
Gene
AP, ER, ECh
TM
AMF_432
AM579
Hypothetical protein
1
0
Gene
AP, ER, ECh, ECa
-
AMF_431
AM580
Hypothetical protein
3
0
Gene
AP, ER, ECh, ECa
TM
AMF_430
AM576
Hypothetical protein
22
0
Gene
AP, ER, ECh, ECa
TM/DS
AMF_429
AM574
Hypothetical protein
15
0
Gene
AP, ER, ECh
TM
AMF_798
AM1055
Hypothetical protein
9
0
Gene
AP, ER, ECh, ECa
TM/SP
AMF_879
AM1165
Hypothetical protein
14
0
Gene
-
TM
AMF_878
AM1164
Outer membrane protein 4
2
0
Gene
AP
TM/SP
AMF_401
AM540
Hypothetical protein
12
0
Gene/Prom
-
TM
AMF_258
AM347
Hypothetical protein
19
0
Gene
AP
TM
Differentially transcribed genes w/SNPs genescarrying candidate SNPs
AMF_474
AM632
Ribosome-associated inhibitor A
1
1
Promoter
ER, ECh, ECa
-
AMF_553
AM748
NADH Dehydrogenase I chain J
1
1
Promoter
AP, ER, ECh, ECa
TM
AMF_793
AM1048
Hypothetical protein
77
2
Gene/Prom
AP, ER, ECh, ECa
TM/SP
AMF_1026
AM1352
Hypothetical protein
18
5
Gene/Prom
AP, ER, ECh
TM/DS
15
Genes carrying candidate SNPs
AMF_051
AM071
Hypothetical protein
5
1
Gene
AP, ER, ECh, ECa
TM
AMF_197
AM265
Hypothetical protein
19
1
Gene
AP, RB
TM
AMF_264
AM354
Hypothetical protein
21
2
Gene
ER, ECh, ECa
TM/DS
AMF_265
AM356
Hypothetical protein
65
1
Gene
AP
-
AMF_269
AM368
Hypothetical protein
43
1
Gene
RB
TM/DS
AMF_480
AM644
DNA gyrase B
14
1
Gene
AP, ER, ECh, ECa, RB, RC, RR
TM
AMF_530
AM712
Hypothetical protein
138
1
Gene
-
TM
AMF_518
AM689
Hypothetical protein
20
1
Gene
AP
-
AMF_547
AM742
Hypothetical protein
1
1
Gene
AP, ER, ECh, ECa
DS
AMF_613
AM823
Hypothetical protein
17
4
Gene
AP, ER, ECh, ECa
-
AMF_703
AM919
Hypothetical protein
1
1
Gene
AP, ER, ECh, ECa, RB, RC, RR
TM
AMF_762
AM1001
methionyl-tRNA synthetase
23
1
Gene
AP, ER, ECh, ECa, RB, RC, RR
TM/DS
AMF_764
AM1005
aspartokinase
13
1
Gene
AP, ER, ECh, ECa
TM/DS
AMF_824
AM1091
D-Ala-D-Ala carboxypeptidase
33
6
Gene
AP, ER, ECh, ECa, RB, RC, RR
TM/SP
AMF_893
AM1183
lipoprotein-releasing transmembrane protein
8
3
Gene
ER, ECh, ECa
TM/DS
AMF_1037
AM345
Hypothetical protein
4
1
Gene
-
-
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 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.
Gene characteristics/bioinformatics
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 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
B26 26
predicted eleven to be deleterious; these substitutions are reported in Additional file
S8 8.
Additional file 8
Candidate non-synonymous changes predicted to be deleterious by the SIFT algorithm. Non-synonymous changes predicted to be deleterious by the SIFT algorithm found between the St. Maries and Florida strains are reported.
1471-2164-13-669-S8.docx
Click here for file
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
B27 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
B28 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 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 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
B29 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
B30 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
B31 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
B32 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
B33 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
B34 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
B35 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
B36 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
B37 37
B38 38
B39 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
B40 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′ 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
B41 41
B42 42
. As 5′ UTRs have been involved in regulation processes in bacteria, further investigation of these elements might reveal translational and transcriptional roles
B43 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 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
B44 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
B45 45
B46 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.
Methods
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.
Strains
The Florida, St. Maries, Virginia, Puerto Rico, South Idaho, EMΦ and 6DE strains used in this study have been described in detail elsewhere
B47 47
B48 48
B49 49
B50 50
B51 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
B52 52
B53 53
. The Florida strain has not been successfully transmitted by any tick species, including the Reynolds Creek stock
15
18
19
.
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
B54 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
B55 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
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
B56 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
B57 57
. Two different tests were applied to evaluate the statistical significance of fold changes: Kal’s and Baggerly’s statistical tests on proportions
B58 58
B59 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
B60 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′ 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
B61 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
B62 62
. Transmembrane domains were predicted using two distinct algorithms: TMpred and Dense Alignment Surface (DAS) methods
B63 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).
Abbreviations
CDS: Coding DNA Sequence; NS: Non-Synonymous; ORF: Open Reading Frame; RPKM: Reads Per Kilobase per Million mapped reads; SNP: Single-Nucleotide Polymorphism; UTR: UnTrasnlated Region.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
SAP, MJD, GHP, KAB conceived the experiments; SAP, MJD, DD performed the experiments; SAP, MJD, DD, GHP, KAB analyzed the data; SAP, MJD, GHP KAB wrote and edited the manuscript. All authors read and approved the final manuscript.
bm
ack
Acknowledgments
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.
refgrp Predicting function: from genes to genomes and backBorkPDandekarTDiaz-LazcozYEisenhaberFHuynenMYuanYJ Mol Biol19982834707lpage 72510.1006/jmbi.1998.2144link fulltext 9790834A cross-genomic approach for systematic mapping of phenotypic traits to genesJimKParmarKSinghMTavazoieSGenome Res2004141109115pmcid 31428714707173Trait-to-gene: a computational method for predicting the function of uncharacterized genesLevesqueMShashaDKimWSuretteMGBenfeyPNCurr Biol200313212913310.1016/S0960-9822(03)00009-512546786Potential genomic determinants of hyperthermophilyMakarovaKSWolfYIKooninEVTrends Genet200319417217610.1016/S0168-9525(03)00047-712683966Measuring genome evolutionHuynenMABorkPProc Natl Acad Sci U S A199895115849585610.1073/pnas.95.11.5849344869600883Assigning protein functions by comparative genome analysis: protein phylogenetic profilesPellegriniMMarcotteEMThompsonMJEisenbergDYeatesTOProc Natl Acad Sci U S A19999684285428810.1073/pnas.96.8.42851632410200254Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: implications for the microbial “pan-genome”TettelinHMasignaniVCieslewiczMJDonatiCMediniDWardNLAngiuoliSVCrabtreeJJonesALDurkinASetal Proc Natl Acad Sci U S A200510239139501395510.1073/pnas.0506758102121683416172379Conservation in the face of diversity: multistrain analysis of an intracellular bacteriumDarkMJHerndonDRKappmeyerLSGonzalesMPNordeenEPalmerGHKnowlesDPsuf JrBraytonKABMC Genomics2009101610.1186/1471-2164-10-16264900019134224Phylogenetic discovery bias in Bacillus anthracis using single-nucleotide polymorphisms from whole-genome sequencingPearsonTBuschJDRavelJReadTDRhotonSDU’RenJMSimonsonTSKachurSMLeademRRCardonMLProc Natl Acad Sci U S A200410137135361354110.1073/pnas.040384410151875815347815Single-nucleotide repeat analysis for subtyping Bacillus anthracis isolatesStratiloCWLewisCTBrydenLMulveyMRBaderDJ Clin Microbiol200644377778210.1128/JCM.44.3.777-782.2006139315116517854Complete genome sequencing of Anaplasma marginale reveals that the surface is skewed to two superfamilies of outer membrane proteinsBraytonKAKappmeyerLSHerndonDRDarkMJTibbalsDLPalmerGHMcGuireTCKnowlesDPJrProc Natl Acad Sci U S A2005102384484910.1073/pnas.040665610254551415618402Multistrain genome analysis identifies candidate vaccine antigens of Anaplasma marginaleDarkMJAl-KhederyBBarbetAFVaccine201129314923493210.1016/j.vaccine.2011.04.131313368521596083Impact of persistent Anaplasma marginale rickettsemia on tick infection and transmissionEriksISStillerDPalmerGHJ Clin Microbiol1993318209120962657028370734Superinfection as a driver of genomic diversification in antigenically variant pathogensFutseJEBraytonKADarkMJKnowlesDPJrPalmerGHProc Natl Acad Sci U S A200810562123212710.1073/pnas.0710333105253888818252822Transmission of Anaplasma marginale by Boophilus microplus: retention of vector competence in the absence of vector-pathogen interactionFutseJEUetiMWKnowlesDPJrPalmerGHJ Clin Microbiol20034183829383410.1128/JCM.41.8.3829-3834.200317981212904396Strain diversity in major surface protein 2 expression during tick transmission of Anaplasma marginaleRurangirwaFRStillerDPalmerGHInfect Immun20006853023302710.1128/IAI.68.5.3023-3027.20009752310769008Relative efficiency of biological transmission of Anaplasma marginale (Rickettsiales: Anaplasmataceae) by Dermacentor andersoni (Acari: Ixodidae) compared with mechanical transmission by Stomoxys calcitrans (Diptera: Muscidae)ScolesGABroceABLysykTJPalmerGHJ Med Entomol200542466867510.1603/0022-2585(2005)042[0668:REOBTO]2.0.CO;216119558Conservation of transmission phenotype of Anaplasma marginale (Rickettsiales: Anaplasmataceae) strains among Dermacentor and Rhipicephalus ticks (Acari: Ixodidae)ScolesGAUetiMWNohSMKnowlesDPPalmerGHJ Med Entomol200744348449110.1603/0022-2585(2007)44[484:COTPOA]2.0.CO;217547235Identification of midgut and salivary glands as specific and distinct barriers to efficient tick-borne transmission of Anaplasma marginaleUetiMWReaganJOJrKnowlesDPJrScolesGAShkapVPalmerGHInfect Immun20077562959296410.1128/IAI.00284-07193285417420231Infectivity of three Anaplasma marginale isolates for Dermacentor andersoniWickwireKBKocanKMBarronSJEwingSASmithRDHairJAAm J Vet Res198748196993826850Pathoadaptive mutations: gene loss and variation in bacterial pathogensSokurenkoEVHastyDLDykhuizenDETrends Microbiol19997519119510.1016/S0966-842X(99)01493-610354593Identification of Anaplasma marginale proteins specifically up-regulated during colonization of the tick vectorRamabuSSUetiMWBraytonKABaszlerTVPalmerGHInfect Immun20107873047305210.1128/IAI.00300-10289737620439479The Anaplasma marginale msp5 gene encodes a 19-kilodalton protein conserved in all recognized Anaplasma speciesVisserESMcGuireTCPalmerGHDavisWCShkapVPipanoEKnowlesDPJrInfect Immun19926012513951442582891280624Standardization of real-time PCR gene expression data from independent biological replicatesWillemsELeynsLVandesompeleJAnal Biochem2008379112712910.1016/j.ab.2008.04.03618485881The rocks and shallows of deep RNA sequencing: Examples in the Vibrio cholerae RNomeRaabeCAHoeCHRandauGBrosiusJTangTHRozhdestvenskyTSRNA (New York, NY20111771357136610.1261/rna.2682311Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithmKumarPHenikoffSNgPCNat Protoc2009471073108119561590A single mutation in chikungunya virus affects vector specificity and epidemic potentialTsetsarkinKAVanlandinghamDLMcGeeCEHiggsSPLoS pathogens2007312e20110.1371/journal.ppat.0030201213494918069894A promoter mutation causes differential nitrate reductase activity of Mycobacterium tuberculosis and Mycobacterium bovisStermannMSedlacekLMaassSBangeFCJ Bacteriol200418692856286110.1128/JB.186.9.2856-2861.200438778215090527Characterization of the membrane domain subunit NuoJ (ND6) of the NADH-quinone oxidoreductase from Escherichia coli by chromosomal DNA manipulationKaoMCDi BernardoSNakamaru-OgisoEMiyoshiHMatsuno-YagiAYagiTBiochemistry20054493562357110.1021/bi047647715736965Rickettsia conorii transcriptional response within inoculation escharRenestoPRoveryCSchrenzelJLeroyQHuygheALiWLepidiHFrancoisPRaoultDPLoS One2008311e368110.1371/journal.pone.0003681257701018997861Cellular and molecular interrelationships between ticks and prokaryotic tick-borne pathogensMunderlohUGKurttiTJAnnu Rev Entomol19954022124310.1146/annurev.en.40.010195.0012537810987Ribosome-associated protein that inhibits translation at the aminoacyl-tRNA binding stageAgafonovDEKolbVASpirinASEMBO Reports200125399402108388511375931Structural basis for the control of translation initiation during stressVila-SanjurjoASchuwirthBSHauCWCateJHNat Struct Mol Biol200411111054105910.1038/nsmb85015502846Natural history of Rickettsia rickettsiiMcDadeJENewhouseVFAnnu Rev Microbiol19864028730910.1146/annurev.mi.40.100186.0014433096192The tripartite associations between bacteriophage, Wolbachia, and arthropodsBordensteinSRMarshallMLFryAJKimUWernegreenJJPLoS pathogens200625e4310.1371/journal.ppat.0020043146301616710453Borrelia burgdorferi in tick cell culture modulates expression of outer surface proteins A and C in response to temperatureObonyoMMunderlohUGFingerleVWilskeBKurttiTJJ Clin Microbiol1999377213721418510110364575Limited transcriptional responses of Rickettsia rickettsii exposed to environmental stimuliEllisonDWClarkTRSturdevantDEVirtanevaKHackstadtTPLoS One200945e561210.1371/journal.pone.0005612268098819440298Deep sequencing-based discovery of the Chlamydia trachomatis transcriptomeAlbrechtMSharmaCMReinhardtRVogelJRudelTNucleic Acids Res201038386887710.1093/nar/gkp1032281745919923228The transcriptional landscape of Chlamydia pneumoniaeAlbrechtMSharmaCMDittrichMTMullerTReinhardtRVogelJRudelTGenome Biol20111210R9810.1186/gb-2011-12-10-r98333378021989159Post-genomic analyses reveal development of infectious Anaplasma phagocytophilum during transmission from ticks to miceMastronunzioJEKurscheidSFikrigEJ Bacteriol201219492238224710.1128/JB.06791-11334707422389475High-throughput identification of transcription start sites, conserved promoter motifs and predicted regulonsMcGrathPTLeeHZhangLIniestaAAHottesAKTanMHHillsonNJHuPShapiroLMcAdamsHHNat Biotechnol200725558459210.1038/nbt129417401361Structure and complexity of a bacterial transcriptomePassalacquaKDVaradarajanAOndovBDOkouDTZwickMEBergmanNHJ Bacteriol2009191103203321110.1128/JB.00122-09268716519304856Transcription attenuation in bacteria: theme and variationsNavilleMGautheretDBrief Funct Genomic Proteomic20098648249210.1093/bfgp/elp02519651704Bacterial transcriptomics: what is beyond the RNA horiz-ome?GuellMYusELluch-SenarMSerranoLNature Reviews20119965866910.1038/nrmicro262021836626Transformation of Anaplasma marginaleFelsheimRFChavezASPalmerGHCrosbyLBarbetAFKurttiTJMunderlohUGVet Parasitol20101672–4167174281778019837516Stability and tick transmission phenotype of gfp-transformed Anaplasma marginale through a complete in vivo infection cycleNohSMUetiMWPalmerGHMunderlohUGFelsheimRFBraytonKAAppl Environ Microbiol201177133033410.1128/AEM.02096-10301971121057014Molecular and biological characterization of a newly isolated Anaplasma marginale strainEriksISStillerDGoffWLPantonMParishSMMcElwainTFPalmerGHJ Vet Diagn Invest19946443544110.1177/1040638794006004067858023Serologic comparisons of 4 Anaplasma isolates as measured by the complement-fixation testKuttlerKLWinwardLDVet Microbiol19849218118610.1016/0378-1135(84)90033-66203210Common and isolate-restricted antigens of Anaplasma marginale detected with monoclonal antibodiesMcGuireTCPalmerGHGoffWLJohnsonMIDavisWCInfect Immun19844536977002633526205996Stochastic transmission of multiple genotypically distinct Anaplasma marginale strains in a herd with high prevalence of Anaplasma infectionPalmerGHKnowlesDPJrRodriguezJLGnadDPHollisLCMarstonTBraytonKAJ Clin Microbiol200442115381538410.1128/JCM.42.11.5381-5384.200452527215528749Methods of immunoprophylaxis against bovine anaplasmosis with emphasis on use of attenuated Anaplasma marginale vaccineRisticCCAAdv Exp Med Biol19779315118810.1007/978-1-4615-8855-9_10596296Specific expression of Anaplasma marginale major surface protein 2 salivary gland variants occurs in the midgut and is an early event during tick transmissionLohrCVRurangirwaFRMcElwainTFStillerDPalmerGHInfect Immun200270111412010.1128/IAI.70.1.114-120.200212763811748171Quantitative differences in salivary pathogen load during tick transmission underlie strain-specific variation in transmission efficiency of Anaplasma marginaleUetiMWKnowlesDPDavittCMScolesGABaszlerTVPalmerGHInfect Immun2009771707510.1128/IAI.01164-08261226518955472Versatile and open software for comparing large genomesKurtzSPhillippyADelcherALSmootMShumwayMAntonescuCSalzbergSLGenome Biol200452R1210.1186/gb-2004-5-2-r1239575014759262Distinctly different msp2 pseudogene repertoires in Anaplasma marginale strains that are capable of superinfectionRodriguezJLPalmerGHKnowlesDPJrBraytonKAGene200536112713216202540Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCRPfafflMWHorganGWDempfleLNucleic Acids Res2002309e3610.1093/nar/30.9.e3611385911972351Mapping and quantifying mammalian transcriptomes by RNA-SeqMortazaviAWilliamsBAMcCueKSchaefferLWoldBNature Methods20085762162810.1038/nmeth.122618516045Differential expression in SAGE: accounting for normal between-library variationBaggerlyKADengLMorrisJSAldazCMBioinformatics (Oxford, England)200319121477148310.1093/bioinformatics/btg173Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sourcesKalAJvan ZonneveldAJBenesVvan den BergMKoerkampMGAlbermannKStrackNRuijterJMRichterADujonBMol Biol Cell1999106185918722538310359602Statistical tests for differential expression in cDNA microarray experimentsCuiXChurchillGAGenome Biol20034421010.1186/gb-2003-4-4-21015457012702200Candidate-gene approaches for studying complex genetic traits: practical considerationsTaborHKRischNJMyersRMNat Rev Genet20023539139710.1038/nrg79611988764Locating proteins in the cell using TargetP, SignalP and related toolsEmanuelssonOBrunakSvon HeijneGNielsenHNat Protoc20072495397110.1038/nprot.2007.13117446895Prediction of transmembrane alpha-helices in prokaryotic membrane proteins: the dense alignment surface methodCserzoMWallinESimonIvon HeijneGElofssonAProtein Eng199710667367610.1093/protein/10.6.6739278280