Introduction

Vector-borne microbial pathogens are transmitted by the bite of arthropods and have evolved sophisticated ways to adapt to vastly different environments as they move between vector and host. Uncovering their adaptive mechanisms can not only open avenues for disrupting pathogen transmission but also provide fundamental insights into vector physiology and microbial symbioses (Shaw and Catteruccia, 2018). Lyme disease, the most reported vector-borne disease in North America, is caused by the bacterial pathogen Borrelia burgdorferi (Bb) (Rosenberg et al., 2018; Steere et al., 2016). Its primary vector, the blacklegged tick Ixodes scapularis, acquires and transmits Bb through two separate multi-day bloodmeals, one in which Bb is acquired from an infected vertebrate host and a second, during the subsequent life stage, in which Bb is transmitted to a new host (Tilly et al., 2008). During the transmission bloodmeal, Bb proliferate in the tick midgut before a subset of these cells disseminate to the salivary glands (Dunham-Ems et al., 2009). Bb is deposited into the new host via the tick saliva extruded into the bite site (Spielman et al., 1987). Given the prolonged nature of I. scapularis feeding, the high specificity of vector-pathogen relationships, and the complicated array of events needed for successful Bb transmission, the tick bloodmeal provides an opportune intervention point for preventing pathogen spread. However, we do not currently have a clear understanding of the molecular mechanisms involved in this process.

Bb must adapt to dramatically different environments as it cycles from tick to vertebrate host, and understanding the genes involved in this process will enable the identification of key interactions to target to prevent transmission. When an infected tick feeds, Bb responds to bloodmeal-induced environmental changes and undergoes cellular modifications driven by key transcriptional circuits, including the RpoN-RpoS sigma factor cascade and the Hk1/Rrp1 two-component system (Radolf et al., 2012). In vitro analyses of Bb cells cultured in tick- or mammalian-like growth conditions have pointed to additional genetic determinants of tick-borne transmission by revealing widespread transcriptome remodeling during host switching (reviewed in Samuels et al., 2021). However, it is not fully clear how these in vitro expression changes correspond to complex in vivo changes over the course of a transmission bloodmeal, which comprises several physiologically distinct stages for the tick vector. Capturing comprehensive, longitudinal data on Bb gene expression from inside its vector has been hampered by technical challenges given the dynamic nature of the bloodmeal and the general low abundance of bacterial cells relative to the tick (Samuels et al., 2021). Although microarray approaches have been used to identify large scale changes in Bb gene expression between the first and second tick bloodmeals (Iyer et al., 2015), we have limited temporal resolution into the molecular transitions that happen across key steps of tick feeding. This problem necessitates novel approaches to capture the transcriptomic changes of Bb within the natural tick environment.

While the full landscape of Bb transmission determinants is not yet known, we do have a growing knowledge of functional processes that are critical during the tick stage, such as motility, metabolism, and immune evasion (Kurokawa et al., 2020; Phelan et al., 2019). These functions often rely on the unique protein-rich Bb outer surface. Notably, several specific tick–Bb protein– protein interactions are important for Bb survival, migration, or transmission to the next host. Bb encodes an extensive outer surface protein (Osp) family with members that are differentially expressed during host switching. One of these proteins, OspA, binds a tick cell surface protein, tick receptor for OspA (TROSPA), which is required for successful Bb colonization of the tick midgut during the first acquisition bloodmeal (Pal et al., 2004). Several other proteins have also been linked to Bb migration within the tick (Pal et al., 2022). For example, BBE31 binds a tick protein TRE31, and disruption of this interaction decreases the number of Bb cells that successfully migrate from the tick gut to salivary glands (Zhang et al., 2011). However, these interactions alone are not sufficient to block Bb growth or migration in ticks, suggesting there are likely additional molecular factors from Bb and ticks at play during tick-borne transmission.

To provide a more comprehensive set of Bb determinants driving tick-borne transmission of this important human pathogen, we longitudinally mapped Bb genome-wide expression changes for bacterial cells isolated from inside ticks during the transmission bloodmeal. We developed a novel sequencing-based strategy for ex vivo transcriptomic profiling of Bb populations within feeding ticks. We identified 192 highly differentially expressed genes, including genes previously implicated in Bb transmission as well as many others. Genes upregulated during tick transmission were enriched for outer surface lipoproteins, suggesting Bb dramatically remodels its cell envelope as it migrates through the tick. Mass spectrometry analyses revealed dramatic changes in the tick environment over feeding, identifying new potential determinants of a more extensive and diverse set of tick–microbe molecular interactions than previously appreciated.

Results

A two-step enrichment process facilitates robust transcriptional profiling of Bb during the tick bloodmeal

To gain a more comprehensive understanding of Bb gene expression throughout the tick phase of the transmission cycle, we developed an experimental approach to characterize the Bb transcriptome of spirochetes isolated from ticks during a bloodmeal. We aimed to establish a longitudinal transcriptional profile encompassing key pathogen transmission events each day of feeding after ticks attached to their mouse bloodmeal hosts (Figure 1A). The major bottleneck for such an RNA sequencing (RNA-seq) approach is capturing sufficient quantities of Bb transcripts from complex multi-organism samples in which pathogen transcripts represent a very small minority. Our initial attempts to uncover Bb mRNA by simply removing tick mRNA with polyA-depletion and removing tick rRNA sequences using Depletion of Abundant Sequences by Hybridization (DASH) (Dynerman et al., 2020; Gu et al., 2016) were unsuccessful. This approach resulted in an average of only 0.09% of RNA-seq reads mapping to Bb mRNA – approximately 10-fold less than we estimated would be needed to feasibly obtain robust transcriptome-wide gene expression analysis (Haas et al., 2012).

Figure 1.

