A Systems Biology Approach to the Characterization of Stress Response in Dermacentor reticulatus Tick Unfed Larvae

Background Dermacentor reticulatus (Fabricius, 1794) is distributed in Europe and Asia where it infests and transmits disease-causing pathogens to humans, pets and other domestic and wild animals. However, despite its role as a vector of emerging or re-emerging diseases, very little information is available on the genome, transcriptome and proteome of D. reticulatus. Tick larvae are the first developmental stage to infest hosts, acquire infection and transmit pathogens that are transovarially transmitted and are exposed to extremely stressing conditions. In this study, we used a systems biology approach to get an insight into the mechanisms active in D. reticulatus unfed larvae, with special emphasis on stress response. Principal Findings The results support the use of paired end RNA sequencing and proteomics informed by transcriptomics (PIT) for the analysis of transcriptomics and proteomics data, particularly for organisms such as D. reticulatus with little sequence information available. The results showed that metabolic and cellular processes involved in protein synthesis were the most active in D. reticulatus unfed larvae, suggesting that ticks are very active during this life stage. The stress response was activated in D. reticulatus unfed larvae and a Rickettsia sp. similar to R. raoultii was identified in these ticks. Significance The activation of stress responses in D. reticulatus unfed larvae likely counteracts the negative effect of temperature and other stress conditions such as Rickettsia infection and favors tick adaptation to environmental conditions to increase tick survival. These results show mechanisms that have evolved in D. reticulatus ticks to survive under stress conditions and suggest that these mechanisms are conserved across hard tick species. Targeting some of these proteins by vaccination may increase tick susceptibility to natural stress conditions, which in turn reduce tick survival and reproduction, thus reducing tick populations and vector capacity for tick-borne pathogens.

Despite its role as a vector of emerging or re-emerging diseases, very little information is available on the genome, transcriptome and proteome of D. reticulatus (115 nucleotide sequences of which only 15 were not of rRNA and 9 protein sequences deposited in the GenBank on June 2013).
This research focused on tick larvae because this is the first developmental stage to infest hosts, acquire infection and transmit pathogens that are transovarially transmitted. Additionally, D. reticulatus larvae hatch at a temperature range of 20-34uC and can survive for 83.5 days at 5uC and 100% relative humidity [4]. However, under natural conditions, larvae are active within 16-20 days after hatching and survive about a month before feeding [5]. D. reticulatus larvae feed on small mammals and are active during the summer [6].
All these facts put tick unfed larvae under extremely stressing conditions. For example, under natural conditions only 5-15% D.
reticulatus larvae produced from a single clutch are activated [5]. In this study, we characterized the transcriptome and proteome of D. reticulatus unfed larvae to get an insight into the mechanisms active at this developmental stage, with special emphasis on stress response.

D. reticulatus Unigenes Identified after Trasncriptomics Analysis of Unfed Larvae
A total of 21,677,414 (,2.1 Gb) Illumina 101 bp paired-end reads (207 bp average insert size) were subjected to analysis. After read assembly, 18,946 transcripts were obtained and annotated (Table S1). Transcripts were clustered by encoded proteins. If two transcripts were annotated as the same protein, then these transcripts were clustered together in the same protein cluster. We considered each set of transcripts annotated by the same protein as a unigene to identify transcripts from the same locus/ gene. This approach identified a set of 3,808 unigenes with 1,2316286 (Ave6S.E) estimated counts per unigene (Table S1).
The analysis of Biological Process (BP) and Molecular Function (MF) gene ontology (GO) showed that the most represented BPs corresponded to unknown process (N = 2,163; 57%), metabolic process (N = 411; 11%) and cellular process (N = 378; 10%) ( Fig. 1A) while proteins with unknown function (N = 2,163; 57%), catalytic activity (N = 658; 17%) and binding activity (N = 628; 16%) were the most represented MFs (Fig. 1B). A closer analysis of the most expressed genes showed that translation and structural constituent of the ribosome were the most represented BP and MF in D. reticulatus unfed larvae, respectively (Figs. 2A and 2B). These genes encoded for 80S ribosomal proteins (Table 1). With the exception of yeast, which lacks L28e, eukaryotic cytoplasmic 80S ribosomes contain the same set of 80 core ribosomal proteins [7]. Thus, the transcripts identified in D. reticulatus larvae encoded for 72% (34/47) and 73% (24/33) of the large and small subunit 80S proteins, respectively (Table 1), representing a high coverage for ribosomal proteins. These results showed that metabolic and cellular processes involved in protein synthesis were the most active in D. reticulatus unfed larvae (Figs. 1A, 1B, 2A, 2B), suggesting that tick metabolism is highly active during this life stage.