To dramatically increase Bb transcript representation in our libraries, we physically enriched Bb cells in tick lysates prior to library preparation by adding an initial step of immunomagnetic separation (Figure 1B). We took advantage of a commercial antibody previously generated against whole Bb cells (αBb, RRID: AB_1016668). By western blot analysis, we confirmed that αBb specifically recognized several Bb proteins, including surface protein OspA, which is highly prevalent on the Bb surface in ticks (Figure S1A). In addition, immunofluorescence microscopy with αBb showed clear recognition of Bb cells from within the tick at each day of feeding (Figure S1B). After collecting infected ticks from mice one, two, three, and four days post-attachment, we used αBb and magnetic beads to enrich Bb cells away from the majority of the tick material in the lysates. We tracked relative Bb enrichment through RT-qPCR of Bb flaB RNA and tick gapdh RNA in the separated samples. Measuring Bb flaB RNA from both Bb-enriched samples and their matched Bb-depleted fractions, we found over 95% of total Bb flaB RNA was present in enriched fractions (Figure 1C), suggesting our approach captured the vast majority of Bb transcripts from the tick.

Total RNA recovered from the enrichment process was used to create RNA-seq libraries that subsequently underwent depletion of rRNA sequences from ticks, mouse, and Bb using DASH, which targets unwanted sequences for degradation by Cas9 (Gu et al., 2016; Ring et al., 2022). DASH reduced unwanted sequences from 94% to 9% of our total libraries, greatly increasing the relative abundance of Bb transcripts (Figure 1D). For resulting libraries generated across feeding (Table S1), between 0.6% and 3.4% of reads mapped to Bb coding regions (Figure 1E), which translated to an average of 4.3 million Bb reads per sample (Figure 1F). Read coverage was sufficient for downstream differential expression analyses for the vast majority of Bb genes, with at least 10 reads mapping to 92% of Bb genes. To evaluate whether our approach introduced any major artifacts, we sequenced and compared RNA-seq libraries from in vitro cultured Bb cells before and after immunomagnetic enrichment. We found minimal expression differences (29 genes with p<0.05, fold changes between 0.83-1.12) (Figure S2, Table S2), suggesting experimental enrichment did not significantly alter global transcriptome profiles for Bb. Thus, we developed an approach to enable genome-wide analysis of Bb population-level expression changes that occur within the feeding tick during transmission.

Global ex vivo profiling of Bb reveals extent and kinetics of transcriptional changes

To provide a broad overview of Bb expression changes during feeding, we performed principal component analysis (PCA) on the Bb transcriptome data from one, two, three, and four days post-attachment (n=4). We reasoned that if many longitudinal expression changes were occurring across Bb populations, we would observe greater data variability between time points than between biological replicates. Indeed, we found replicates from each day grouped together, whereas distinct time points were largely non-overlapping. The first principal component, which explained 64% of the variance in our data, correlated well with day of feeding (Figure 2A). The global pattern suggested that Bb gene expression changes generally trended in the same direction over the course of feeding with the most dramatic differences between flanking timepoints on day 1 and day 4.

Figure 2.

Using day 1 (early attachment) as a baseline, we performed differential expression analysis for all Bb genes at subsequent time points and examined changes above a twofold threshold (Figure 2B, Table S3). These analyses mirrored the global longitudinal expression pattern predicted by our PCA. The total number of differentially expressed (DE) genes increased steadily over time and peaked at day 4, which when compared to day 1 had 186 DE genes, including 153 upregulated and 33 downregulated (Figure 2C-E). Across time points, DE genes were highly overlapping and largely changed in the same directions. In all three comparisons, we found 192 DE genes in total (Table S4). However, there were some notable differences between genes, such as the overall timing and kinetics of expression changes. Transcript levels for some DE genes changed suddenly over the course of feeding, while others were more gradual. To our knowledge, this is the first comprehensive report of global Bb expression changes within feeding ticks.

To help assess the integrity of our dataset, we examined expression profiles for previously characterized targets of major transcriptional programs activated at the onset of the bloodmeal: the Hk1/Rrp1 two-component system and the RpoN/RpoS sigma factor cascade (Caimano et al., 2019, 2015). Between day 1 and day 4, the majority of genes previously reported to be activated or repressed by Rrp1 in vitro (Caimano et al., 2015) trended significantly in the expected direction ex vivo (111/148 upregulated, 37/57 downregulated, p<0.05, Wald tests, Figure S3A). The majority of the genes previously reported to be activated by RpoS in a mammalian context, including canonical targets ospC and dbpA, were also significantly upregulated during tick feeding (33/55, p<0.05, Wald tests) (Figure S3B), along with rpoS itself (Figure S3C), as expected (Hübner et al., 2001). However, very few genes repressed by RpoS in mammals were similarly repressed in ticks (2/33, p<0.05, Wald tests). These data support the hypothesis that while RpoS-dependent gene activation occurs across both vertebrate and vector stages of Bb transmission, RpoS-dependent repression is specific to spirochetes infecting a mammalian host (Caimano et al., 2019, 2007). We also compared our DE genes to those regulated by Rel_Bbu as part of the stringent response, another major transcriptional program active in Bb in the tick during nutrient starvation (Drecktrah et al., 2015). About half of Rel_Bbu-regulated genes changed in the direction expected if the stringent response was active during this time (Figure S3D); however, more of the genes that were downregulated over feeding have been reported to be regulated by Rel_Bbu than either RpoS or Rrp1, suggesting it may play some role during feeding.

By looking at these previously reported transcriptional regulons in our data, we were able to verify that our data captured expected transcriptional programs active during tick feeding. Nevertheless, 25% of the twofold DE genes were not previously identified through these RNA-seq studies as dependent upon RpoS, Rrp1, or Rel_Bbu (Caimano et al., 2019, 2015; Drecktrah et al., 2015), which are the three main Bb regulatory programs active in the tick (Samuels et al., 2021). These additional genes highlight the necessity of measuring transcription in the tick environment and suggest we uncovered gene expression changes specific to the tick-stage of the Bb enzootic cycle. The nature and dynamics of these changes provide insights into potential genetic determinants of Bb survival, proliferation, and dissemination in the tick during transmission.