D. reticulatus Proteins Identified after Proteomics Analysis of Unfed Larvae
Proteomics analysis was replicated using two different experimental approaches to increase the probability of identifying low abundant proteins such as those involved in stress response. In both approached, mass spectra were searched against Ixodida protein database. The first approach used protein concentration and resulted in the identification of 74 proteins while the second approach analyzed proteins separated by SDS-PAGE and resulted in 239 proteins identified (Table S2), suggesting that for nonquantitative analysis protein fractionation provides better resolution. Of 74 proteins identified with the first approach, 26 (35%) were identified by both methods.
A recently described technique named proteomics informed by transcriptomics (PIT) [8] was used against data generated by the first proteomics approach to validate this method in ticks. This approach uses a database created from transcriptomics data to search mass spectra and has been reported to increase the number of identified proteins [8]. PIT approach resulted in 104 proteins identified in unfed tick larvae (Table S2), representing a 40% increase with respect to the search against Ixodida protein database. The analysis of de novo sequences increased the number of identified proteins using both approaches for proteomics data analysis (Table S2). However, while de novo protein sequences represented 4% (N = 3) of the identified proteins searching against Ixodida protein database, the number of identified proteins increased in 47% (N = 49) using PIT (Table S2). These results support the use of PIT for the analysis of proteomics data, particularly for organisms such as D. reticulatus with little sequence information available.
After removing proteins with unknown BP and MF, transcriptomics and proteomics data correlated well with respect to the most represented BPs (Figs. 3A-3C) and MFs (Figs. 4A-4C). These results were similar for both proteomics approaches, showing a good correlation in the proteomics analysis and providing additional support for the identified mechanisms active in D. reticulatus unfed larvae.

Rickettsia sp. Identified in D. reticulatus Unfed Larvae
Although these ticks were obtained from a colony considered to be free of tick-borne rickettsiae, some reads matching Rickettsia spp. were identified in D. reticulatus unfed larvae resulting in 16 unigenes (Table S3). These transcripts were probably wrongly annotated as I. scapularis proteins in Uniprot when they are likely Rickettsia proteins. In these cases, the Uniref representative protein of the cluster to which belongs the I. scapularis protein is a Rickecttsia protein and the rest of the members of the Uniref90 cluster are also from Rickettsia. Proteomics analysis corroborated the presence of Rickettsia proteins in D. reticulatus unfed larvae with the identification of 14 proteins searching against Rickettsiae database (Table S3).
The Rickettsia sp. identified in unfed larvae could be a commensal bacterium that has been described in Dermacentor and other tick species, but not in D. reticulatus [9][10][11][12] or a pathogen [13]. The Rickettsia proteins identified in D. reticulatus unfed larvae are highly conserved among Rickettsia spp. and thus not suitable to characterize these organisms at the species level.
To gain further information on this Rickettsia sp., the PCR amplification and sequencing of gene markers that have been previously used for species classification was conducted [14][15][16]. The results showed .99% pairwise nucleotide sequence identity to Rickettsia sp. sequences, especially to R. raoultii ( Table 2). As previously shown [15], the in silico PstI and RsaI restriction analysis of ompA sequences was highly informative and corroborated that the Rickettsia sp. identified in this study is similar to R. raoultii. These results suggested, as in previous studies in tick cell culture [17], that the Rickettsia sp. identified in D. reticulatus unfed larvae is closely related to the tick-borne pathogen, R. raoultii. However, until this Rickettsia sp. is fully characterized, we cannot exclude the possibility of a symbiont closely related to R. raoultii. These results suggested that the pathogen could be an additional stress factor in D. reticulatus unfed larvae, which correlated with the activation of immune response in these ticks (Figs. 1A, 3A and 3C). Rickettsia sequences were deposited in the GenBank with accession numbers [GenBank: KF478838, KF478839].