Bb genes upregulated during feeding are found predominantly on plasmids

Bb has a complex, highly fragmented genome (Barbour, 1988) (Figure 3A), including numerous plasmids that are necessary during specific stages of the enzootic cycle (reviewed in Schwartz et al., 2021) suggesting they contain genes that are crucial for pathogen transmission and survival. In fact, many genes found on the plasmids have been previously shown to alter expression upon environmental changes or in different host environments (Iyer et al., 2015; Ojaimi et al., 2005; Revel et al., 2002; Tokarz et al., 2004). Thus, we reasoned that many of the 192 Bb genes that change expression over the course of the tick bloodmeal would reside on the plasmids, and we examined their distribution throughout the genome. Consistent with these previous reports, we found that most of the upregulated genes were located on the plasmids (143/158; 90%), while fewer were found on the chromosome (15/158; 10%) (Figure 3B), which is home to the majority of metabolic and other housekeeping genes. In contrast, the majority of the downregulated genes were found on the chromosome (27/34, 79%) (Fig 3C).

Figure 3.

Several plasmid-encoded genes that were longitudinally upregulated in our datasets have known roles during the tick bloodmeal or in mammalian infection. For example, linear plasmid 54 (lp54), which is an essential plasmid present in all Bb isolates (Casjens et al., 2012), contained the largest number of upregulated genes. Many of the genes on lp54 are known to be regulated by RpoS during feeding, including those encoding adhesins DbpA and DbpB, which are important for infectivity in the host (Blevins et al., 2008). This set also included five members of a paralogous family of outer surface lipoproteins BBA64, BBA65, BBA66, BBA71, and BBA73. BBA64 and BBA66 are necessary for optimal transmission via the tick bite (Gilmore et al., 2010; Patton et al., 2013). These findings indicate our dataset captures key Bb transcriptional responses known to be important for survival inside the tick during a bloodmeal.

Many upregulated genes were encoded by cp32 plasmid prophages. Bb strain B31-S9 harbors seven cp32 isoforms that are highly similar to each other (Casjens et al., 2012). When cp32 prophages are induced, phage virions called ϕBB1 are produced (Eggers and Samuels, 1999). In addition to phage structural genes, cp32 contain loci that encode various families of paralogous outer surface proteins (Stevenson et al., 2000). Amongst the cp32 genes that increased over feeding were members of the RevA, Erp, and Mlp families, which are known to increase expression during the bloodmeal (Gilmore et al., 2001). We also found several phage genes that were upregulated, including those encoding proteins annotated as phage terminases on cp32-3, cp32-4, and cp32-7 (BBS45, BBR45, and BBO44). Some cp32 genes have been shown to change expression in response to temperature changes (Tokarz et al., 2004) and as a part of the stringent response regulated by Rel_Bbu (Drecktrah et al., 2015). Our data suggest that some prophage genes are upregulated over the course of tick feeding, raising the possibility that cp32 prophage are induced towards the end of feeding. Overall, our data support the long-held idea that the Bb plasmids, which house many genes encoding cell envelope proteins, proteins of unknown function, and prophage genes, play a critical role in the enzootic cycle during the key transition period of tick feeding.

Bb genes encoding outer surface proteins are enriched among upregulated genes

To gain a better overall sense of the types of genes that changed over feeding and the timing of those changes, we grouped DE genes into functional categories. Since a high proportion of plasmid genes encode lipoproteins within the unique protein-rich outer surface of Bb, genes of unknown function, and predicted prophage genes (Casjens et al., 2000; Fraser et al., 1997), we expected that many of the DE genes would fall into these categories. We classified the genes as related to either: cell envelope, bacteriophage, cell division, DNA replication and repair, chemotaxis and motility, metabolism, transporter proteins, transcription, translation, stress response, protein degradation, or unknown (as in Drecktrah et al., 2015)(see Table S4).

Of the genes that increased at least twofold over feeding, the clear majority each day of feeding fell into two broad categories, cell envelope (55 of 158 total genes) and proteins of unknown function (69 of 158 total genes), with fewer genes related to metabolism, chemotaxis and motility, transporters, bacteriophage, cell division, and transcription (Figure 4A). In contrast to the upregulated genes, genes downregulated during feeding were more evenly distributed among the functional categories – including translation, protein degradation, transcription, and metabolism – consistent with many of them being located on the chromosome.

Figure 4.

When looking at changes across these functional categories, the overrepresentation of cell envelope proteins was striking, while not unexpected. The Bb outer surface is covered with lipoproteins (Figure 4B), and these proteins are critical determinants in Bb interactions with the various environments encountered during the enzootic cycle (Kurokawa et al., 2020). We found that more than half of these outer surface lipoproteins changed expression during this time, with 46 of 83 total outer surface lipoproteins (Dowdell et al., 2017) differentially expressed twofold over feeding (Figure 4C). These data suggested widespread changes may be occurring on the Bb outer surface during feeding.

To understand the functional implications of expression changes in a majority of outer surface lipoproteins, we also compared their relative expression. The magnitude of expression varied greatly, with ospA, ospB, ospC, and bba59 being the most highly expressed outer surface protein transcripts. Many of the outer surface protein genes that we found increased expression over feeding were much less abundantly expressed (Table S5). However, even the genes that appeared to have low expression in these population level measurements could play important roles in transmission if they are highly expressed in a small number of crucial cells, such as those that ultimately escape the midgut. While bulk RNA-seq cannot distinguish what is happening at the single cell level, our data suggest that during the bloodmeal, Bb are undergoing a complex outer surface transformation driven by increases in a majority of the genes encoding these crucial lipoproteins.

Identification of candidate tick interaction partners of Bb cells ex vivo