Stress Response in D. reticulatus Unfed Larvae
The results showed that metabolic processes and translation in particular were highly represented at the transcriptional level by genes encoding 80S ribosomal proteins in D. reticulatus unfed larvae (Table 1). Stress regulates ribosomal protein expression in other organisms, but no information is available in ticks [18][19][20][21]. Furthermore, a growing body of evidence suggests that the ribosome serves as a hub for co-translational folding, chaperone interaction, degradation, and stress response [22]. These results suggested a connection between transcription of ribosomal protein genes and stress response in ticks that deserves further investigation.
Transcripts and proteins mapped to stress response BP in D. reticulatus unfed larvae were selected for further analysis. Transcriptomics results showed that heat shock, cold shock and other stress responses were active in unfed larvae, represented by 39 unigenes (1% of all identified unigenes) and 27,937 counts (Table 3). Of them, the most represented functions corresponded to heat shock response (Figs. 5A and 5B). In general, protein identification has a lower resolution when compared to transcriptomics, a limitation that is more evident when working with species such as D. reticulatus for which sequence information is very scarce in the databases [23]. The search of MS data against the Ixodida database resulted in 8 stress response proteins identified (Table 4). However, when a database of transcripts identified as encoding for stress response proteins was generated and used for targeted PIT analysis, the results showed that 16 new stress response proteins were identified (Table 4). Additionally, while only 1% of the unigenes corresponded to stress response proteins, over 7% of the identified proteins were involved in this BP, supporting that stress response is active in tick unfed larvae. Furthermore, in agreement with transcriptomics data, the most represented function corresponded to heat shock response (Figs. 5C and 5D).
Some transcripts mapped to stress response BP were selected for the characterization of mRNA levels in D. reticulatus tick unfed larvae and guts and salivary glands from adult ticks incubated at 4, 37 or 19uC by real-time RT-PCR (Figs. 6A-6E). The results showed that all selected genes encoding for stress response proteins were more expressed in unfed larvae than in adult tissues, thus reinforcing the significance of this BP in D. reticulatus tick unfed larvae (Fig. 6A). In adult ticks, some genes were differentially expressed in response to temperature changes in guts or salivary glands (Figs. 6B-6E). The differential expression of selected genes encoding for stress response proteins was more evident in female salivary glands than in female guts and male tissues (Fig. 6E), suggesting differences between female and male ticks and between tissues in stress response to temperature changes. Additionally, at least for the genes characterized in this experiment, differential expression was more pronounced at 4uC than at 37uC (Fig. 6E), suggesting that D. reticulatus ticks respond differently to different temperatures. The sequences of D. reticulatus genes encoding for stress response proteins were deposited in the GenBank with accession numbers [GenBank: SRR950367; Bioproject: PRJNA214849].
Ticks are very sensitive to temperature and their life cycle is dependent on a complex combination of climate variables for development and survival [24]. In particular, D. reticulatus tick unfed larvae are exposed to extremely stressing conditions that affect their survival and development [5]. The heat-shock and other stress responses are a conserved reaction of cells and organisms to different stress conditions such as extreme temperatures, toxicity and pathogen infection [25]. Crucial to cell survival is the sensitivity of proteins and enzymes to heat inactivation and denaturation. Therefore, adaptive mechanisms exist that protect cells from the proteotoxic effects of stress. The heat-shock proteins and other stress response proteins protect cells and organisms from damage, providing higher levels of tolerance to environmental stress. Recent studies demonstrated that the stress response is activated in ticks and cultured tick cells in response to Anaplasma spp. infection and heat shock [26,27]. These results showed that at high temperatures and during blood feeding, when hsp20, hsp70 and subolesin are overexpressed, Ixodes scapularis ticks are protected from stress and pathogen infection and have a higher questing speed. Herein, genes encoding for stress response proteins were also differentially expressed in D. reticulatus in response to cold or heat shock. These results suggested a connection between tick stress response, questing behavior and pathogen infection [26,27], which may be present also in D. reticulatus tick unfed larvae.
Experiments characterizing the mRNA and protein levels of genes identified in this study in D. reticulatus ticks exposed to blood feeding and pathogen infection would add additional support to the importance of these proteins in tick stress response.