Our RNA-seq data suggested that the outer surface of Bb transforms over the course of feeding as Bb are primed for transmission to a vertebrate host. At the same time, the tick midgut environment is changing as the tick begins to digest its bloodmeal. Tick-Bb interactions are likely crucial at the beginning of the tick bloodmeal, as Bb adhere to the tick gut epithelium before becoming motile and migrating out of the midgut and into the salivary glands (Dunham-Ems et al., 2009). Some Bb outer surface proteins, such as BBE31 and BBA52 play key roles in pathogen migration through interactions with the tick environment (Kurokawa et al., 2020). We wanted to explore the changing tick environment to identify tick proteins with which Bb could interact throughout the tick bloodmeal. Since our Bb enrichment process retained some tick material, we reasoned that tick proteins that interact with Bb would be present in these samples.

To outline the changes occurring in the tick during feeding and to identify candidate Bb-interacting tick proteins, we used mass spectrometry to survey the content of the tick material that was enriched along with the Bb cells we sequenced at early and late stages during feeding. We purified proteins from the αBb-enriched fraction of crushed infected ticks on day 1 after attachment and day 4 after attachment in triplicate. As controls for each day, we also performed the Bb enrichment process on lysate from uninfected ticks mixed with in vitro cultured Bb to help rule out proteins that were not pulled down through in vivo tick-Bb interactions (Figure 5A). We identified between 414 and 2240 protein groups per sample replicate. To identify proteins of interest, we looked for those that were detected in at least two of three replicates with infected ticks and had a mean average coverage twice that of uninfected ticks mixed with cultured Bb. We found 256 proteins that were enriched with Bb from infected ticks one day after attachment and 226 proteins that were enriched with Bb from infected ticks four days after attachment (Table S6). Of these proteins, only 27 were detectable on both days, suggesting the tick proteins present upon Bb enrichment change dramatically over the course of feeding. The vast majority of all detected proteins were from ticks (98%). Amongst the small number of Bb proteins, we identified OspC in the samples from day 4 after attachment but not day 1 after attachment, confirming the expression change we saw in our RNA-seq. The distinct sets of tick proteins we identified at each timepoint suggest dramatic changes occur in the Bb-infected tick midgut environment during feeding that may alter the landscape of tick-Bb interactions.

Figure 5.

Some of the proteins enriched with Bb from the tick may be good candidates for key Bb-interacting partners during feeding, especially if they are localized to the surface of tick cells where they may encounter Bb. The scarcity of both predicted and experimentally validated functions and localizations for tick proteins makes it difficult to fully assess the potential for tick protein interactions with extracellular Bb. Nevertheless, of the proteins found exclusively one day after attachment, 10 were categorized as extracellular matrix proteins using the PANTHER gene database (Thomas et al., 2022), and this category was statistically enriched (Fisher’s exact test, FDR=0.000093) (Figure 5B, Table S7). 30 additional proteins were annotated with a cellular component as plasma membrane. Four days after attachment, we did not detect any annotated extracellular matrix proteins; however, we identified 31 proteins that are likely to be found at the membrane, including two proteins annotated as putative low-density lipoprotein receptors (Figure 5C, Table S8). These extracellular matrix and membrane proteins may be the most likely to directly interact with Bb during this timeframe and are candidates for tick proteins important in the Bb dissemination process. The proteins present in the changing tick environment may be key determinants of pathogen transmission as Bb remodels its outer surface while preparing to migrate through the tick to a new host.

Discussion

Vector-borne pathogens must adapt to distinct environments as they are transmitted by arthropods to colonize new bloodmeal hosts. The tick-borne Lyme disease pathogen Bb undergoes transcriptional changes inside its vector over the course of a single bloodmeal that are important for a number of key transmission events. For example, expression changes enable survival in the feeding tick (He et al., 2011), facilitate transmission across several internal compartments into the bloodmeal host (Kurokawa et al., 2020), and prime bacteria for successful infection of the next host (Kasumba et al., 2016). Capturing these changes as they occur in vivo has been challenging due to the low relative abundance of Bb material inside of the rapidly growing, blood-filled tick. To date, many advances in our understanding of Bb expression during transmission have come from tracking changes in small subsets of genes during tick feeding (Gilmore et al., 2001), leveraging Bb culture conditions that approximate environmental changes across its lifecycle (Ojaimi et al., 2005; Revel et al., 2002; Tokarz et al., 2004), and defining transcriptional regulons outside of the tick (Caimano et al., 2019, 2015; Drecktrah et al., 2015). Here, we developed an experimental RNA-seq-based strategy to more directly and longitudinally profile gene expression for Bb populations isolated from ticks over the course of a transmitting bloodmeal.

Our longitudinal collection of transcriptomes for Bb cells isolated from ticks serves as a resource and starting point for delineating the functional determinants of Bb adaptation and transmission. Our transcriptome profiles reveal more extensive expression changes for Bb outer surface proteins than previously appreciated. These results raise the possibility that Bb is actively remodeling its outer coat to navigate the dynamic tick environment during feeding, akin to wardrobe changes across seasons. These outer surface changes may play roles in the persistence of Bb in the tick, cell adhesion, or in immune evasion, either inside of the tick or later in the vertebrate host (Kenedy et al., 2012). It has long been observed that molecular interactions between pathogens and the midgut of their vectors are key determinants of transmission (Barillas-Mury et al., 2022). Modifications to the Bb surface could facilitate different tick–pathogen interactions that are critical for its physical movement through tick compartments. Notably, we also identified numerous genes of unknown function that change expression during feeding. Although challenging to address, efforts focused on uncovering the contribution of these genes to the Bb life cycle could be groundbreaking for understanding the unique aspects of tick-borne transmission.

More comprehensive knowledge about Bb outer surface proteins and their functional consequences could lead to new avenues for curtailing the spread of Lyme disease. Protective vaccines against vector-borne pathogens have often targeted surface proteins encoded by pathogens (Kovacs-Simon et al., 2010). While previous Lyme disease vaccine efforts effectively targeted the highly expressed Bb outer surface protein OspA (Steere et al., 1998), more targets could increase the likelihood of a successful vaccine. Our study provides a catalog of new Bb candidates that could be explored. In addition, our biochemical pull-downs using tick-isolated Bb cells as bait unearthed a preliminary list of potential tick proteins that could be involved in tick–Bb molecular interactions. Blocking molecular interactions that are functionally critical for transmission could also be explored as a therapeutic strategy (Barillas-Mury et al., 2022; Manning and Cantaert, 2019).

Our transcriptomic resource provides critical insights into Bb population-level changes during the vector stage of its lifecycle, an important starting point for understanding the primary drivers of tick-borne transmission. Strikingly few Bb cells out of the total pathogen population in ticks are ultimately transmitted to the next host during feeding. There are several major bottlenecks to infection, including dissemination to the tick salivary glands and survival upon inoculation of the host (Dunham-Ems et al., 2009; Rego et al., 2014). There may be important physiological or molecular variations across Bb cells within the population residing in the tick that contribute to these differential outcomes. An important clue is the fact that Bb proteins are indeed heterogeneously expressed across cells within the feeding tick midgut (Ohnishi et al., 2001). Conducting transcriptomic analyses at the single-cell level will be key. Given that only a small number of Bb escape the tick midgut, genes that appear unaltered or show low expression levels across the bulk population could still play an outsized role in infection for a minority of the cells. Our work provides a foundational methodology that can be leveraged to greatly improve the resolution of tick–microbe studies. Technological advances stemming from this work will provide molecular details that potentiate new questions and may unearth surprising mechanistic insights into the unique lifestyles of tick-borne pathogens.

Methods

B. burgdorferi culture

Bb strain B31-S9 (Rego et al., 2011) was provided by Dr. Patricia Rosa (NIAID, NIH, RML) and cultured in BSK II media at 35°C, 2.5% CO₂.

Tick feeding experiments

I. scapularis larvae were purchased from the Tick Lab at Oklahoma State University (OSU) (for RNA-seq experiments) or provided by BEI Resources, a division of the Center for Disease Control (for mass spectrometry experiments). Before and after feeding, ticks were maintained in glass jars with a relative humidity of 95% (saturated solution of potassium nitrate) in a sealed incubator at 22°C with a light cycle of 16h/8h (light/dark). Animal experiments were conducted in accordance with the approval of the Institutional Animal Care and Use Committee (IACUC) at UCSF. Ticks were fed on young (4–6-week-old) female C3H/HeJ mice acquired from Jackson Laboratories. Mice were anesthetized with ketamine/xylazine before placement of ≤ 100 larval or ≤ 30 nymphal ticks. Ticks were either pulled off isoflurane anesthetized mice at various times during feeding (1-3 days after placement) or allowed to feed to repletion and collected from mouse cages. Replete larval ticks were placed in the incubator to molt before being used as nymphs in experiments.

Western blot with anti-Bb antibody

To determine whether the anti-Bb antibody targeted ospA, we cultured wildtype Bb (B31-A3), ospA1-mutants (ospA1) and ospA-restored Bb (ospA⁺B1) (Battisti et al., 2008) to approximately 5×10⁷ Bb/mL and centrifuged 3 mLs of culture for 7 minutes at 8000 × g, washing twice with PBS. Pelleted cells were lysed in 50 μL of water, and 25 ng per sample was mixed with 5X loading dye (0.25% Bromophenol Blue, 50% Glycerol, 10% Sodium Dodecyl Sulfate, 0.25M Tris-Cl pH 6.8, 10% B-Mercaptoethanol), run on a Mini-PROTEAN TGX 4-15% gel (Biorad), and transferred using the Trans-Blot Turbo Transfer System (Biorad). After transfer, the blot was blocked for 30 minutes at 4°C in TBST (Tris buffered saline with .1% tween) with 5% milk, then treated with anti-Bb antibody (Invitrogen: PA1-73004; RRID: AB_1016668) diluted 1:10,000 for 1 hour at room temperature, followed by anti-rabbit HRP secondary antibody (Advansta) diluted 1:5,000 for 45 minutes at room temperature with short PBST washes between each step. Blots were exposed using Clarity Western ECL Substrate (Biorad) and imaged using the Azure C400 imaging system (Azure Biosystems). This experiment was repeated three times.

Enrichment of Bb from feeding ticks

To sequence RNA from Bb from feeding ticks, we enriched Bb to increase the ratio of Bb to tick material. Larval ticks were fed to repletion on three mice that were infected with Bb through intraperitoneal and subcutaneous injection with 10⁴ total Bb. Approximately five months later, the molted nymphal ticks were fed on eight mice. We estimated that 83% of the ticks were infected with Bb by crushing 12 unfed nymphs in BSK II media and checking for viable Bb days later. Ticks were pulled from all mice and pooled into four replicates 1 day after placement (14 ticks per replicate), 2 days after placement (12 ticks per replicate), and 3 days after placement (6 ticks per replicate) and collected from cages 4 days after placement (7 ticks per replicate). Shortly after collection, ticks were washed with water and placed in a 2 mL glass dounce grinder (Kimble) in 500 μL of phosphate-buffered saline (PBS). Ticks were homogenized first with the large clearance pestle and then the small clearance pestle. The homogenate was transferred to a 1.5 mL Eppendorf tube and 500 μL of PBS was added to total 1 mL. At this stage, 50 μL of homogenate was removed as an input sample and mixed with 500 μL of TRIzol (Invitrogen) for RNA extraction. 2 μL of anti-Bb antibody (Invitrogen: PA1-73004; RRID: AB_1016668) was added to the homogenate, which was then placed on a nutator at 4°C for 30 minutes. During incubation, 50 μL of Dynabeads Protein G (Invitrogen) per sample were washed twice in PBS. After incubation with the antibody, the homogenate and antibody mixture were added to the beads. This mixture was placed on a nutator at 4°C for 30 minutes. Tubes were then placed on a magnet to secure beads, and the homogenate was removed and saved to create depleted samples. The depleted homogenate was centrifuged at 8,000 × g for 7 minutes, 900 μL of supernatant was removed, and 500 μL of TRIzol was added to the pellet to create depleted samples. The beads were washed twice with 1 mL of PBS, resuspending the beads each time. The second wash was removed and 500 μL of TRIzol was added to the beads to create enriched samples. RNA was extracted from all input, enriched, and depleted samples using the Zymo Direct-zol RNA Microprep Kit with on-column DNase treatment (Zymo Research). The step-by-step Bb enrichment protocol is available at: dx.doi.org/10.17504/protocols.io.36wgqjrbovk5/v1.