Conclusions
The characterization of the transcriptome and proteome of D. reticulatus unfed larvae showed that stress response was active in this developmental stage. Although descriptive in its nature, these results showed that combination of transcriptomics and proteomics approaches provide strong support for the characterization of biologically relevant pathways in ticks. The activation of stress responses in D. reticulatus unfed larvae likely counteracts the negative effect of temperature and other stress conditions such as Rickettsia infection and favors tick adaptation to environmental conditions to increase tick survival. These results are relevant to understand how D. reticulatus ticks have evolved mechanisms to  survive under stress conditions and suggest that these mechanisms are conserved across hard tick species. Targeting some of these proteins by vaccination may increase tick susceptibility to natural stress conditions, which in turn reduce tick survival and reproduction, thus reducing tick populations and vector capacity for tick-borne pathogens [28].

Experimental Design and Rationale
In this research, we completed the analysis of the transcripts and proteins present in D. reticulatus unfed larvae, which are described in Tables S1 and S2. This information, which was not available for this species, was then used to characterize stress response by focusing on the relevant genes and proteins. Individual variability, which certainly exists in ticks as in other organisms, was considered by pooling a large number of larvae for transcriptomics (N = 500) and proteomics (N = 200) studies. As in previous studies [29][30][31], we did not use biological replicates for RNA-Seq but the algorithm used to quantitate transcriptomics data allows the use of non-replicated samples [32]. Proteomics analysis, although also used for a non-comparative study that does not require replicates [31], was replicated using a different experimental approach to increase the probability of identifying low abundant proteins such as those involved in stress response. The statistical significance of reads and peptide assignments is supported by the application of eXpress and SEQUEST (FDR,0.01) algorithms described bellow for the analysis of transcriptomics and proteomics data, respectively.  Basel, Switzerland). Sample was sonicated for 1 min in an ultrasonic cooled bath followed by 10 sec vortex. After 3 cycles of sonication-vortex, the homogenate was centrifuged at 206g for 5 min at room temperature to remove cellular debris. The supernatant was collected and protein concentration was determined using the BCA Protein Assay (Thermo Scientific, San Jose, CA, USA) using BSA as standard.