Enrichment of Bb from culture

To test whether the Bb enrichment process altered gene expression levels, we performed the enrichment protocol on cultured Bb. We grew three tubes of Bb in BSK II media to 9×10⁴ Bb/mL at 35°C. 1 mL of culture was spun down at 8,000 × g for 7 minutes, media was removed, Bb were washed in 1 mL of PBS and spun again. Pelleted Bb were resuspended in 1 mL of fresh PBS. These samples were used as starting homogenate for the Bb enrichment protocol and input, enriched, and depleted fractions were collected as above. RNA-seq libraries from these samples were prepared and sequenced as below.

RNA-seq library preparation and sequencing

To make RNA-seq libraries from enriched Bb RNA, we used 50 ng of total RNA as input into the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England BioLabs). We prepared libraries following the manufacturer’s protocol for use of the kit with purified mRNA or rRNA-depleted RNA, despite starting with total RNA. Libraries were barcoded using NEBNext Multiplex Oligos for Illumina Dual Index (New England BioLabs).

To deplete reads from the libraries that were from tick, Bb, or mouse rRNA, we used Depletion of Abundant Sequences by Hybridization (DASH) (Gu et al., 2016), which targets Cas9 to unwanted reads in RNA-seq libraries using custom dual-guide RNAs (dgRNAs). Our dgRNAs targeted short sequences within tick rRNA, mouse rRNA, and Bb rRNA that were designed using DASHit software (Dynerman et al., 2020). We ordered crRNAs that targeted our sequences of interest (Ring et al., 2022). To transcribe the dgRNAs, we followed our protocol for In Vitro Transcription for dgRNA V2 (Lyden et al., 2019b). Both tracrRNA and pooled crRNA DNA templates were annealed to equimolar amounts of T7 primer by heating to 95°C for 2 minutes and slowly cooling to room temperature. Annealed templates were used in 1 mL in vitro transcription reactions with: 120 μL10X T7 buffer (400 mM Tris pH 7.9, 200 mM MgCl2, 50 mM DTT, 20 mM spermidine (Sigma)), 100 μL of T7 enzyme (custom prepped enzyme from E. Crawford, diluted 1:100 in T7 buffer, final concentration: 100 μg/mL), 300 μl NTPs (25 mM each, Thermo Fisher Scientific), 4 μg of annealed crRNA template or 8 μg of annealed tracrRNA template, and water to 1 mL. In vitro transcription was performed for 2 hours at 37°C. Reactions were purified twice with the Zymo RNA Clean & Concentrator-5 Kit (Zymo Research). To form the dgRNA complex for DASH, crRNA and tracrRNAs were diluted to 80 μM, mixed in equimolar amounts and annealed by heating to 95°C for 30 seconds and cooling slowly to room temperature.

After transcription of dgRNAs, we performed DASH Protocol Version 4 (Lyden et al., 2019a) on each library individually. Cas9 and transcribed dgRNAs were prepped by mixing: 2.5 μL 10X Cas9 buffer, 5 μL Cas9 (New England BioLabs), and 5 μL of 40 μM transcribed dgRNAs. The mixture was incubated at 37°C for 5 minutes before 7.5 μL of RNA-seq library (2.8 nM) was added. The mixture was incubated at 37°C for 1 hour and then purified with Zymo DNA Clean & Concentrator-5 (Zymo Research) following the PCR product protocol and eluting DNA into 10.5 μL of water.

During cleanup, Cas9 was again mixed with buffer and dgRNAs and incubated at 37°C for 5 minutes. Following cleanup, eluted DNA was added to the second Cas9-dgRNA mixture and incubated at 37°C for 1 hour for a second time. Then, 1 μL of proteinase K (New England Biolabs) was added, and the mixture was incubated at 50°C for 15 min. The libraries were then purified with 0.9x volume of sparQ PureMag Beads (QuantaBio) following the standard protocol, eluting in 24 μL of water. rRNA-depleted RNA-seq libraries were then amplified in a BioRad CFX96 using the Kapa HiFi Real-Time Amplification Kit in a 50 μL reaction with 25 μL master mix, 23 μL of the DASHed library pool, and 2 μL of 25 μM mix of Illumina P5 (5’-AATGATACGGCGACCACCGAGATCT) and P7 (5’-CAAGCAGAAGACGGCATACGAGAT) primers. The qPCR program for amplification was as follows: 98°C for 45 sec (1 cycle), (98°C for 15 sec, 63°C for 30 sec, 72°C for 45 sec, plate read, 72°C for 20 sec) for 10 cycles (day 3 and day 4 samples) or 11 cycles (day 1 and day 2 samples). The libraries were removed from cycling conditions before leaving the exponential phase of amplification and then purified with 0.9X volume of sparQ PureMag Beads according to the standard protocol.

Following DASH, RNA-seq libraries were sequenced on an Illumina NovaSeq S2 (2 lanes) with paired-end 100 base pair reads. Libraries from in vitro cultured control experiment were sequenced on an Illumina NextSeq with paired-end 75 base pair reads. FASTQ files and raw Bb read counts for in vitro control experiment (GSE217146) and ex vivo experiment (GSE216261) have been deposited in NCBI’s Gene Expression Omnibus (Edgar et al., 2002) under SuperSeries accession number GSE217236.

RNA-seq data analysis

To measure the success of our rRNA depletion through DASH, we used DASHit software (Dynerman et al., 2020) to determine the percentage of reads that would be DASHable by our guide RNAs (Ring et al., 2022). For pre-DASH data, we sequenced the input of each of our RNA-seq libraries before performing DASH on a MiSeq V2 Micro (Illumina). We tested DASHability on a random subset of 200,000 paired-end reads chosen by seqtk (https://github.com/lh3/seqtk) from each pre- and post-DASH library.

To map our RNA-seq data to Bb, we wanted to optimize for mapping reads that came from the many paralogous gene families found across the plasmids of the genome. We used the pseudoalignment tool Salmon v1.2.1 (Patro et al., 2017), which is used to accurately map reads coming from different isoforms of the same gene, for this reason. While using Salmon to map to CDS sequences may improve mapping to paralogous genes, it may also have a tradeoff of reduced mapping of reads that fall on the ends of genes that reside in operons. Nevertheless, all samples should be similarly affected, and any undercounting should not change differential expression results. Reads were first trimmed of bases with quality scores less than 20 using Trim Galore v0.6.5 (https://github.com/FelixKrueger/TrimGalore). We mapped our reads to Bb CDS sequences as a reference transcriptome: NCBI Genbank GCA_000008685.2 ASM868v2 (with plasmids lp5, cp9, and lp56 removed as they are not present in B31-S9) using Salmon with the following parameters: --validateMappings --seqBias --gcBias. Before mapping, the transcriptome was first indexed using the Salmon index command with the whole genome as decoys and the parameter --keepDuplicates to keep all duplicate genes.

Read counts from Salmon were used as input into DESeq2 v1.24.0 (Love et al., 2014) for differential expression analysis in R version 3.6.1. DESeq2 function PlotPCA() was used to create a PCA plot from read counts after running the varianceStabilizingTransformation() function. For differential expression analysis between days, a DESeq object was created from count data using the DESeq() function. The lfcShrink() function with the apeglm method (Zhu et al., 2018) was used to calculate fold changes between days. Code used for differential expression analysis is available at https://github.com/annesapiro/Bb-tick-feeding.

Gene classification

To classify genes into functional groups, functional categories were sourced from Drecktrah et al., 2015 where available. Other gene functions were sourced from Fraser et al., 1997. Genes found within the co-transcribed “late” bacteriophage operon (Zhang and Marconi, 2005) were considered “bacteriophage” even if their function is unknown. Outer surface lipoproteins were those found in Dowdell et al., 2017 plus additional lipoproteins listed in Iyer et al., 2015 that were found in Dowdell et al. 2017 categories SpII and SpI and evidence of outer surface localization. Outer surface and periplasmic lipoproteins were classified as “cell envelope” in the absence of other classifications. Gene family information from Casjens et al., 2000 was considered to aid in classification.

RT-qPCR measuring Bb enrichment

To test the efficacy of our enrichment protocol, we used RT-qPCR to quantify Bb flaB and I. scapularis gapdh transcript levels in enriched and depleted fractions. We synthesized cDNA from 8 μL of RNA extracted from Bb enrichment samples and their matched depleted samples from feeding day 2, day 3, and day 4 using the qScript cDNA Synthesis Kit (Quantabio). cDNA was diluted 2X before use in qPCR. To measure flaB copies, we created standards of known concentration from purified PCR products. These standards were made from PCR with primers with the following sequences: 5’-CACATATTCAGATGCAGACAGAGGTTCTA and 5’-GAAGGTGCTGTAGCAGGTGCTGGCTGT. We ran a dilution series with ten-fold dilutions between 10⁶ copies and 10¹ copies of this PCR template alongside our enriched and depleted samples. We performed qPCR using Taqman Universal PCR Master Mix (Applied Biosystems). The primers used to amplify flaB were: 5’-TCTTTTCTCTGGTGAGGGAGCT and 5’-TCCTTCCTGTTGAACACCCTCT (used at 900 nM) and the probe was /56-FAM/AAACTGCTCAGGCTGCACCGGTTC/36-TAMSp (used at 250 nM). For tick gapdh RT-qPCR, the cDNA samples were diluted an additional 2X. Standards of known concentration were created using the qPCR primer sequences. Primer sequences were: 5’-TTCATTGGAGACACCCACAG and 5’-CGTTGTCGTACCACGAGATAA (used at 900 nM). qPCR was performed using PowerUp SYBR Green Master Mix (Applied Biosystems). For both flaB and gapdh, the number of copies in each sample was calculated based on the standards of known concentration. Three technical replicates were averaged from each of four biological replicates at each time point tested. We totaled the number of copies in each matched enriched and depleted fraction to calculate the fraction of flaB or gapdh that was found in either sample. All qPCR was performed on the QuantStudio3 Real-Time PCR System (Applied Biosystems).

Immunofluorescence microscopy

To test whether the anti-Bb antibody recognized Bb inside of the tick, ticks at each day of feeding were crushed in 50 μL of PBS. 10 μL of lysate was spotted onto slides and allowed to air dry before slides were heated briefly three times over a flame. Heat fixed slides were then treated with acetone for one hour. Slides were incubated with anti-Bb primary antibody (1:100 diluted in PBS + 0.75% BSA) for 30 minutes at 37°C in a humid chamber. A control without primary antibody was also used for each day. Slides were washed once in PBS for 15 minutes at room temperature, then rinsed in distilled water and air dried. Anti-rabbit IgG Alexa 488 (Invitrogen) diluted 1:100 in PBS + 0.75% BSA was added for 30 minutes at 37°C in a humid chamber. Slides were washed in PBS for 15 minutes at room temperature three times, adding 1:100 Propidium Iodide (Invitrogen) during the second wash. Slides were then rinsed with distilled water and air dried before the addition of mounting media (Fluoromount-G, SouthernBiotech) and cover slips. Fluorescence imaging was performed on a Nikon Ti2 inverted microscope for widefield epifluorescence using a 100x/1.40 objective. Images were captured with NIS-Elements AR View 5.20 and then processed with ImageJ software (Schneider et al., 2012). No strong florescence signal was observed on the control slides without primary antibody.

Mass spectrometry of Bb-enriched samples

To identify which tick proteins were found in samples after Bb enrichment across feeding, we fed both uninfected and infected ticks on mice. We collected three replicates of uninfected ticks one day after attachment (11 ticks per replicate), infected ticks one day after attachment (27 ticks per replicate), uninfected ticks four days after attachment (8 ticks per replicate), and infected ticks four days after attachment (16 ticks per replicate). Before anti-Bb enrichment, the uninfected tick samples were mixed with Bb grown in culture that was washed with PBS (3×10⁴ Bb one day after attachment and 3×10⁶ Bb four days after attachment) and mixed lysates were rotated at room temperature for 30 minutes. Infected tick samples underwent the Bb enrichment process immediately. The enrichment process followed the same protocol used for RNA-seq. Sample volumes were increased to 1 mL as needed, and then 2 μL of anti-Bb antibody (Invitrogen: PA1-73004; RRID: AB_1016668) was added, and samples were rotated at 4°C for 30 minutes. 50 μL of Dynabeads Protein G per sample were washed in PBS during this incubation and added to the lysates, which were rotated at 4°C for 30 minutes. The beads were washed twice with 1 mL of PBS, and then placed into 50 μL of lysis buffer (iST LYSE, PreOmics). Samples were boiled at 95°C for 5 minutes, and lysates were removed from beads and frozen for mass spectrometry preparation.

For mass spectrometry, a nanoElute was attached in line to a timsTOF Pro equipped with a CaptiveSpray Source (Bruker). Chromatography was conducted at 40°C through a 25 cm reversed-phase C18 column (PepSep) at a constant flowrate of 0.5 μL min−1. Mobile phase A was 98/2/0.1% water/MeCN/formic acid (v/v/v) and phase B was MeCN with 0.1% formic acid (v/v). During a 108 min method, peptides were separated by a 3-step linear gradient (5% to 30% B over 90 min, 30% to 35% B over 10 min, 35% to 95% B over 4 min) followed by a 4 min isocratic flush at 95% for 4 min before washing and a return to low organic conditions. Experiments were run as data-dependent acquisitions with ion mobility activated in PASEF mode. MS and MS/MS spectra were collected with m/z 100 to 1700 and ions with z = +1 were excluded.

Raw data files were searched using PEAKS Online Xpro 1.6 (Bioinformatics Solutions Inc.). The precursor mass error tolerance and fragment mass error tolerance were set to 20 ppm and 0.03 respectively. The trypsin digest mode was set to semi-specific and missed cleavages was set to 2. The I. scapularis reference proteome (Proteome ID UP000001555, taxon 6945) and Bb reference proteome (Proteome ID UP000001807, strain ATCC 35210/B31) was downloaded from Uniprot, totaling 21,774 entries. The I. scapularis proteome was the primary search reference and the Bb was used as a secondary to identify any bacterial proteins present. Carbamidomethylation was selected as a fixed modification. Oxidation (M) and Deamidation (NQ) was selected as a variable modification. Experiments were performed in biological triplicate, with samples being a single run on the instrument. Proteins present in a database search (−10 log(p-value) ≥ 20, 1% peptide and protein FDR) were subjected to the following filtration process. We first filtered to include proteins found in 2 out of 3 biological replicates, within each respective day (one or four days after attachment). We calculated the mean area of proteins found in uninfected and infected samples. Proteins with missing values (i.e. not identified in a sample) were set to 1. We calculated the ratio of mean area for each protein as infected/uninfected, and enriched proteins were identified by having a infected:uninfected ratio greater than 2 within their respective feeding day (one or four days after attachment).

To classify the proteins we identified into functional groups, we used a PANTHER Overrepresentation Test (released 07/12/2022) (Mi et al., 2019) with PANTHER version 17.0 (released 02/22/2022) (Thomas et al., 2022). We used all Ixodes scapularis genes in the database as our reference list and analyzed the lists of proteins enriched on each day for their PANTHER Protein Class. The test type used was Fisher’s Exact, calculating a false discovery rate as the correction.

Raw data files and searched datasets are available on the Mass Spectrometry Interactive Virtual Environment (MassIVE), a full member of the Proteome Xchange consortium under the identifier: MSV000090560.

Supplemental tables

Table S1. Overview of mapping statistics from 16 Bb sequencing samples.

Table S2. Differential expression analysis results between in vitro cultured Bb before and after Bb enrichment.

Table S3. Transcriptome-wide differential expression analysis results from Bb across tick feeding timepoints.

Table S4. Twofold differentially expressed Bb genes from across tick feeding timepoints.

Table S5. Transcripts per million (TPM) for all Bb genes across feeding timepoints.

Table S6. Mass spectrometry analysis for Bb-enriched samples.

Table S7. Annotation and GO term enrichment for tick proteins enriched on feeding day 1.

Table S8. Annotation and GO term enrichment for tick proteins enriched on feeding day 4.

Acknowledgements

We are grateful to all members of the Chou lab for their feedback throughout the project and specifically to Fauna Yarza and Patrick Rockefeller Grimes for assistance with tick feeding and Ethel Enoex-Godonoo for administrative assistance. We thank Amy Lyden and Emily Crawford for assistance with and reagents for DASH along with Olga Botvinnik and The Chan Zuckerberg Biohub sequencing team for their input and sequencing assistance. We thank Patricia Rosa, Jenny Wachter, Scott Samuels, and Meghan Lybecker for their feedback on the project. We also thank William Hatleberg for assistance with figure schematics. This work was funded by: a Life Sciences Research Foundation fellowship from the SVCF-Wave Fund to Anne Sapiro, a Beckman Young Investigator award from the Arnold and Mabel Beckman Foundation to Balyn Zaro, grants from CZ Biohub, the Pew Biomedical Research Foundation, and NIH funding 1R01AI132851 to Seemay Chou, and NIH INBRE funding P20GM103474 to Patrick Secor and Margie Kinnersley.

Competing interests

Seemay Chou is president and CEO of Arcadia Biosciences.