Transcriptomics Data Acquisition
The RNA purified from unfed tick larvae was used for library preparation using the TruSeq RNA sample preparation kit v.1 and the standard low throughput procedure (Illumina, San Diego, CA, USA). Briefly, 0.7 mg total RNA was used as starting material for library preparation. Messenger RNA was captured using poly-dT magnetic beads and purified polyA+ RNA was chemically fragmented and reverse-transcribed. Remaining RNA was enzymatically removed and the second strand generated following an end repair procedure and preparation of double-stranded cDNA for adaptor ligation. Adaptor oligonucleotides containing the signals for subsequent amplification and sequencing were ligated to both ends and the cDNA was washed using AMPure SPRIbased magnetic beads (Beckman Coulter, IZASA, Barcelona, Spain). Adaptors contained identifiers, which allow multiplexing in the sequencing run. An enrichment procedure based on PCR was then performed to ensure that all molecules in the library conserved the adapters at both ends. The number of PCR cycles was adjusted to 15. The final amplified library was checked again on a BioAnalyzer 2100 (Agilent, Santa Clara, CA, USA) and titrated by quantitative PCR using a reference standard to characterize molecules concentration in the library (12.44 nM). The library was denatured and seeded on the lane of the flowcell at a final concentration after re-naturalization of 10-14 pM. After binding, clusters were formed in the flowcell by local amplification using a Cluster Station apparatus (Illumina). Following sequencing primer annealing, flowcell was loaded into a GAIIx equipment (Illumina) to perform sequencing using the TruSeqH system (Illumina). The sample was run under a pair-end 26100 bp protocol for de novo sequencing. After sequencing and quality filtering, reads were split according to their different identifiers and fastq files were generated to proceed to quality analysis and de novo transcript assembly and gene expression analysis.

Bioinformatics for Transcriptomics Data
Sequence reads were trimmed at the error probability higher than 0.05 and assembled only when two members of the pair remained after filtering at trimming. Oases [33] was used for read assembly in the mode of single (not merged) assembly because results were better in this mode. A K value of 79 was chosen, which was very close to the total length of the read (,100 bp) to avoid misassemblies since the higher the overlapping required the more accurate the transcript is. Final assembly was explored in detail using Tablet (http://bioinf.scri.ac.uk/tablet/download. shtml) [34].
Functional annotations were inferred by similarity to Uniprot reference proteins using Blast E values ,10E-10. We selected a set of 34,095 reference proteins downloaded from Uniprot on March 7, 2013, including all proteins that were representative of Uniref90 clusters belonging to the taxonomic node Chelicerata, which are 8 levels above D. reticulatus taxon. In the Uniref90 clusters, each protein belongs to only one cluster with a 90% similarity to the representative protein for all members of the cluster. It provides a more homogeneous and uniform distance between reference proteins. Reference proteins were used for transcript clusterization to obtain a protein-centred analysis of gene expression that is more useful for functional analysis in a de novo transcriptome.
The eXpress algorithm was used for mapping reads to multiple targets to quantify gene expression levels [32]. The eXpress algorithm [32] for quantifying the abundances of the transcripts addresses multi-mapping based on an on-line expectationmaximization algorithm (online-EM) [35] that is used to estimate transcript abundances in multi-isoform genes and gene families, and that does not require a reference genome. The underlying model is based on previously described probabilistic models developed for RNA-seq and allows the use of parameters for fragment length distributions, errors in reads, and sequencespecific fragment bias [36]. The algorithm alternates between assigning fragments to targets with a probability according to abundance parameters (expectation step) and updating abundances to the maximum-likelihood solution on the basis of the expectation-step assignments (maximization step). At the beginning the abundances are set to a uniform initial value. Then, for the fragments that map to multiple sites, eXpress calculates probabilities for each site, considering previous estimates of targetsequence abundances. As fragments are processed, they are assigned increasing 'mass' to improve the estimation of abundance according to the assignment probability. Parameters for fragmentlength (L) distribution, sequence bias and sequence read errors are updated and used in the next round of assignment. While relative  abundance and count estimates are updated, uncertainties in assignment are propagated so that posterior count distributions can be estimated. The probabilistic model is described in detail in the online methods section in Roberts and Pachter [32]. eXpress was also used to analyze the read mapping results. The mapper tool used was Bowtie setting the mapping parameters following the eXpress recommendations. Bowtie is an ultrafast, memory-efficient short read aligner that indexes the reference with a Burrows-Wheeler index to have low memory requirements [37]. Bowtie indexes the reference genome using a scheme based on the Burrows-Wheeler transform (BWT) [38], that is a reversible permutation of the characters in a text developed for data compression and the Ferragina and Manzini (FM) index [39]. Bowtie adopts the exact-matching algorithm of Ferragina and Manzini for searching in the FM index but introduces a qualityaware backtracking algorithm that allows mismatches and 'double indexing', to avoid excessive backtracking.
The script used for mapping the reads to the transcripts with Bowtie and for the final quantification with eXpress (File S1) was performed using cloud computing (Amazon Web Services). The process took 100 minutes in an Amazon EC2 m2.4xlarge instance. This kind of instances has 8 virtual CPUs and 68.4 GiB of RAM. desalted onto OMIX Pipette tips C 18 (Agilent Technologies, Santa Clara, CA, USA), dried-down and stored at 220uC until mass spectrometry analysis.
The desalted protein digest was resuspended in 0.1% formic acid and analyzed by RP-LC-MS/MS using an Easy-nLC II system coupled to an ion trap LCQ Fleet mass spectrometer (Thermo Scientific). The peptides were concentrated (on-line) by reverse phase chromatography using a 0.1620 mm C18 RP precolumn (Thermo Scientific), and then separated using a 0.0756100 mm C18 RP column (Thermo Scientific) operating at 0.3 ml/min. Peptides were eluted using a 180-min gradient from 5 to 35% solvent B (Solvent A: 0,1% formic acid in water, solvent B: 0,1% formic acid in acetonitrile). ESI ionization was done using a Fused-silica PicoTip Emitter ID 10 mm (New Objective, Woburn, MA, USA) interface. Peptides were detected in survey scans from 400 to 1600 amu (1 mscan), followed by three data dependent MS/MS scans (Top 3), using an isolation width of 2 mass-to-charge ratio units, normalized collision energy of 35%, and dynamic exclusion applied during 30 sec periods.
Proteins separated by SDS-PAGE. The protein extract (150 mg) was precipitated following the methanol/chloroform procedure [40], resuspended in 100 ml Laemmli sample buffer and applied onto 1.2-cm wide wells on a 12% SDS-PAGE gel. The protein bands were visualized by staining with GelCode Blue Stain Reagent (Thermo Scientific) and sliced each gel lane into 25 slices as previously described [42]. Protein digestion and RP-LC-MS/MS analysis was performed as described before for proteins concentrated by SDS-PAGE.

Bioinformatics for Proteomics Data
The MS/MS raw files were searched against Ixodida (40,849 entries in June 2013) and Rickettsieae (58,899 entries in June 2013) Uniprot databases and against a database created from transcriptomics data (PIT) [8] using the SEQUEST algorithm (Proteome Discoverer 1.3, Thermo Scientific) with the following constraints: tryptic cleavage after Arg and Lys, up to two missed cleavage sites, and tolerances of 1 Da for precursor ions and 0.8 Da for MS/MS fragment ions and the searches were performed allowing optional Met oxidation and Cys carbamidomethylation. For peptide validation, the Percolator node present in the Proteome Discoverer 1.3 software was used. Percolator is a machine-learning supplement of the Sequest algorithm that uses a decoy database search strategy to learn to distinguish between correct and incorrect peptide identifications increasing the sensitivity and specificity of peptide identification [43,44]. The filtering criteria applied in this case are based on the q-value generated by Percolator that is defined as the minimal false discovery rate at which the identification is deemed correct [43]. These q-values are estimated using the distribution of scores from the decoy database search. A false discovery rate (FDR) ,0.01 was considered as condition for successful peptide assignments, including only peptides with q-values #0.01 and delta Cn . 0.05. De novo peptide sequencing was conducted with Peaks Studio 6.0 software (Bioinformatics Solutions Inc., Waterloo, ON Canada).

Gene and Protein Ontology Assignments
Functional data for each protein were obtained from Uniprot and included GO annotations, EC number and Interpro motifs. Assignment of GO terms to identified proteins was done by Blast2GO software (version 2.6.6; http://www.blast2go.org/) in three main steps: blasting to find homologous sequences, mapping to collect GO-terms associated to blast hits and annotation to assign functional terms to query sequences from the pool of GO terms collected in the mapping step [45]. Sequence data of identified proteins were uploaded as FASTA file to the Blast2GO software and the function assignment was based on GO database. The blast step was performed against NCBI public databases through blastp. Other parameters were kept at default values: evalue threshold of 1e-3, recovery of 20 hits per sequence, minimal alignment length (hsp filter) 33 (to avoid hits with matching region smaller than 100 nucleotides) and Blast mode was set to QBlast-NCBI. Configuration for annotation was an e-value-Hit-filter of 1.0E-6, annotation cut off of 55 and GO weight of 5. For visualizing the functional information (GO categories: Molecular Function and Biological process), the analysis tool of the Blast2GO software was used.
The GO analysis for the 500 more represented unigenes was based on the GO annotations included in the Uniprot entry of the representative protein of each cluster. The GO analysis was done using Bio4j Go Tools developed by Era7 Bioinformatics and available at http://gotools.bio4j.com:8080/Bio4jTestServer/ Bio4jGoToolsWeb.html. Bio4j Go Tools is a set of GO related Web Services using the open source graph bioinformatics platform Bio4j as back-end. Bio4j is a graph-based database including most data available in UniProt KB (SwissProt+Trembl), Gene Ontology (GO), UniRef (50,90,100), RefSeq, NCBI taxonomy, and Expasy Enzyme (http://bio4j.com/). Specifically designed java programs were used for the generation of the GO frequency chart data.

Effect of Temperature on Gene Expression
Three groups of 9 D. reticulatus unfed female or male adults each were incubated for 4.5 h at 4uC, 19uC or 37uC and 20% relative humidity. After incubation, ticks were dissected and salivary glands and guts were separated, pooled in groups of three (3 groups for each temperature and sex) and immediately stored in TriReagent (Sigma, St. Louis, MO, USA) for RNA extraction.

Analysis of mRNA Levels by Real-time RT-PCR
For real-time RT-PCR, tick larvae (three pools of 100 larvae each), unfed female and male adult ticks (3 ticks each) and guts and salivary glands from unfed female and adult adults incubated at different temperatures (3 groups for each temperature and sex) were used for RNA extraction using TriReagent (Sigma, St. Louis, MO, USA) following manufacturer's recommendations. Real-time RT-PCR was performed on tick RNA samples (5 ng) with gene specific primers (20 pmol each) and conditions (Table 5) using the iScript One-Step RT-PCR Kit with SYBR Green and the iQ5 thermal cycler (Bio-Rad, Hercules, CA, USA) following manufacturer's recommendations. A dissociation curve was run at the end of the reaction to ensure that only one amplicon was formed and that the amplicons denatured consistently in the same temperature range for every sample [46]. The mRNA levels were normalized against tick ribosomal protein S4 using the genNorm method (ddCT method as implemented by Bio-Rad iQ5 Standard Edition, Version 2.0) [47,48]. Normalized Ct vales were compared Table 5. Primer sequences used for real-time RT-PCR. between larvae and adult samples and between samples from ticks incubated at 4 or 37uC and 19uC by Student's t-test with unequal variance (P = 0.05).

PCR and Sequence Analysis of Rickettsia Amplicons
Rickettsia sp. DNA was characterized by PCR, cloning and sequence analysis of the amplicons. At least three clones were sequenced for each amplicon. Genes targeted by PCR included fragments of ATP synthase alpha subunit (atpA), heat-shock protein 70 (dnaK), outer membrane protein A (ompA), outer membrane protein B (ompB), 16S rRNA, and recA [14][15][16]. Nucleotide sequence identity to reference strains and in silico PstI and RsaI restriction analysis of ompA sequences was used to characterize Rickettsia sp. [15,16].

Supporting Information
Table S1 Tick transcripts and encoded proteins identified in transcriptomics analysis of D. reticulatus unfed larvae.
(XLSX)  File S1 Script used for mapping the reads to the transcripts with Bowtie and for the final quantification with eXpress. (DOCX)