Karasuyama H, Tabakawa Y, Ohta T, Wada T and Yoshikawa S (2018) Crucial Role for Basophils in Acquired Protective Immunity to Tick Infestation. Front. Physiol. 9:1769. doi: 10.3389/fphys.2018.01769 Ticks are blood-sucking arthropods that can transmit various pathogenic organisms to host animals and humans, causing serious infectious diseases including Lyme disease. Tick feeding induces innate and acquired immune responses in host animals, depending on the combination of different species of animals and ticks. Acquired tick resistance (ATR) can diminish the chance of pathogen transmission from infected ticks to the host. Hence, the elucidation of cellular and molecular mechanism underlying ATR is important for the development of efficient anti-tick vaccines. In this review article, we briefly overview the history of studies on ATR and summarize recent findings, particularly focusing on the role for basophils in the manifestation of ATR. In several animal species, including cattle, guinea pigs, rabbits and mice, basophil accumulation is observed at the tick re-infestation site, even though the frequency of basophils among cellular infiltrates varies in different animal species, ranging from approximately 3% in mice to 70% in guinea pigs. Skin-resident, memory CD4+ T cells contribute to the recruitment of basophils to the tick re-infestation site through production of IL-3 in mice. Depletion of basophils before the tick re-infestation abolishes ATR in guinea pigs infested with Amblyomma americanum and mice infested with Haemaphysalis longicornis, demonstrating the crucial role of basophils in the manifestation of ATR. The activation of basophils via IgE and its receptor FcεRI is essential for ATR in mice. Histamine released from activated basophils functions as an important effector molecule in murine ATR, probably through promotion of epidermal hyperplasia which interferes with tick attachment or blood feeding in the skin. Accumulating evidence suggests the following scenario. The 1st tick infestation triggers the production of IgE against tick saliva antigens in the host, and blood-circulating basophils bind such IgE on the cell surface via FcεRI. In the 2nd infestation, IgE-armed basophils are recruited to tickfeeding sites and stimulated by tick saliva antigens to release histamine that promotes epidermal hyperplasia, contributing to ATR. Further studies are needed to clarify whether this scenario in mice can be applied to ATR in other animal species and humans.

Keywords: basophil, mast cell, tick resistance, IgE, histamine

## INTRODUCTION

fphys-09-01769 December 5, 2018 Time: 12:34 # 2

Ticks, particularly ixodid family members, are blood-sucking ectoparasites of vertebrates and can transmit various pathogens to animals and humans during blood feeding for days, causing serious infectious diseases, including Lyme disease, babesiosis, Rocky Mountain spotted fever, human monocytic ehrlichiosis and severe fever with thrombocytopenia syndrome (Gratz, 1999; Parola and Raoult, 2001; de la Fuente et al., 2008; Embers and Narasimhan, 2013; Wikel, 2013; Yamaji et al., 2018). Besides tickborne infectious diseases, some people with the experience of tick bites show recurrent episodes of anaphylaxis, a life-threatening systemic allergic reaction, after eating red meat or treating with anticancer monoclonal antibodies (Platts-Mills and Commins, 2013; Steinke et al., 2015). Thus, tick infestation is of medical and veterinary public health importance.

Host defense mechanism is a threat to successful blood feeding by ticks and hence must be counteracted. To this end, ticks inject saliva containing various bioactive substances into the host during tick infestation, including vasodilator and antihemostatic, antiinflammatory and immunosuppressive reagents (Wikel, 2013). On the other hand, some animals, such as mice, guinea pigs, rabbits and cattle, have been shown to develop the resistance to tick feeding after single or multiple tick infestations, depending on the combination of animal species/strains and tick species (Trager, 1939; Wikel, 1996). This acquired tick resistance (ATR) is commonly assessed by several parameters, including the reduction in the number and/or body weight of engorged ticks or tick death when sensitized animals are re-infested with ticks. ATR was first described in 1938 by Trager who found that after infestation with Dermacentor variabilis, guinea pigs develop resistance to subsequent tick infestations (Trager, 1939). Since then, ATR has been further characterized by using cattle and laboratory animals including guinea pigs (Wikel, 1996). ATR is not restricted to the skin lesion of previous tick bites and can be observed in un-infested skin of sensitized animals, indicating the contribution of systemic responses rather than a localized response at the previously infested skin lesion. Moreover, ATR can be transferred to naive animals with sera or cells isolated from previously infested animals (Wikel and Allen, 1976; Brown and Askenase, 1981; Askenase et al., 1982), suggesting that ATR is a type of immune reaction. Importantly, ATR can diminish the chance of pathogen transmission from infected ticks to host animals and humans (Bell et al., 1979; Wikel et al., 1997; Nazario et al., 1998; Burke et al., 2005; Dai et al., 2009). Therefore, the elucidation of cellular and molecular mechanisms underlying ATR is important for developing efficient antitick vaccines that can minimize the transmission of pathogens causing serious infectious diseases.

Basophils are the least abundant granulocytes and account for less than 1% of peripheral blood leukocytes (Galli, 2000). They are named after basophilic granules in the cytoplasm that stain with basic dye, as first documented by Paul Ehrlich in 1879. In addition to basophilic granules, blood-circulating basophils share some phenotypic properties with tissue-resident mast cells, such as the surface expression of the high-affinity IgE receptor FcεRI and the release of allergy-inducing chemical mediators, including histamine, in response to various stimuli (Galli, 2000; Stone et al., 2010). Therefore, basophils have often been erroneously considered as minor and redundant relatives or blood-circulating precursors of tissue-resident mast cells (Falcone et al., 2000). It is now accepted well that basophils and mast cells are distinct cell lineages, and that basophils play crucial and non-redundant roles distinct from those played by mast cells (Voehringer, 2017; Karasuyama et al., 2018; Varricchi et al., 2018). Basophils contribute to protective immunity, particularly to parasitic infections while they are involved in the pathogenesis of various disorders, including allergic and autoimmune disorders.

In this review article, we summarize recent advances in our understanding of the cellular and molecular mechanisms underlying ATR, particularly focusing on the role of basophils identified mainly in mouse models of tick infestation.

### BASOPHILS ARE KEY EFFECTOR CELLS IN THE MANIFESTATION OF ATR

#### In Guinea Pigs

An early study described that cutaneous reactions at tickfeeding sites in tick-resistant guinea pigs were characterized by granulocytic inflammatory infiltrates, edema, and epidermal hyperplasia whereas the 1st tick-feeding site in previously uninfested guinea pigs showed minimal skin reactivity (Trager, 1939). Accumulation of numerous basophils and eosinophils, with basophils comprising up to 70% of cellular infiltrates, was detected at tick-feeding sites of guinea pigs that manifested ATR (Allen, 1973). Such basophil-rich cutaneous reaction was referred as cutaneous basophil hypersensitivity (CBH) and extensively studied in 1970s and early 1980s (Katz, 1978). Basophil depletion in A. americanum-infested guinea pigs by using antiserum raised against basophils abolished ATR (Brown et al., 1982), demonstrating the important role for basophils in ATR. Basophil infiltration at the site of tick re-infestation was also observed in cattle and rabbits (Allen et al., 1977; Brossard and Fivaz, 1982), even though the frequency of basophils among cellular infiltrates varied, and the functional role of basophils in these animals has not yet been determined to our knowledge. Thus, it remained elusive whether the important finding on basophils in guinea pig ATR can be generalized to other animal species and humans.

#### In Mice

A previous study reported that basophil infiltration was hardly detected at the tick-feeding site of WBB6F1-+/+ mice during re-infestation with H. longicornis, in spite of the fact that the mice showed ATR (Matsuda et al., 1990). Mast cell-deficient WBB6F1-W/Wv mice failed to manifest ATR, and adoptive transfer of mast cells conferred ATR on these mice (Matsuda et al., 1985, 1987, 1990), suggesting that mast cells in place of basophils contributed to ATR in mice, unlike in guinea pigs. On the contrary, other studies reported that the same mast celldeficient strain of mice showed ATR to another tick species Dermacentor variabilis (denHollander and Allen, 1985; Steeves and Allen, 1991). Murine basophils had been notoriously difficult to identify owing to their fewer basophilic granules compared

to those in other animals and humans, and therefore, electron microscopic examination was needed to identify them in tissue sections (Urbina et al., 1981; Dvorak et al., 1982; Dvorak, 2000). Notably, the infiltration of basophils, along with eosinophils and neutrophils, was detected by electron microscopy at the tickfeeding site in the 3rd infestation with D. variabilis in both mast cell-sufficient and -deficient mice (Steeves and Allen, 1991). Thus, the mechanism underlying ATR in mice, including the distinct roles played by basophils and mast cells, and the influence of different genetic background of both mice and ticks remained to be clarified.

Recent characterization of cell surface markers on murine basophils (Min et al., 2004; Voehringer et al., 2004) and the identification of murine basophil-specific serine protease, mouse mast cell protease-8 (mMCP-8) (Poorafshar et al., 2000; Ugajin et al., 2009) have enabled us to identify and isolate murine basophils much more easily. Taking the advantage of a mMCP-8-specific mAb TUG8 (Ugajin et al., 2009), we demonstrated that mMCP-8-expressing basophils are recruited to the tick-feeding site and make a cluster around the tick mouthpart during the 2 nd but rarely the 1st infestation with H. longicornis in C57BL/6 mice (Wada et al., 2010). Intravital fluorescence microscopic analysis, using Mcpt8GFP (green basophil) mice in that only basophils express green fluorescent protein (GFP), confirmed the basophil accumulation at the 2nd but not 1st tick-feeding site (Ohta et al., 2017; **Figure 1**). Basophils represented less than 5% of leukocytes at the 2nd tick-feeding site in mice, much fewer than in guinea pigs, while monocytes/macrophages, neutrophils and eosinophils were abundant. Importantly, we found that basophil depletion by treating mice with basophil-depleting mAbs, either anti-FcεRIα (MAR-1) or anti-CD200R3 (Ba103), just before the 2nd tick infestation completely abolished ATR with no apparent effect on the number of other types of cells, including monocytes/macrophages, neutrophils and eosinophils (Wada et al., 2010). The essential role of basophils in ATR was further demonstrated by diphtheria toxin-mediated ablation of basophils in genetically engineered Mcpt8DTR mice in that only basophils expressed diphtheria toxin receptors (Wada et al., 2010). Of note, we also demonstrated that mast cell-deficient KitW−sh/W−sh C57BL/6 mice failed to manifest ATR, confirming the importance of mast cells in ATR reported previously (Matsuda et al., 1985, 1987, 1990). Thus, mast cells, in addition to basophils, appear to contribute to ATR in C57BL/6 mice infested with H. longicornis, in contrast to ATR in D. variabilis-infested WBB6F1-+/+ mice in that mast cells are dispensable (denHollander and Allen, 1985; Steeves and Allen, 1991).

#### In Humans

Basophil infiltration was detected in humans at the tick-feeding sites and in the skin lesions of scabies (Ito et al., 2011; Nakahigashi et al., 2013; Kimura et al., 2017). Of note, a patient lacking basophils and eosinophils reportedly suffered from widespread scabies (Juhlin and Michaelsson, 1977). These observations suggest the possible involvement of basophils in protective immunity to ectoparasites, including ticks.

#### BASOPHIL ACTIVATION THROUGH IgE AND ITS RECEPTOR FcεRI IS ESSENTIAL FOR ATR

It was shown in guinea pigs that transfer of serum from previously infested animals conferred ATR on naive animals (Wikel and Allen, 1976; Brown and Askenase, 1981; Askenase et al., 1982). Similarly, in mice, transfer of serum from tickinfested but not un-infested mice conferred ATR on naive mice (Matsuda et al., 1990), suggesting the involvement of tick-specific antibodies in ATR. Of note, the heat treatment of the serum at 56◦C for 2 h abolished the ATR transfer activity (Matsuda et al., 1990), indicating that antibodies of IgE isotype contribute to the manifestation of ATR. Consistent with this observation, we demonstrated that both antibody-deficient µMT mice and Fcer1g−/− mice, that lack the expression of high affinity IgE receptor FcεRI, failed to show ATR (Wada et al., 2010). This

suggested the following scenario (**Figure 2**). The 1st infestation triggers the production of IgE against tick saliva antigens, and basophils and mast cells bind IgE on the cell surface via FcεRI. In the 2nd infestation, tick saliva antigens delivered into the tickfeeding site bind to IgE on these cells, leading to the cross-linking of FcεRI and hence activation of these cells that may contribute to ATR.

Intriguingly, mast cell-deficient KitW−sh/W−sh C57BL/6 mice reconstituted with mast cells derived from Fcer1g−/− mice could manifest ATR as did mice reconstituted with wild-type mast cells (Wada et al., 2010), indicating that FcεRI on mast cells is dispensable for IgE-mediated ATR. In contrast, adoptive transfer of basophils isolated from previously infested wild-type, but not Fcer1g−/−, mice conferred ATR on naive mice (Wada et al., 2010). These results suggested that basophils rather than mast cells play a critical role in IgE-dependent ATR through FcεRImediated activation, even though both types of cells contribute to ATR.

Ticks inject a plethora of substances, including proteins, into the host during feeding (Wikel, 2013). However, it remains ill-defined which components among tick saliva injected are the major targets of IgE that is involved in ATR, even though a series of tick saliva antigens recognized by sera from tickinfested animals and humans have been identified (Brown et al., 1984; Brown, 1988; Mayoral et al., 2004). It was demonstrated that infestation with A. americanum can induce a strong IgE response to tick saliva antigens including the carbohydrate α-gal, which is also present in red meats such as beef and pork. The production of such anti-α-gal IgE in tick-infested people can lead to anaphylaxis after ingestion of red meats (Platts-Mills and Commins, 2013; Steinke et al., 2015). It remains to be investigated whether anti-α-gal IgE is involved in ATR.

Molecular characterization of tick salivary components has demonstrated that different members among the same multigene family are expressed at distinct time points during tick feeding (Karim and Ribeiro, 2015). For example, two cystatin genes from Ixodes scapularis change their expression reciprocally during feeding (Kotsyfakis et al., 2007; Karim and Ribeiro, 2015). Such antigenic variation or sialome switch during tick feeding is considered as a possible mechanism by which ticks avoid host

immune responses. It remains to be determined whether such variation can affect the production of anti-tick IgE and hence IgE-mediated ATR, and whether IgE raised against one family member is cross-reactive to other members of the same family.

Host-derived IgG molecules containing blood meal pass through the midgut barrier of Rhipicephalus appendiculatus into the hemolymph and are excreted via the saliva back into the host during feeding. IgG binding proteins detected in the tick hemolymph and salivary glands are thought to contribute to this excretion of IgG, as a strategy by which ticks evade the damage caused by host antibodies (Wang and Nuttall, 1999). IGBP-MA, a member of IgG binding proteins has been shown to bind to IgE (Wang and Nuttall, 2013). Further studies are needed to examine whether such IgG binding proteins can interfere with IgE-mediated ATR in the host and whether the host raises antibodies against them to neutralize their activity.

#### BASOPHIL-DERIVED HISTAMINE IS AN IMPORTANT EFFECTOR MOLECULE IN ATR

Biologically active molecules, such as histamine and proteases, stored in the secretary granules in basophils and mast cells have been implicated as effectors of ATR. It was reported in cattle that the tick resistance is correlated with hypersensitivity to tick antigens and the amount of histamine at the tickfeeding site (Willadsen et al., 1979). Moreover, administration of antihistamine in cattle resulted in higher tick numbers (Tatchell and Bennett, 1969) whereas the injection of histamine into the cattle skin promoted tick detachment (Kemp and Bourne, 1980). Similar observations were reported in guinea pigs (Wikel, 1982), suggesting the possible involvement of histamine to ATR. However, the cellular source of histamine responsible for ATR and the mechanism underlying histamine-mediated ATR remained ill-defined.

We have recently addressed these questions by analyzing C57BL/6 mice infested with H. longicornis, in that both basophils and mast cells contribute to ATR (Wada et al., 2010). Treatment of mice with histamine H1 antagonist during the 2nd infestation abolished ATR (Tabakawa et al., 2018). Consistent with this, mice deficient for histamine production due to the lack of histidine decarboxylase (HDC) failed to show ATR (Tabakawa et al., 2018). Moreover, repeated injection of histamine or histamine H1 receptor agonist beneath the tick-infested site during the 1st infestation inhibited the tick feeding in wild-type mice (Tabakawa et al., 2018). These observations illustrated the important role of the histamine-histamine H1 receptor axis in the manifestation of ATR in mice, consistent with previous studies in guinea pigs and cattle (Tatchell and Bennett, 1969; Willadsen et al., 1979; Wikel, 1982).

Both basophils and mast cells are well-known producers of histamine, and therefore supposed to contribute to ATR through histamine release. Unexpectedly, however, adoptive transfer of histamine-deficient mast cells reconstituted ATR in mast cell-deficient KitW−sh/W−sh C57BL/6 mice as did that of wild-type mast cells (Tabakawa et al., 2018), indicating that mast cell-derived histamine is dispensable for ATR. In contrast, adoptive transfer of wild-type but not histamine-deficient basophils conferred ATR on basophil-depleted Mcpt8DTR mice (Tabakawa et al., 2018), demonstrating the crucial role of basophil-derived histamine in the manifestation of ATR.

Intravital imaging analysis of cells at the 2nd tick feeding site demonstrated that basophils make a cluster within the epidermis and surround a tick mouthpart. In contrast, mast cells are mostly scattered in the dermis rather than epidermis and localized more distantly from the tick mouthpart (Tabakawa et al., 2018). It is well known that histamine has a short half-life. Therefore, basophil-derived histamine may be much more effective than mast cell-derived in the manifestation of ATR, considering the fact that higher numbers of basophils are localized closer to a tick mouthpart, compared to mast cells.

Previous studies reported that histamine promotes itching and grooming response in the skin, resulting in removal of ticks in host animals (Koudstaal et al., 1978). In the mouse mode of tick infestation, ticks are placed inside of a small tube attached to the skin. Therefore, the effect of host grooming on tick feeding is minimized, implying other mechanisms underlying histaminemediated ATR. Mice deficient for histamine H1 receptor failed to manifest ATR (Tabakawa et al., 2018), indicating that histamine acts on host cells rather than ticks. We detected the thickening of the epidermis and the formation of basophil cluster within the thickened epidermis at the 2nd but not 1st tick-feeding site in mice (Tabakawa et al., 2018) as reported previously in guinea pigs (Trager, 1939; Allen, 1973). This epidermal hyperplasia was absent in histamine-deficient or basophil-deficient mice (Tabakawa et al., 2018), suggesting that basophil-derived histamine is involved in epidermal hyperplasia. Considering that keratinocytes express functional H1 receptor (Ohsawa and Hirasawa, 2014) and that histamine promotes the proliferation of keratinocytes (Maurer et al., 1997; Albrecht and Dittrich, 2015), histamine released from basophil localized in the epidermis perhaps induces the thickening of the epidermis that may interfere with tick attachment or blood-sucking in the skin during the 2nd infestation (**Figure 2**).

Histamine-binding proteins (HBPs) have been identified in tick saliva (Paesen et al., 1999; Sangamnatdej et al., 2002; Mans et al., 2008). They show high-affinity binding to histamine and can efficiently compete for histamine with its native receptor. Thus, they may interfere with histamine-induced host responses at tick feeding sits, including itching and grooming. However, it remains to be determined whether tick HBPs can give any impact on histamine-mediated ATR in the host and whether the host raises antibodies against them to neutralize their activity. Mast cells and basophils are the major source of histamine at tick feeding sites. Basophils accumulate at tick feeding sites during the 1 st but not 2nd infestation while mast cells always reside there. Given that higher numbers of basophils are localized closer to a tick mouthpart, compared to mast cells, during the 2nd infestation (Tabakawa et al., 2018), the concentration of histamine near tick mouthparts should be much higher during the 2nd infestation compared to the 1st infestation. Therefore, one may assume that HBPs might be less effective in sequestering histamine at the 2nd tick-feeding site in which ATR is executed. The influence of HBPs

on histamine-mediated ATR could be explored by generating HBP-deficient ticks in future studies.

It has been reported that H. longicornis, Dermacentor andersoni, and Boophilus microplus larval ticks are highly reactive to histamine in the induction of tick resistance while A. americanum and Ixodes holocyclus ticks are less responsive to histamine (Bagnall, 1975; Kemp and Bourne, 1980; Wikel, 1982; Brown and Askenase, 1985). The former tick species have shorter mouthparts than the latter (Suppan et al., 2017), suggesting the possibility that histamine-induced thickening of the epidermis prevents the former's but not the latter's mouthparts from penetrating into the dermis in order to form blood pools. This may explain the differential responsiveness to histamine among tick species in terms of ATR induction. Alternatively, but not mutually exclusively, it is possible that the presence or absence (or differential amounts) of HBPs in different tick species is correlated in part with differential reactivity to histamine in the induction of tick resistance.

#### SKIN-RESIDENT MEMORY CD4+ T CELLS ARE RESPONSIBLE FOR BASOPHIL RECRUITMENT TO THE 2nd TICK-FEEDING SITE

Basophils circulate in the peripheral blood under homeostatic conditions, and they infiltrate the skin at the tick-feeding site during the 2nd but not 1st infestation. Importantly, the recruitment of basophils can be observed in previously uninfested skin, far from the 1st infestation site, of sensitized animals, implicating that the 1st tick infestation may induce systemic alteration in the skin throughout the body, so that basophils can readily infiltrate the tick re-infestation site anywhere in the body. We have recently demonstrated in mice that skin CD4+ memory T cells play an important role in basophil recruitment to the 2nd tick-feeding site, leading to ATR (Ohta et al., 2017). Tick antigen-specific CD4+ effector T cells are generated during the 1st tick infestation and distributed to the skin all over the body, and some of them are retained as skin-resident memory T cells (**Figure 2**). In the 2nd tick infestation, tick saliva antigens delivered into the skin stimulate these memory T cells present in the skin to produce IL-3 that is required for basophil recruitment to the 2nd tick-feeding site (Ohta et al., 2017). Even though the exact mechanism underlying IL-3-mediated basophil recruitment remains to be clarified, IL-3 might promote basophil adhesion to endothelium (Bochner et al., 1990; Korpelainen et al., 1996; Lim et al., 2006), leading to transendothelial migration of basophils and their accumulation in the skin.

#### AN UNSOLVED ISSUE: THE ROLE OF MAST CELLS IN ATR

As described above, mast cells contribute to ATR in mice infested with H. longicornis (Matsuda et al., 1985, 1987, 1990) whereas they are dispensable for ATR in D. variabilis-infested mice (denHollander and Allen, 1985; Steeves and Allen, 1991). As far as we are aware, the involvement of mast cells to ATR has not yet been documented in other animal species. In the case of mice infested with H. longicornis, histamine derived from basophils but not mast cells is essential for the manifestation of ATR (Tabakawa et al., 2018), even though both basophils and mast cells are involved in ATR (Wada et al., 2010). The deficiency of either basophils or mast cells almost completely abolishes ATR (Wada et al., 2010), suggesting that the role of these cells may not be additive. Of note, the number of basophils accumulating at the 2nd tick-feeding site is comparable between mast cellsufficient and -deficient mice (Wada et al., 2010), indicating that mast cells are not prerequisite for basophil recruitment. Nevertheless, closer examination with intravital imaging revealed that basophils accumulating at the 2nd tick-feeding site are more motile and less-clustered around a tick mouthpart in mast cell-deficient mice than in mast cell-sufficient mice (Tabakawa et al., 2018). Therefore, one may assume that mast cells may contribute to ATR by directly or indirectly regulating basophil behavior. Further studies are needed for elucidating how mast cells contribute to the manifestation of ATR.

## CONCLUSION

Recent development of a series of analytical tools in laboratory animals has advanced our understanding of the cellular and molecular mechanism underlying ATR. In several animal species, basophil accumulation is observed at the tick re-infestation site (**Figure 1**), and basophil depletion abolishes ATR in guinea pigs and mice, demonstrating the crucial role of basophils in the manifestation of ATR. The 1st tick infestation triggers the production of IgE against tick saliva antigens. In the 2 nd infestation, IgE-armed basophils are recruited to the tickfeeding site and stimulated by tick saliva antigens to release histamine that functions as a key effector in ATR, probably through promotion of the epidermal hyperplasia that in turn interferes with tick attachment or blood feeding in the skin (**Figure 2**). Further studies on the detailed mechanism underlying ATR, including the role of mast cells, may help develop the strategy to prevent tick infestation and tick-borne diseases.

#### AUTHOR CONTRIBUTIONS

HK has substantially contributed to the literature review and drafted the manuscript. SY contributed to the discussion and drafting, and editing of the manuscript. YT, TO, and TW reviewed and edited the manuscript.

## FUNDING

This work was supported by research grants from Japanese Ministry of Education, Culture, Sports, Science and Technology [15H05786 (HK) and 17K15719 (SY)] and TMDU President's Young Researchers Award (SY).

Frontiers in Physiology | www.frontiersin.org

#### REFERENCES

fphys-09-01769 December 5, 2018 Time: 12:34 # 7

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Karasuyama, Tabakawa, Ohta, Wada and Yoshikawa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Transcriptome and Proteome Response of Rhipicephalus annulatus Tick Vector to Babesia bigemina Infection

Sandra Antunes1 \*, Joana Couto1 , Joana Ferrolho1 , Gustavo Seron Sanches1 , José Octavio Merino Charrez2 , Ned De la Cruz Hernández2 , Monica Mazuz3 , Margarita Villar4 , Varda Shkap3 , José de la Fuente4,5 and Ana Domingos1

1 Global Health and Tropical Medicine, Instituto de Higiene e Medicina Tropical, Universidade Nova de Lisboa, Lisbon, Portugal, 2 Facultad de Medicina Veterinaria y Zootecnia, Universidad Autónoma de Tamaulipas, Ciudad Victoria, Mexico, 3 Kimron Veterinary Institute, Bet Dagan, Israel, 4 SaBio, Instituto de Investigación en Recursos Cinegéticos, IREC, CSIC-UCLM-JCCM, Ciudad Real, Spain, 5 Department of Veterinary Pathobiology, Center for Veterinary Health Sciences, Oklahoma State University, Stillwater, OK, United States

#### Edited by:

Itabajara Da Silva Vaz Jr., Federal University of Rio Grande do Sul, Brazil

#### Reviewed by:

Mauro Mandrioli, University of Modena and Reggio Emilia, Italy Claudio Mafra, Universidade Federal de Viçosa, Brazil Tae Kim, Texas A&M University, United States

> \*Correspondence: Sandra Antunes<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Invertebrate Physiology, a section of the journal Frontiers in Physiology

Received: 10 January 2019 Accepted: 11 March 2019 Published: 02 April 2019

#### Citation:

Antunes S, Couto J, Ferrolho J, Sanches GS, Merino Charrez JO, De la Cruz Hernández N, Mazuz M, Villar M, Shkap V, de la Fuente J and Domingos A (2019) Transcriptome and Proteome Response of Rhipicephalus annulatus Tick Vector to Babesia bigemina Infection. Front. Physiol. 10:318. doi: 10.3389/fphys.2019.00318 A system biology approach was used to gain insight into tick biology and interactions between vector and pathogen. Rhipicephalus annulatus is one of the main vectors of Babesia bigemina which has a massive impact on animal health. It is vital to obtain more information about this relationship, to better understand tick and pathogen biology, pathogen transmission dynamics, and new potential control approaches. In ticks, salivary glands (SGs) play a key role during pathogen infection and transmission. RNA sequencing obtained from uninfected and B. bigemina infected SGs obtained from fed female ticks resulted in 6823 and 6475 unigenes, respectively. From these, 360 unigenes were found to be differentially expressed (p < 0.05). Reversed phase liquid chromatography–mass spectrometry identified a total of 3679 tick proteins. Among them 406 were differently represented in response to Babesia infection. The omics data obtained suggested that Babesia infection lead to a reduction in the levels of mRNA and proteins (n = 237 transcripts, n = 212 proteins) when compared to uninfected controls. Integrated transcriptomics and proteomics datasets suggested a key role for stress response and apoptosis pathways in response to infection. Thus, six genes coding for GP80, death-associated protein kinase (DAPK-1), bax inhibitor-1 related (BI-1), heat shock protein (HSP), heat shock transcription factor (PHSTF), and queuine trnaribosyltransferase (QtRibosyl) were selected and RNA interference (RNAi) performed. Gene silencing was obtained for all genes except phstf. Knockdown of gp80, dapk-1, and bi-1 led to a significant increase in Babesia infection levels while hsp and QtRibosyl knockdown resulted in a non-significant decrease of infection levels when compared to the respective controls. Gene knockdown did not affect tick survival, but engorged female weight and egg production were affected in the gp80, dapk-1, and QtRibosylsilenced groups in comparison to controls. These results advanced our understanding of tick–Babesia molecular interactions, and suggested new tick antigens as putative targets for vaccination to control tick infestations and pathogen infection/transmission.

Keywords: proteomics, transcriptomics, Rhipicephalus annulatus, Babesia bigemina, vector–pathogen interactions, apoptosis, stress response

**52**

## INTRODUCTION

fphys-10-00318 March 30, 2019 Time: 16:7 # 2

Babesiosis is a worldwide tick-borne hemoprotozoa disease caused by intra-erythrocytic parasites of the genus Babesia (Beugnet and Moreau, 2015; Lempereur et al., 2017). The disease may range from asymptomatic carrier to more severe states, characterized by hemolytic anemia, fever, hemoglobinuria, and occasionally death, affecting a large variety of mammals, including pets, farm, and wild animals and also humans (Schnittger et al., 2012). Human infection is mainly due to Babesia microti or Babesia divergens, and although it rarely occurs, it has been considered an emerging zoonosis due to the growing number of fatal cases (Leiby, 2006; Schnittger et al., 2012). In contrast, cattle babesiosis, caused by either Babesia bovis or B. bigemina, is an important disease causing high morbidity and mortality, thus leading to severe economic losses to the cattle industry (Chauvin et al., 2009). The tick species R. annulatus and Rhipicephalus microplus, the most important ticks of cattle in tropical and subtropical regions, conduce a negative impact on meat, milk, and leather productions and are considered as the main vectors and reservoirs of B. bovis and B. bigemina (Bock et al., 2004).

The application of chemical acaricides is the method of choice for tick control. However, it results in environmental contamination, selection of resistant ticks, and presence of residues in meat and milk, potentially harmful for animals and humans (Corson et al., 2001; Ghosh et al., 2007). To reduce these negative impacts, much attention has been directed to the development of new approaches that are efficient and at the same time environmentally friendly. Vaccination, as a prophylactic measure, stands out representing a promising and sustainable alternative for the control of ticks and tick-borne pathogens (Almazan et al., 2018). A vaccine targeting both tick fitness and pathogen competence is an attractive choice requiring the identification of tick molecules with a dual effect.

Several studies, based on omics approaches, have been conducted to understand the tick–pathogen interactions, identifying possible new vaccine candidates (Villar et al., 2017). System biology approaches have been efficiently used to characterize vector and pathogen interactions: in Drosophila melanogaster, one of the most well studied genetic model organisms (McTaggart et al., 2015), and Anopheles spp. mosquitoes (Domingos et al., 2017). The characterization of tick organs response to infection using technologies such as transcriptomics, proteomics, and functional genomics improved current knowledge on tick–pathogen interactions and allowed the development of new strategies and/or the identification of targets for tick and disease control. Also, research has shown that during the long-lasting tick–pathogen co-evolution, microorganisms have developed important strategies to manipulate or modulate tick response to infection, without impairing tick survival, enhancing their capacity of infection, replication, and transmission guaranteeing the survival of both (de la Fuente et al., 2016; Šimo et al., 2017). In Ixodes scapularis, the effect of Anaplasma phagocytophilum infection on tissue-specific responses and on cellular pathways, such as apoptosis or stress response, which can be activated in a certain tissue in response to infection in order to contain parasite evasion and or improve immune response, have been extensively studied (Alberdi et al., 2015; Ayllon et al., 2015; Villar et al., 2015; Cabezas-Cruz et al., 2017). However, there is limited knowledge on tick vector responses to Babesia spp. infection. Babesia parasites invade several tick tissues including the midgut, salivary glands (SGs), and ovaries, affecting tick fitness. Some tick molecules such as calreticulin, a calcium binding protein, kunitztype serine protease inhibitors, lachesin, vitellogenins, among others, have been identified as having a role on tick–Babesia interface (Zhou et al., 2006; Rachinsky et al., 2008; Antunes et al., 2017, 2018). Scarce information is available regarding the mechanisms used by B. bigemina to infect, develop, multiply, and survive inside the tick vector. Additionally, the impact of parasite infection at tick transcriptome and proteome levels, particularly on the molecular pathways affected by B. bigemina, is still to be investigated.

The overall objective of this study was to deepen the understanding on the complex R. annulatus response to B. bigemina infection. To conclude on mRNA and protein levels of R. annulatus in response to B. bigemina infection, the present research combined transcriptomics and proteomics analysis to obtain a sialotranscriptome and proteome by RNA sequencing (RNA-seq) and reversed phase liquid chromatography–mass spectrometry (RP-LC–MS/MS). Research focused on genes and proteins that were found differentially expressed or represented after infection. Six genes were further studied by RNAi-mediated gene silencing including genes related to apoptosis and stress response. This work represents the first report concerning the effect of B. bigemina infection on the sialotranscriptome and proteome of R. annulatus, and in the influence of the presence of the parasite on specific and crucial cellular pathways, constituting an important step further on the development of new measures for ticks and parasite control.

## MATERIALS AND METHODS

## Ethics Statement

Animal experiments were conducted according with the "Guide for Care and Use of Laboratory Animals" of the institutions involved in the study, following protocols approved by the Committee on the Ethics of Animal Experiments and the principle of the three Rs, to replace, reduce, and refine the use of animals for scientific purposes.

## Rhipicephalus annulatus and Babesia bigemina

R. annulatus adult ticks were obtained from a laboratory colony free of tick-transmissible infections maintained at the Kimron Veterinary Institute, Israel. Six-month-old Friesian calves were tested for the presence of antibodies against Babesia spp. using an immunofluorescence assay, as described previously (Shkap et al., 2005), and were kept under strict tick-free conditions. To obtain B. bigemina-infected ticks, one calf was splenectomized and inoculated intravenously with 1 × 106 B. bigemina (Moledet strain) cryopreserved parasites a fortnight later. Once the peak of

parasitemia was reached, ticks were fed for 1 week inside cotton bags attached with non-toxic silicone glue GE Advantage Silicone Sealant (General Electric, New York, NY, United States) to a shaved area on the dorsal region of the animal. Similarly, a naïve calf was used to obtain uninfected ticks. Detached adult female ticks from the infected and naïve calf were maintained at 28◦C and 80% humidity.

#### Sialome

#### Sample Preparation and Probe-Based qPCR Detection of Babesia bigemina

A total of 30 uninfected and 30 B. bigemina-infected female ticks were used for tissue isolation. Ticks were washed twice in distilled water, once in 75% (v/v) ethanol, and one last wash in water. Ticks were dissected and the SGs extracted under a Motic SMZ-171B stereomicroscope (Motic Instruments Inc., Xiamen, China), using ice-cold phosphate buffer saline (PBS). Total RNA from 10 SGs obtained from each group was extracted using TriReagent (Sigma–Aldrich, St. Louis, MO, United States), according to the manufacturer's instructions. RNA concentration and integrity were evaluated by analysis of rRNA band integrity using Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, United States). cDNA was synthesized from 100 ng of total RNA using the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, United States), following the manufacturer's protocol. To determine the presence of B. bigemina in the SGs, probe-based qPCR reactions were carried out as described before for the detection of 18S rRNA of B. bigemina (Kim et al., 2007). Briefly, triplicate 20 µl reactions were prepared with 10 µl of Probe Xpert Fast Probe 2× Mastermix (GRISP Research Solutions, Porto, Portugal), 400 nM of each primer, 100 nM of probe, 1 µl of cDNA template, and nuclease-free water up to the final volume. The qPCR was carried out in a CFX ConnectTM Real-Time PCR Detection System thermal cycler (Bio-Rad, Hercules, CA, United States) with a thermal cycling profile of 10 min at 95◦C, followed by 45 cycles of 20 s at 95◦C and 1 min at 55◦C. Negative controls were prepared with no template. Positive controls were prepared with B. bigemina Moledet strain purified DNA. A standard curve was constructed with fivefold serial dilutions of the positive control DNA to determine reaction efficiency. Data were analyzed using the Bio-Rad CFX Manager Software version 3.1 (Bio-Rad, Hercules, CA, United States). Samples with quantification cycle (Cq) values above 39 were considered negative for the presence of the pathogen. After B. bigemina infection confirmation, the remaining 20 SGs from infected or uninfected ticks were grouped in two pools, and the total RNA, for RNA-seq was extracted as described above.

#### RNA Sequencing and Analysis

Library preparation was performed using TruSeq RNA kit (Illumina, San Diego, CA, United States), following the manufacturer's instructions. Shortly, prior to cDNA library construction magnetic beads with Oligo (dT) were used to enrich poly (A) mRNA from 1 µg of total-RNA. Purified mRNA was disrupted into short fragments. Two different fragmentation conditions were applied so that a "shorter" and a "longer" preparation were made for both the control (uninfected) and the B. bigemina infected RNA samples. Next, the purified mRNA was disrupted into short fragments, and double-stranded cDNA was immediately synthesized. cDNA was subjected to end-repair and adenylation, then connected with sequencing adapters. Suitable fragments, purified by size selection protocol with AMPure XP beads (Beckman Coulter, Pasadena, CA, United States), were selected as templates for PCR amplification. The final library sizes and qualities were evaluated by electrophoresis using an Agilent High Sensitivity DNA Kit (Agilent Technologies Inc.); fragment size range was between 287 and 296 bp for the short insert and 397 and 436 bp for the longer insert preparations. Subsequently, libraries were pooled and titrated using qPCR to obtain an accurate estimation of concentration. Cluster generation was performed in a Cluster Station (Illumina, San Diego, CA, United States) and the libraries sequenced using a GAIIx equipment (Illumina, San Diego, CA, United States), with a 2 × 100 cycle sequencing reads separated by a pairedend turnaround. Image analysis was performed using the HiSeq control Software version 1.8.4. (Illumina, San Diego, CA, United States) at Parque Científico de Madrid, Spain.

#### Assembly and Analysis of Transcripts

Reads were trimmed where the error probability was higher than 0.05. To get an optimal assembly, the sequence reads included were only those that the two members of the pair remained after filtering and trimming. Oases Software (Velvet, version 1.2.10) suitable for short paired-end reads assemblies (Schulz et al., 2012) was used selecting the mode of single assembly, i.e., not merged. For the two conditions, a kmer of 83 was selected for being close to the total length of the read (115 bp), avoiding misassembles. The Transcriptome Shotgun Assembly (TSA) project was submitted to the DDBJ/EMBL/GenBank under accession numbers GBJT00000000 and GBJS00000000 for uninfected and infected ticks, respectively. Functional annotation of each transcript was conducted based on the Basic Local Alignment Search Tool (BLAST1 ) and the results inferred by similarity to UniProt database2 reference proteins. The minimum similarity threshold required for annotating a transcript was a BLAST E-value < 10−10. A total of 34,095 reference proteins, considered to be representative on the UniRef90 clusters belonging to the taxonomic node Chelicerata, were downloaded from the UniProt database (January 2015). Chelicerata is eight levels above the Rhipicentor reticulatus taxon (Madder et al., 2014). Each selected protein belonged only to one cluster, showing 90% similarity with all members of that cluster. A protein-centred analysis of the differential expression was performed for the de novo transcriptome comparison. For the UniGene cluster formation, in each sample, for a set of transcripts annotated with the same UniProt, the longer nucleotide sequence was chosen to be the representative of each UniGene. Functional data were curated using Blast2GO platform version 4.0.7 available at https://www.blast2go.com (Conesa et al., 2005; Götz et al., 2008). Manual annotation was done with sequences retrieving no hits needed by compiling information from UniProt, RefSeq, GO,

1https://blast.ncbi.nlm.nih.gov/Blast.cgi

2https://www.uniprot.org/

Panther, KEGG, Pfam, and NCBI databases. For the uninfected and infected condition comparison, transcripts were clustered by proteins. When two transcripts were annotated under the same UniProt they were included in the same cluster; and when the same protein cluster was present in both conditions, the protein cluster expression levels were compared. The software edgeR, for Empirical Analysis of Digital Gene Expression Data in R [Bioconductor version: Release (3.7)] (Robinson et al., 2010), was used to compare the mRNA expression levels detected in samples.

#### Transcriptomics Data Normalization and Validation

For transcriptomics data normalization, four reference genes were selected including tubulin-beta-2B chain, elongation factor 1, cyclophilin, and transcription factor TFIID (Robinson et al., 2010). Since no biological replicates were used, a biological coefficient of variation (BCV) was firstly selected to 0.4. Considering that if the BCV was perfectly approximated, the p-value for the four genes would have to be 1, the lowest p-value obtained for any of these four genes was selected. This corresponded to the transcription factor TFIID that had a p = 0.465 (for BCV = 0.4). All the genes with a p-value higher than 0.465 were selected as genes without differential expression to approximate the dispersion. The dispersion and the BCV were analyzed for this set of 2719 theoretically housekeeping genes. A BCV = 0.6 was obtained after approximating the dispersion of the housekeeping genes and this BCV value was used to select the differentially expressed genes in response to infection. For RNA-seq data validation, 10 paired SGs obtained individually from each uninfected and B. bigeminainfected R. annulatus female ticks (previously screened for infection), respectively, were used in triplicate. **Supplementary Table S1** shows the primer sets used for qPCR and their respective final concentrations and annealing temperatures. Briefly, 10 µl triplicate reactions were prepared with 5 µl iTaqTM Universal SYBR R Green Supermix (Bio-Rad, Hercules, CA, United States), specific primer concentration, 1 µl of cDNA template, and nuclease-free water up to the final volume. The qPCR was carried out in a CFX ConnectTM Real-Time PCR Detection System thermal cycler (Bio-Rad, Hercules, CA, United States), with a thermal cycling profile of 10 min at 95◦C, followed by 45 cycles of 15 s at 95◦C and 45 s at specific-annealing temperature. Negative controls were prepared with no template. A standard curve was constructed with fivefold serial dilutions of R. annulatus cDNA to determine reaction efficiency. To ensure that only one amplicon was formed, a melting curve was performed at the end of every reaction (55–95◦C; 0.5◦C/s melt rates). To confirm correct fragment amplification, gel electrophoresis was performed, the fragments were sequenced by the chain termination method at StabVida (Lisbon, Portugal) and the resulting sequence was analyzed against the sequences available at the NCBI nucleotide database3 . The CFX ManagerTM Software (Bio-Rad, Hercules, CA, United States) was used to analyze the gene expression data between conditions and reference gene validation was based on the geNorm algorithm (Vandesompele et al., 2002) and on the expression stability value M of each reference gene (M < 1).

# Proteome

#### Sample Preparation

A total of 16 uninfected and 16 B. bigemina-infected engorged female ticks obtained as described in the section "Rhipicephalus annulatus and Babesia bigemina" were equally divided into two different biological samples for each group. Ticks were dissected, cuticle removed, and washed in 10 mM PBS to eliminate the maximum host blood possible. Internal tick tissues were homogenized with a glass homogenizer (20 strokes) in lysis buffer [0.25 M de sucrose, 10 mM Tris–HCl pH 7.5, 1 mM MgCl2, 0.7% DDM, 0.5% ASB detergent, supplemented with Complete protease inhibitor cocktail (Roche, Basel, Switzerland)] in a ratio of 1 ml of buffer per gram of sample. Samples were sonicated for 1 min in an ultrasonic cooled bath (JP Selecta, Barcelona, Spain) followed by 10 s of vortex. After three cycles of sonication-vortex, the homogenates were centrifuged at 200 × g for 5 min at 4◦C, to remove cellular debris. The supernatants were then collected and protein concentration was determined using the Bicinchoninic Acid (BCA) Assay (Thermo Scientific, San Jose, CA, United States), according to the manufacturer's instructions.

#### Proteomics

One-hundred and fifty grams of protein extracts per group were on-gel concentrated by SDS–PAGE and trypsin digested as described previously (Villar et al., 2014). The desalted protein digest was re-suspended in 0.1% formic acid and analyzed by RP-LC–MS/MS technique using an Easy-nLC II system coupled to an ion trap LTQ mass spectrometer (Thermo Scientific, San Jose, CA, United States). Peptides were concentrated (online) by RP chromatography using a 0.1 × 20 mm C18 RP precolumn (Thermo Scientific, San Jose, CA, United States), and then separated using a 0.075 × 100 mm C18 RP column (Thermo Scientific, San Jose, CA, United States) operating at 0.3 ml/min. Peptides were eluted using a 180 min gradient from 5 to 40% solvent B (Solvent A: 0.1% formic acid in water, solvent B: 0.1% formic acid in acetonitrile). Electrospray (ESI) ionization was performed using a Fused-silica PicoTip Emitter ID 10 mm (New Objective, Woburn, MA, United States) interface. Peptides were detected in survey scans from 400 to 1600 amu (1 mscan), followed by 15 data dependent MS/MS scans (Top 15), using an isolation width of two mass-to-charge ratio units, normalized collision energy of 35%, and dynamic exclusion applied during 30 s periods.

#### Data Analysis

The MS/MS raw files were searched against a compiled database containing all sequences from Ixodidae, Bos taurus, and Babesia (135,071, 32,156, and 23,215 Uniprot entries in September 2017, respectively)4 and with a database created from transcriptomics data (Holt et al., 2002; Evans et al., 2012), according to the Proteomics Informed by Transcriptomics (PIT) technique using the SEQUEST algorithm (Proteome Discoverer 1.4, Thermo Scientific, San Jose, CA, United States). The following constraints were used for the searches: tryptic cleavage after

3https://blast.ncbi.nlm.nih.gov/Blast.cgi

4http://www.uniprot.org

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. A false discovery rate (FDR) < 0.05 was considered as a condition for successful peptide assignments and at least two peptides per protein in at least one sample analyzed were the necessary condition for protein identification (**Supplementary Table S2**). Two biological replicates were used for each of uninfected and B. bigeminainfected ticks. For the quantitative analysis of tick proteins, after discarding Babesia and host proteins, the total number of peptide-spectrum matches (PSMs) for each tick protein was normalized against the total number of PSMs in ticks and compared between uninfected and infected ticks by χ 2 -test statistics with Bonferroni correction in the IDEG6 Software5 (p = 0.05). Only proteins that showed no significant differences between replicates and significant differences between the two conditions in at least three of the uninfected versus infected pair comparisons were considered as differentially represented (**Supplementary Table S3**). The mass spectrometry proteomics data have been deposited at the PeptideAtlas repository6 with the dataset identifier PASS013397 .

#### Gene and Protein Ontology Analysis

Gene ontology (GO) terms were assigned to differentially expressed genes and represented proteins with the Software Tool for Rapid Annotation of Proteins (STRAP) version 1.5 (Cardiovascular Proteomics Centre of Boston University School of Medicine at http://www.bumc.bu.edu/ cardiovascularproteomics/cpctools/strap/). For the analysis, two categories of GO terms were evaluated that included the Biological Process (BP) and Molecular Function (MF). Chord diagrams were generated to display GO and the expression/representation of apoptotic and stress responserelated UniProt IDs by using the GOplot R package in RStudio (Version 1.1.453) (Walter et al., 2015).

#### Proteomics Data Validation

Proteomics data validation was conducted by Western blot using two commercial mouse monoclonal antibodies targeting homologous amino acids sequences between tick and mouse and an in-house rabbit anti-subolesin polyclonal antibody (**Supplementary Table S4**). Ten engorged female ticks of each condition (uninfected and infected) were dissected, internal organs were rinsed in PBS and used to extract proteinaceous content using TriReagent (Sigma–Aldrich, St. Louis, MO, United States), according to the manufacturer's instructions. Elution was performed in 1× PBS, supplemented with 1% sodium dodecyl sulfate (SDS) and 1% protease inhibitor cocktail (Amresco, Solon, OH, United States). Protein concentration was assessed by spectrophotometry using a NanoDrop ND-1000 (Thermo Fisher Scientific). Fifty micrograms of protein extracts from each condition were re-suspended in Laemmli buffer (Bio-Rad, Hercules, CA, United States) containing 5% (v/v) of 2-β-mercaptoethanol, separated on a 12.5% discontinuous SDS–PAGE gel and transferred overnight at a constant 25 V to a nitrocellulose membrane, with a pore size of 0.2 µm (Bio-Rad, Hercules, CA, United States), using the Mini Trans-Blot R Electrophoretic Transfer Cell (Bio-Rad, Hercules, CA, United States). Membranes were stained with Ponceau S and blocked with 5% (w/v) nonfat dry milk (Bio-Rad, Hercules, CA, United States) in PBS containing 0.05% (v/v) Tween 20 (PBS-T) at room temperature (RT), for 1 h. In parallel, polyacrylamide gels were stained with BlueSafe (NZYTech, Lisbon, Portugal) for 1 h at RT. Membranes were washed with Tris-buffered saline complemented with 0.05% (v/v) Tween 20 (TBS-T) and incubated with the primary antibody diluted according to the manufacturer or literature (**Supplementary Table S4**), for 90 min at RT. After washing, membranes were incubated with a goat anti-mouse Polyvalent Immunoglobulins (G, A, M) – Alkaline Phosphate secondary antibody (1:3000; Sigma–Aldrich, St. Louis, MO, United States) for 1 h at RT, and afterward washed with TBS-T. The antigen–antibody complexes were detected using the alkaline phosphatase (AP) conjugate substrate kit (Bio-Rad, Hercules, CA, United States) and the protein band intensities were estimated using ImageJ Software (version 1.51K) (Schneider et al., 2012).

#### RNA Interference Assay Ticks and dsRNA Synthesis

For the RNAi assay, R. annulatus pathogen-free adult ticks were obtained from the laboratory colony established at Knipling-Bushland U.S. Livestock Insects Research Laboratory and Veterinary Pest Genomics Center, Texas, United States, and maintained on cattle at the tick rearing facilities at the University of Tamaulipas, Mexico, as described in the section "Rhipicephalus annulatus and Babesia bigemina." For the RNAi-mediated gene-silencing assays, gene-specific double stranded (ds) RNA was synthesized using as template the sequences assembled in this study and sequences publicly available in GenBank. In particular, Rhipicephalus pulchellus TSA GACK01000273, GACK01006332, GACK01002259, and GACK01006752, Amblyomma maculatum transcript JO843858 and R. microplus transcript U49934. Fragments of interest were amplified using the iProofTM High-Fidelity PCR Kit (Bio-Rad, Hercules, CA, United States), with specific primers containing T7 promoter (**Supplementary Table S5**) and previously synthesized cDNA. Briefly, PCR reactions of 50 µL were prepared with 1× iProof HF Buffer, 10 mM of dNTP mix, 0.5 µM of each primer, 1 mM of MgCl2, 0.02 U/µl proof DNA polymerase, and 1 µl of template cDNA. The cycling conditions were as follows: initial denaturation at 98◦C for 3 min, 35 cycles of denaturation at 98◦C for 10 s, annealing at specific temperature (**Supplementary Table S5**) for 30 s and extension at 72◦C for 15 s; and a final extension at 72◦C for 10 min. To confirm correct fragment amplification, gel electrophoresis was performed on a 0.5× TBE 1.2% (w/v) agarose gel, purified using the ZigmocleanTM Gel DNA Recovery Kit (Zymo Research, Irvine, CA, United States), and the fragments sequenced at StabVida (Caparica, Portugal).

5http://compgen.bio.unipd.it/bioinfo/software/

6http://www.peptideatlas.org/

7http://www.peptideatlas.org/PASS/PASS01339

The resulting sequence was compared against the sequences available at the NCBI nucleotide database8 . After correct confirmation of the amplified sequences, dsRNA was synthesized using the MEGAscript RNAi Kit (Ambion, Austin, TX, United States) accordingly with instructions provided by the manufacturer, purified, and analyzed by spectrophotometry and agarose gel.

#### Tick Inoculation and Cattle Infestation

Thirty pathogen-free R. annulatus unfed female ticks per group were injected with 1 × 1011 to 1 × 1012 molecules of specific dsRNA in the right spiracular plate, using a Hamilton syringe (Sigma) with a 1 in. × 33-gauge needle. The control group was injected with the same volume of elution buffer, similar to what has been described in other studies (Hussein et al., 2015; Ferrolho et al., 2017). While the use of unrelated dsRNA is now a well-established control in species with full genome annotation, the R. annulatus genome is not annotated thus increasing the possibility of off targets in silencing assays. After injection, ticks were kept in an acclimatized room at 22–25◦C and 95% relative humidity, for 6 h before infestation. Two 6 month-old calves were inoculated intravenously with 2 ml of cryopreserved B. bigemina-infected blood (1 × 106 infected erythrocytes) (field strain from Chiapas, Mexico). Cattle infection was confirmed by visual examination of blood smears. The

8https://blast.ncbi.nlm.nih.gov/Blast.cgi

cattle infestation was conducted as described in **Figure 1**. Ticks were monitored daily, at the second day unattached ticks were removed and the remaining were allowed to feed on the infected calves for 7 days. After this period, attached ticks were manually removed and stored in a humidity chamber until further use. SGs from 10 detached ticks per silenced-group were dissected and stored in RNAlater (Ambion, Austin, TX, United States). The remaining 20 female ticks of each group were kept in individual pierced 1.5 ml microfuge tubes, under the laboratory conditions described above, for oviposition.

#### Efficiency of Gene Silencing, B. bigemina Infection Levels, and Tick Fitness

Gene silencing efficiency in SG was determined by comparing mRNA levels between specific dsRNA-injected and control groups by qPCR as described in the Section "Transcriptomics Data Normalization and Validation." Quantitative PCR was carried out in both SG and ovaries for the 18S rRNA of B. bigemina, following a previously described protocol (Antunes et al., 2012), and the effect of gene silencing in infection was evaluated by a Student's t-test (p = 0.05). Engorged female weight (EFW), egg mass weight (EMW), and egg production efficiency (EPE) were determined for each collected tick and analyzed between silenced and control groups by Student's t-test with unequal variance (p > 0.05). The EPE was calculated according to the formula EPE = (EMW/EFW) × 100. Tick ability to feed was determined by the ratio of engorged female ticks collected

cloth chambers were glued to the calf skin (one day before infestation). In day 7, thirty specific-dsRNA inoculated female ticks were added to its respective feeding chambers together with thirty males. Two days after infestation (day 9) unattached ticks were removed from the patches. Ticks were monitored daily and allowed to feed on the infected calves for 7 days. After this period, attached ticks were manually removed and stored in a humidity chamber until further use (analysis of biological parameters and tissue dissection).

from the feeding chambers to the total number of female ticks added initially, and statistically compared with the Chi-square test (p > 0.05) with the null hypothesis that tick mortality was independent of gene silencing.

#### RESULTS

## R. annulatus Female Sialotranscriptome: Assembly and Annotation

All the SG tested from R. annulatus female ticks fed on a B. bigemina infected calve were found to be positive for infection, while no infection was found on the SG of female ticks fed on the naïve calve. This enabled the production of a pool of infected SG and an uninfected SG pool to proceed with RNA-seq. With 40,573 988 high-quality paired-end reads, two transcriptomes were assembled, from the control and the infected samples. The assembly resulted in 33,379 putative transcripts in the control sample (33,118 with length < 200 bp) and 30,435 for the infected one (30,301 with length < 200 bp). After the first step of annotation, only the alignments (BLAST) with e-value minor than 10−10 were selected to infer functional annotation, resulting in 16,564 transcripts for the control and 15,037 for the infected sample. From these, 6823 unigenes with 1,374,576 ± 1,787,133 (Ave ± SE) estimated counts per protein were identified in the control population and 6475 unigenes with 1,182,677 ± 1,494,052 (Ave ± SE) in the infected group. Moreover, 884 unigenes were found to be exclusive to the infected samples while 1232 to the control samples (**Supplementary Table S2**). GO analysis of the obtained sialotranscriptome resulted in 69 sequences without BLAST hits, 877 sequences with BLAST hits without annotation, 473 with mapping without annotation, and 6286 sequences were functionally annotated. **Supplementary Figures S1, S2** summarizes the functional annotation of the sialotranscriptome of R. annulatus female ticks. These multi-level charts show relevant functions from the general levels to the specific ones. The GO annotation highlights BPs such as single-organism cellular process, regulation of cellular process, gene expression or transport, and MF related to protein binding, hydrolase activity, or transferase activity. Using a BCV of 0.6, 360 unigenes were found to be differentially expressed (p < 0.05), whereas 123 were up regulated and 237 down regulated (**Figure 2**). Automated functional annotation was performed using Blast2GO and unigenes that did not retrieve GO terms were manually annotated. GO level 2 charts are shown in **Figures 3**, **4**. The most represented BP were cellular and metabolic processes with 125 and 133 sequences, respectively. These categories include all processes carried out at the cellular level (e.g., apoptosis or homeostasis) and chemical reactions and pathways, including anabolism and catabolism, by which living organisms transform chemical substances. Biological regulation and regulation of BPs were represented by 48 and 44 sequences. Response to stimulus, including processes such as stress, immune, or redox state responses, was represented by 30 sequences. BP, such as growth or detoxification with a single representative, were included in the "other" BP category. Regarding the assignments of MF, the classes with higher representation, binding and catalytic activity,

with 148 and 129 sequences, respectively, followed by transporter activity, with 26 sequences were found.

## Validation of RNA-Seq Results

To confirm differential gene expression, qPCR was performed on individual SGs. The levels of mRNA genes related to stress response encoding for the following proteins were analyzed: (1) small heat shock protein (HSP) II (UniProt ID:E0YPC0), (2) putative 60 kDa HSP (UniProt ID:L7M6W4), (3) putative heat shock hsp20 protein (UniProt ID:L7M7M1), (4) HSP 90 (UniProt ID:F1CGQ9), (5) Hsp70, putative (UniProt ID:L7M6A7), (6) putative selenoprotein k (UniProt ID:L7M5G4), (7) glutathione peroxidase (UniProt ID:Q2XW15), (8) thioredoxin peroxidase (UniProt ID:L7M1W2), (9) adenylate kinase isoenzyme 6 (UniProt ID:L7M323), (10) tumor rejection antigen-Gp96 (UniProt ID:B7QC85), (11) putative HSP (UniProt ID:L7MFS4), (12) serum amyloid A protein-like (UniProt ID:A6N9S3), (13) dual oxidase 1 (UniProt ID:B7PVC0), (14) putative HSP (UniProt ID:L7M4B9), and (15) putative heat shock transcription factor (PHSTF) (UniProt ID:L7MFL0). In addition, genes encoding for proteins related to apoptosis metabolic pathway were selected for validation of RNAseq by qPCR: (16) cytochrome c oxidase Subunit 1 (UniProt ID:A0A059VIA9), (17) putative death-associated protein kinase dapk-1 (DAPK-1) (UniProt ID:L7MKM3), (18) putative apoptosis inhibitor 5 (UniProt ID:L7LVN8), (19) bax inhibitor-1 related (BI-1) (UniProt ID:G3MPQ7), (20) putative apoptosis antagonizing transcription factor (UniProt ID:L7MAJ6), (21) mitochondrial ribosome small subunit component mediator of apoptosis dap3 (UniProt ID:L7M383) and (22) queuine trna-ribosyltransferase (QtRibosyl) (UniProt ID:L7M340), and

FIGURE 3 | R. annulatus SGs transcriptional response to B. bigemina infection based on GO BP (level 2) of encoded proteins. Red bars represent down regulated genes, green bars represent up regulated genes with statistical significance (p < 0.05).

(23) GP80 (UniProt ID:Q17174) related to proteolysis and lipid transport, respectively (**Supplementary Figure S3**). A moderate positive correlation between the mRNA levels by RNA-Seq and qPCR was obtained (Pearson's correlation coefficient 0.6387, p = 0.001).

#### Proteomics of the R. annulatus Female Ticks in Response to B. bigemina Infection

Proteomics analysis identified a total of 4,594 proteins in R. annulatus female ticks (**Supplementary Table S2**), where

80.08% (n = 3,679) were peptide sequences related to ticks while the remaining were mainly associated with the bovine host. **Supplementary Figure S2** summarizes the functional annotation of the proteome of R. annulatus in multi-level charts. When comparing uninfected and B. bigemina-infected samples, 406 proteins were found differently represented (p < 0.05) and those were annotated using Blast2GO and UniProt-related databases (**Supplementary Table S3**). The proteomics data showed that Babesia infection leads to a higher number of under-represented proteins (n = 212), in comparison to over-represented (n = 194) (**Figure 2**). GO analysis contributed for the annotation of 89.41% differentially represented proteins (n = 363), while only 43 proteins were classified as "Unknown" due to the absence of annotation. The level 2 GO terms for BP and MF are shown in **Figures 5**, **6**. The most represented GO terms were cellular, metabolic, and single-organism processes (BP, **Figure 5**), binding, catalytic, and nucleic acid-binding transcription factor activity (MF, **Figure 6**). Focusing on BPs identified herein, only cellular component organization or biogenesis, and cellular and metabolic processes, showed a higher number of overrepresented proteins than under-represented. Regarding the MFs, with the exception of structural molecule activity, all were found to be under-represented.

## Validation of Proteomics by Western Blot

For the validation of proteomics results, three antibodies against apoptotic and stress-related proteins, Ran, clathrin, and subolesin (4D8), were used for western blot analysis using proteinaceous extracts of fed-uninfected and fed-Babesia-infected ticks (**Supplementary Table S4**). The results corroborated proteomics data, showing that in response to B. bigemina infection Ran protein (UniProt ID: B7QIB6) is over-represented (Proteomics fold change = +∞; band intensity: NI = non-determinated, I = 5,849,430) while Clathrin (UniProt ID: B7PUK8) was down-regulated (Proteomics fold change = −3.7; band intensity: NI = 11,493.492, I = 5,064.865). The 4D8 protein was found to be over-represented in Babesia-infected ticks in comparison to the uninfected ticks (Proteomics fold change = +∞; band intensity: NI = 5,8119.247, I = 11,2021.586) (**Supplementary Figure S4**).

## Stress Response and Apoptosis: A System Biology Approach

All apoptotic and stress response-related molecules from transcriptomics and proteomics results were selected to further generate a chord diagram that displays the modulation of gene expression or protein production and its function (**Figure 7**). Eighteen and seven UniProt ID's were found to be related to stress response and apoptosis, respectively. GO clusters such as regulation of BPs, response to stimulus, biological regulation, detoxification, and antioxidant activity were intertwined with these two pathways. Consistent and significant results between transcriptomics and qPCR approaches allowed the identification of three molecules, of which two were up-regulated HSP and putative heat shock-related protein (UniProt IDs: L7MFS4 and L7M4B9, respectively) and DAPK-1 was found down-regulated (UniProt ID: L7MKM3). From this dataset, only two targets were

over-represented proteins and red bars under-represented proteins with statistical significance (p < 0.05).

identified in both transcriptomics and proteomics databases, a putative heat shock-related protein (UniProt ID: L7M4B9) and serum amyloid A protein-like (UniProt ID: A6N9S3), with specific modulation of their expression and representation. Four UniProt IDs demonstrated an opposite correlation between cellular protein levels and mRNA abundance, glutathione peroxidase putative heat shock-related protein, PHSTF, and putative apoptosis inhibitor 5 (UniProt IDs: Q2XW15, L7M4B9, L7MFL0, and L7LVN8, respectively).

#### Selection of Genes for RNA Interference

To gain insight into the functional role of identified tick vector genes/proteins in response to Babesia infection, targets for RNAi were selected based on their role in highly represented BPs such as apoptosis and stress response, and differential mRNA/protein levels in response to infection (**Table 1**). Two genes encoding for stress response-related proteins were selected (UniProt ID: L7MFS4 and L7MFL0), identified as being more expressed/represented in infected ticks. The first is a member of the HSP90 family, which is described as being able to regulate a specific subset of cellular signaling proteins that have been implicated in disease processes, particularly in ticks where they have been shown to have a role in virus replication (Weisheit et al., 2015). Transcriptional activation of heat shock genes relies on specific regulators such as the HSTF. Subsequently, a reduction of expression of such regulators may have a higher impact on activation of cellular stress response. Evidence of apoptosis manipulation in tick host cells by infecting pathogens has been demonstrated previously (Ayllon et al., 2015; Alberdi et al., 2016), which lead to the parallel selection of DAPK-1 and BI-1 genes (UniProt ID: L7MKM3 and G3MPQ7). The first, a pro-apoptotic kinase initially identified in vertebrate models, is involved in autophagy and tumor suppression (Inbal et al., 2002) and correspondent mRNA levels were found to be down regulated. The second, with an anti-apoptotic role (Lisak et al., 2015), appears to be under-represented in the current study. In the invertebrates Drosophila and C. elegans the role of DAPK-1 and related kinases has been briefly addressed revealing additional functions (Chuang and Chisholm, 2014). In ticks no functional studies have been performed so far with these two targets. The protein GP80 (UniProt ID: Q17174) is acknowledged to have a role in lipid transport and vitellogenesis (Tellam et al., 2002) thus correspondent reduction in mRNA levels is expected to deeply affect tick physiological parameters. QtRibosyl (UniProt ID: L7M340) is implicated in transfer RNAs (tRNAs) hypermodification, which in its turn are central players in nuclei acid translation (Lorenz et al., 2017). Particularly, queuine or Q-tRNA participate in many cellular functions, such as cell proliferation inhibition and stress (Vinayak and Pathak, 2009). In ticks, information of the impact of QtRibosyl depletion in tick physiological parameters, as well as in pathogen invasion, is still missing.

## Functional Studies

#### Assessment of Gene Silencing and B. bigemina Infection

The results of gene knockdown efficiency in SG and its effect on B. bigemina infection in both SG and ovaries are shown in **Table 2**. Silencing was achieved for target genes gp80 (98.6%), bi-1 (92%), hsp (88.5%), and QtRibosyl (71.4%), as opposed to

expressed/represented; gray square: non-significant (p > 0.05) and white square: not applicable.

TABLE 1 | Targets selected for RNAi assays, BPs in which are involved, and their representation on transcriptomics and proteomics analysis.

Red color indicates a significant down regulation or representation of the genes/proteins, the green color indicates a significant up regulation or representation of the genes/proteins while the gray color corresponds to non-significant regulation of the genes/proteins. ∗p < 0.05.

phstf. Despite testing different conditions, levels of dapk-1 mRNA were not detected on the dsdapk-1 inoculated ticks, suggesting an efficient knockdown. Silencing of gp80, dapk-1, and bi-1 lead to a significant increase of infection levels in tick SG (p < 0.05) with a ratio of 2.29, 19.75, and 4.39, respectively. Knockdown of the remaining targets resulted in a non-significant decrease of infection levels of SG in comparison to the control group. Regarding the ovaries, a significant reduction of B. bigemina levels (p < 0.05) (ratio 0.16) was observed in response to gp80 knockdown and in the remaining groups a non-significant increase was detected.

#### Biological and Reproductive Parameters of Tick After RNAi

Biological and reproductive parameters of R. annulatus female ticks injected with specific dsRNA before feeding are shown in **Table 3**. In all cases, the results obtained demonstrate that dsRNA inoculation did not affect the tick feeding ability (Chi-square, p > 0.05). Knockdown of gp80, dapk-1, bi-1, and QtRibosyl resulted in significantly lower female weights after feeding, in comparison to their respective controls (p < 0.05), and also in lower weight of laid eggs (p < 0.05). Females belonging to the gp80, dapk-1, and QtRibosyl silenced groups showed a significant (p < 0.05) reduction in EPE in comparison with their respective controls, which was not observed in the bi-1 silenced group. Four fed females belonging to the QtRibosyl silenced group were not able to lay eggs. phstf dsRNA inoculation resulted in lower female weight, egg mass, and EPE, but no statistical analysis was performed since silencing was not demonstrated.

#### DISCUSSION

Phylogenetic studies on certain pathogens and their tick vectors revealed a deep-rooted co-evolutionary relationship (de la Fuente et al., 2010; Antunes et al., 2017). The origin of piroplasmids could be traced back over millions of years to the time when tick species emerged, suggesting that Babesia parasites co-evolved with early tick lineages and their vertebrate hosts (Silva et al., 2011). Some pathogens may have a modest impact

TABLE 2 | Efficiency of gene knockdown by RNAi and its influence on B. bigemina infection levels in R. annulatus female ticks.

Thirty female R. annulatus ticks per group were injected with specific dsRNA or elution buffer as control. Ticks were allowed to feed in separated patches on calves experimentally infected with B. bigemina. Gene knockdown was analyzed by qPCR in 10 individual SGs per group by comparing mRNA levels between specific dsRNAinjected and control ticks using the CFX ManagerTM Software by means of the Pfaff method, ∗p < 0.05. The B. bigemina infection levels were determined in tick SG and ovaries by qPCR of the pathogen 18S rRNA gene and normalized against tick 16S rRNA using the ddCq method (2−Cq target −Cqcontrol). Infection rate was calculated by the ratio of silenced per control groups. The mRNA levels and B. bigemina infection in ticks were compared between specific dsRNA injected and control ticks by a Student's t-test (∗p < 0.05). ds, double-stranded; ND, not demonstrated; a, mRNA levels not detected in the silenced group.

TABLE 3 | Evaluation of biological and reproductive parameters of R. annulatus female ticks injected with specific-dsRNA.

Thirty R. annulatus female ticks per group were injected with dsRNA or elution buffer as control. Ticks were allowed to feed in separated patches on calves experimentally infected with B. bigemina. All attached ticks were collected after 7 days of feeding, weighed, and held in an acclimatized room. Ten ticks per group were randomly selected and SG dissected. The remaining were allowed to lay eggs. Tick ability to feed was evaluated as the ratio of removed ticks to the total number of ticks placed on the lamb using Chi-square test (∗p < 0.05). Female tick weight after feeding, EMW, and EPE were compared between dsRNA and unrelated dsRNA ticks by Student's t-test. The EPE was calculated according to the formula EPE = (EMW/EFW) × 100. N, number of samples used in the statistical analysis of the subsequent parameters; EFW, engorged female weight; EMW, egg mass weight; EPE, egg production efficiency; ∗ , statistical difference when compared to the control group. aNo statistical analysis was performed due to no gene knockdown. bFour females didn't laid eggs.

upon their vector and/or host, while others greatly influence host fitness and may manipulate host gene expression to favor infection, dissemination, and transmission (Cen-Aguilar et al., 1998; Mercado-Curiel et al., 2011; Antunes et al., 2017, 2018; Cornejo et al., 2017). We have found in our previous studies that the presence of Babesia in ticks induces a transcriptional shift (Antunes et al., 2012), particularly in the SGs (Antunes et al., 2018). Herein, a systems biology approach was applied to gain a more refined understanding on the complex interplay between B. bigemina and its tick vector R. annulatus. Our results identified dapk-1, bi-1, gp80, and QtRibosyl genes, that have important roles in apoptosis and vitellogenesis, and demonstrated that Babesia infection fundamentally affects key processes in SG tick cells.

It has been well documented that apoptosis is a highly regulated form of cell death. This process ensures a cell or organ to respond adequately to stress, thus limiting the risk of a disease or an infection that may compromise the ability of an organism to survive. Since its discovery, apoptosis has been extensively studied and its molecular regulation is now well understood. One of the critical components of this cell death pathway is dapk-1, a pro-apoptotic kinase, which can also regulate autophagy. Typically, dapk-1 has been studied in the context of its suppression in tumor growth and metastasis but additional roles emerged in last years (Farag and Roh, 2019). In this study, we have shown that exposure of ticks to B. bigemina significantly down regulated the transcription of dapk-1 in tick SG cells. This might have a significant impact upon the survival and propagation of Babesia parasites in its vector host. Inhibition of apoptosis in ticks and tick cells toward infection and survival is a strategy previously described in Anaplasma infection (Ayllon et al., 2015; Alberdi et al., 2016). Furthermore, silencing dapk-1 led to a substantial increase in B. bigemina infection in the SG cells of ticks that fed on an infected calf. After a tick takes a blood meal, the SGs enlarge, then degenerate, and become atrophied. Recent evidence has demonstrated that SG transformation into a vestigial state is due to caspase-1-mediated apoptosis, limiting the potential for pathogen transmission (Yu et al., 2017). It is therefore interesting to observe the ability of B. bigemina infection to selectively restrict a vital gene within the apoptotic pathway. By blocking the apoptosis process, the parasite might very likely survive and persist in the tissue and further facilitate dissemination of the infection throughout the tick, and subsequent transmission to another host. In the event of apoptosis occurring in SGs, the parasite might be eliminated, likely via its destruction following the engulfment of the apoptotic bodies it is occupying. The dapk-1 silenced ticks had significant lower weight than the control ticks, after a blood meal which might be attributed to the increased virulence observed, as this knockdown was not a spurious consequence of elevated biting rates. Such impact in female ticks has been previously described after silencing apoptosis-related genes in A. phagocytophilum infection in the tick vector I. scapularis (Ayllon et al., 2013, 2015). To the best of our knowledge, our study is the first ever reported to identify dapk-1 from a high-resolution multi-omic screen of any Rhipicephalus tick species infected with Babesia parasites. Further research will seek to elucidate the molecular mechanism underlying dapk-1 suppression in these cells and its functional significance.

The BI-1-related protein (UniProt ID: G3MPQ7) belongs to the Bax inhibitory protein-like family, exhibiting a highly conserved BI-1 domain throughout eukaryotes and prokaryotes (Lisak et al., 2015). Sequence analyses showed a high homology to an Amblyomma triste putative growth hormone-induced protein (UniProt ID: A0A023GIS5) and a growth hormone-inducible transmembrane protein-like from Rhipicephalus zambeziensis (UniProt ID: A0A224YYC2), also known as transmembrane BAX inhibitor motif (TMBIM) containing protein 5 supporting its anti-apoptotic role (Reimers et al., 2007; Rojas-Rivera and Hetz, 2015). TMBIM5 is recognized to be the only member from the TMBIM family of proteins localized mainly in the inner mitochondrial membrane and directly involved in the outer mitochondrial membrane permeability. Down regulation of such a protein leads to the release of pro-apoptotic proteins from the mitochondria, whereas overexpression results in the stabilization of cytochrome c at the inner membrane (Oka et al., 2008). Herein, the R. annulatus putative ortholog protein was found to be down regulated in response to infection, suggesting that under invasion there is an induction of cell apoptosis. However, reinforcement of gene down regulation resulted in the increase of B. bigemina levels in the SGs. Interestingly, the expression of this gene was also found to be down regulated in many types of cancer suggesting that the functional impact of its regulation is not fully understood yet (Rojas-Rivera and Hetz, 2015). Our results point to interactions between Babesia parasites and the tick apoptosis-related molecules, which requires further investigation with a special reference to this particular metabolic pathway.

In ticks, vitellogenesis is the process of yolk formation with the material accumulating in the developing oocyte and characterized by the synthesis of vitellogenins (Vg) during the reproductive period (Coons et al., 1989). GP80 is a processed product from vitellogenin, encoded by one Vg gene (Thompson et al., 2007; Taheri et al., 2014), and for simplicity we will maintain, herein, the GP80 nomenclature. Multiple Vgs have been described in ticks (Boldbaatar et al., 2010; Rodriguez et al., 2016; Xavier et al., 2018). Sequence analysis demonstrate a high similarity of gp80 to Vg-1 (UniProt ID: A0A034WTV5). In this study, B. bigemina infection led to a significant down regulation of gp80 transcription, which is in contrast to our previous study results (Antunes et al., 2012). A possible explanation for this discrepancy relies on the different tissues analyzed in the two studies. In the present one, only SGs were targeted, whereas Vgs are considered to be absent and are thought to be functionally replaced by the hemelipoglyco-carrier protein (CP) (Donohue et al., 2009). Encountering transcription of the gp80 in this organ suggests that Rhipicephalus SGs possess genes highly similar to Vgs, sharing many molecular features to since similar results were attained in SGs of Rhipicephalus bursa (Antunes et al., 2018). Our data revealed a significant decline in engorged female R. annulatus weight, EMW, and EPE, when gp80 was silenced, upon feeding on an infected calf in accordance to previous studies (Boldbaatar et al., 2010; Antunes et al., 2018). Keeping

in mind that the silencing methodology is not tissue specific, our results suggest that a decreased expression of Vgs leads to a reduction in lipid transport and normal production of energy (ATP) provided by lipid disruption, leading to an abnormal or delayed development of ovaries and eggs. A consequent increase in B. bigemina infectivity in SG was observed by silencing gp80 and a parallel decrease was detected in the ovaries. This parasite sustains transmission via both horizontal and vertical routes, through feeding and eggs production, respectively, which imply that a parasite-driven strategy to restrict vitellogenin synthesis is less likely. This data suggest that the host may redirect resources away from cell defense against pathogens toward processes with a more urgent demand such as egg production. Therefore, once facing lack of Vgs transcription, which will ultimately affect progeny, cells efforts will go to compensate this process, leaving SG cells more prone to infection and pathogen dissemination. However, the opposite scenario was observed in the ovaries, whereas a decrease of the levels of B. bigemina was detected suggesting that parasite transmission is hampered by the GP80 deficiency. A similar effect has been observed in Laodelphax striatellus whereas the transmission of a virus is facilitated by the expression of an insect-specific Vg (Huo et al., 2018). In fact, it seems that the process of vitellogenesis is critical for efficient pathogen transovarial transmission and is has been suggested that Babesia may directly interact with a Vg receptor (Boldbaatar et al., 2008; Hussein et al., 2019). Targeting vitellogenesis to control tick infestation is a very attractive approach, since the inhibition or disruption of such process may result in tick female mortality, decreased egg production, and viability and, ultimately, have an effect on pathogens with transovarial transmission such as Babesia spp.

Transfer RNAs are primarily involved in translation, functioning as adaptors from DNA to proteins (Pathak et al., 2007; Lorenz et al., 2017). One of the most important tRNA hypermodification is the replacement of guanosine by an analog, the nucleoside queuosine (Vinayak and Pathak, 2009), resulting in queuine-tRNA (Q-tRNA) and implicating a cascade of enzymatic reactions (Betat et al., 2010). This Q-tRNA have been shown to correlate with diverse cell processes including stress tolerance, cell proliferation, tumor growth, and protein translation (Fergus et al., 2015; Tuorto et al., 2018). QtRibosyl, herein designated as QtRibosyl, is an enzyme involved in the chemical reaction that results in queuosine base modification (Deshpande and Katze, 2001; Fergus et al., 2015). In the present study, the expression of QtRibosyl was found to be up regulated while protein levels were found to be down regulated in response to infection. This may be due to several biological factors such as the occurrence of post-transcriptional and post-translational modifications but also due to methodological constraints (Maier et al., 2009; Ayllon et al., 2015; Villar et al., 2015). Interestingly, the human ortholog enzyme gene (Queuine tRNAribosyltransferase catalytic subunit 1-QTRT1) has been found differentially regulated in some cancer tissues. It has been shown that, as a consequence, the Q-tRNA deficiency is associated with different neoplastic tissues and may even correlate with tumor grade (Vinayak and Pathak, 2009; Fergus et al., 2015). A high expression of such enzyme suggests that the cell needs to produce Q-tRNA either to promote the activity of antioxidant enzymes, secure protein folding, or to maintain high glycolytic rate (Pathak et al., 2008; Fergus et al., 2015; Tuorto et al., 2018), processes which have been demonstrated to occur in infected ticks (Villar et al., 2015; Cabezas-Cruz et al., 2017). Moreover, the present study demonstrated that efficient knockdown of QtRibosyl affected tick parameters and led to a marked reduction of female weight, egg mass, and EPE. These findings are clearly supporting the role of Q-tRNA in cell proliferation and glycolytic metabolism. A non-significant reduction of Babesia infection in tick SG was also observed, which can be related to a deficit on cellular ATP pools, which in its turn results in a shortage of nutrients essential for parasite multiplication (Trager, 1974).

Stimulation of a cell stress response by pathogens in ticks has been previously reported (Mulenga et al., 2003; Villar et al., 2010; Busby et al., 2012; Weisheit et al., 2015). Herein, two heat shock-related genes were targeted: the hsp which was found to be up regulated upon infection and the phstf that showed a similar expression in SG on both conditions. Knockdown of hsp did not significantly affected the studied parameters, although further studies centered in this biologically relevant pathway in ticks should be performed to deepen our understanding on tick–Babesia interface.

Our study has pinpointed interesting tick proteins interacting with Babesia. Still, due to its transovarial and transstadial transmission capacity, future studies focusing on early feeding time points, other tick stages, or other tissues, such as ovaries, are needed to fully understand the molecular dynamics at the tick–pathogen interface.

## CONCLUSION

For many intracellular parasites, the host cell is no longer an opponent but rather an assistant, evolving to efficiently and precisely manipulate their host. Accordingly, in ticks, Anaplasma bacteria are able to subordinate cell BPs such as apoptosis and glycogenesis. The present study shed some light regarding the Rhipicephalus–Babesia interface by making use of transcriptomics and proteomics allied to gene functional analysis. The results obtained point out a strong parasitic pressure over cellular functions, of which processes like apoptosis and stress response stand out. The genes dapk-1, QtRibosyl, and gp80 were characterized with regard to their role on tick reproductive parameters and interaction with Babesia infection, and showed to be highly influential functional molecules. The importance of understanding the tick biology and tick–pathogen molecular interactions is of the outmost importance in the discovery of suitable vaccine candidates to control tick populations and tickborne diseases.

## DATA AVAILABILITY

The datasets generated for this study can be found in DDBJ/EMBL/GenBank and PeptideAtlas repository, Accession numbers GBJT00000000, GBJS00000000, and identifier PASS01339, respectively.

#### AUTHOR CONTRIBUTIONS

fphys-10-00318 March 30, 2019 Time: 16:7 # 15

SA, JdlF, and AD designed the study. VS, MM, and JMC were responsible for tick rearing. SA, JC, JF, MV, and JdlF performed transcriptomics and proteomics analyses. SA and JC performed qPCR assays and western blots. JMC, ND, SA, JC, and JF performed the RNA interference studies. SA, JC, JdlF, GS, and AD performed data analysis and wrote the manuscript. All authors edited and approved the final manuscript.

#### FUNDING

RHIBAB - PTDC/CVT/112050/2009 "Differential expression and functional characterization of tick (Rhipicephalus annulatus) genes in response to pathogen infection (B. bigemina)." SA is the recipient of a post-doctoral grant supported by FCT

#### REFERENCES

Bock, R., Jackson, L., De Vos, A., and Jorgensen, W. (2004). Babesiosis of cattle. Parasitology 129, S247–S269. doi: 10.1017/s0031182004005190

(SFRH/BPD/108957/2015), JC and JF are the recipients of Ph.D. grants supported by the FCT (SFRH/BD/121946/2016 and SFRH/BD/122894/2016, respectively).

#### ACKNOWLEDGMENTS

The authors would like to acknowledge Fundação para a Ciência e Tecnologia (FCT) for funds to GHTM – UID/Multi/04413/2013.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2019.00318/full#supplementary-material

tick are crucial for ovarian development. J. Insect. Physiol. 56, 1587–1598. doi: 10.1016/j.jinsphys.2010.05.019

Cornejo, E., Schlaermann, P., and Mukherjee, S. (2017). How to rewire the host cell: a home improvement guide for intracellular bacteria. J. Cell Biol. 216, 3931–3948. doi: 10.1083/jcb.201701095

Corson, M. S., Teel, P. D., and Grant, W. E. (2001). Influence of acaricide resistance on cattle-fever tick (Boophilus spp.) infestations in semi-arid thornshrublands: a simulation approach. Exp. Appl. Acarol. 25, 171–184. doi: 10.1023/A:<PHONE_NUMBER>962

de la Fuente, J., Estrada-Peña, A., Cabezas-Cruz, A., and Kocan, K. M. (2016). Anaplasma phagocytophilum uses common strategies for infection of ticks and vertebrate hosts. Trends Microbiol. 24, 173–180. doi: 10.1016/j.tim.2015.12.001

responsible for heme sequestration. Insect. Biochem. Mol. Biol. 37, 363–374. doi: 10.1016/j.ibmb.2007.01.004

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Antunes, Couto, Ferrolho, Sanches, Merino Charrez, De la Cruz Hernández, Mazuz, Villar, Shkap, de la Fuente and Domingos. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Modeling Modulation of the Tick Regulome in Response to Anaplasma phagocytophilum for the Identification of New Control Targets

Sara Artigas-Jerónimo1† , Agustín Estrada-Peña2† , Alejandro Cabezas-Cruz3 , Pilar Alberdi1 , Margarita Villar1 and José de la Fuente1,4 \*

1 SaBio, Instituto de Investigación en Recursos Cinegéticos IREC-CSIC-UCLM-JCCM, Ciudad Real, Spain, 2 Facultad de Veterinaria, Universidad de Zaragoza, Zaragoza, Spain, 3 UMR BIPAR, INRA, ANSES, Ecole Nationale Vétérinaire d'Alfort, Université Paris-Est, Maisons-Alfort, France, 4 Department of Veterinary Pathobiology, Center for Veterinary Health Sciences, Oklahoma State University, Stillwater, OK, United States

#### Edited by:

Abid Ali, Abdul Wali Khan University Mardan, Pakistan

#### Reviewed by:

Shahid Karim, University of Southern Mississippi, United States Snorre Stuen, Norwegian School of Veterinary Science, Norway

> \*Correspondence: José de la Fuente<EMAIL_ADDRESS>
†These authors have contributed equally to this work

#### Specialty section:

This article was submitted to Invertebrate Physiology, a section of the journal Frontiers in Physiology

Received: 21 February 2019 Accepted: 04 April 2019 Published: 18 April 2019

#### Citation:

Artigas-Jerónimo S, Estrada-Peña A, Cabezas-Cruz A, Alberdi P, Villar M and de la Fuente J (2019) Modeling Modulation of the Tick Regulome in Response to Anaplasma phagocytophilum for the Identification of New Control Targets. Front. Physiol. 10:462. doi: 10.3389/fphys.2019.00462 Ticks act as vectors of pathogens affecting human and animal health worldwide, and recent research has focused on the characterization of tick-pathogen interactions using omics technologies to identify new targets for developing novel control interventions. The regulome (transcription factors-target genes interactions) plays a critical role in cell response to pathogen infection. Therefore, the application of regulomics to tickpathogen interactions would advance our understanding of these molecular interactions and contribute to the identification of novel control targets for the prevention and control of tick infestations and tick-borne diseases. However, limited information is available on the role of tick regulome in response to pathogen infection. In this study, we applied complementary in silico approaches to modeling how Anaplasma phagocytophilum infection modulates tick vector regulome. This proof-of-concept research provided support for the use of network analysis in the study of regulome response to infection, resulting in new information on tick-pathogen interactions and potential targets for developing interventions for the control of tick infestations and pathogen transmission. Deciphering the precise nature of circuits that shape the tick regulome in response to pathogen infection is an area of research that in the future will advance our knowledge of tick-pathogen interactions, and the identification of new antigens for the control of tick infestations and pathogen infection/transmission.

Keywords: regulome, transcription, tick, Ixodes scapularis, Anaplasma phagocytophilum, ISE6 cells, vaccine

## INTRODUCTION

Ticks (Acari: Ixodida) are major vectors of pathogens affecting human and animal health worldwide, and consequently the focus of research for developing novel control interventions (de la Fuente, 2018). Among tick-transmitted pathogens, Anaplasma phagocytophilum (Alphaproteobacteria: Rickettsiales) is mainly transmitted by Ixodes spp. and the causative agent of human and animal anaplasmosis and tick-borne fever in small ruminants (Severo et al., 2015).

**69**

Recent developments in tick genomics have advanced research using latest omics technologies for the characterization of tickhost-pathogen interactions and the identification of candidate protective antigens (de la Fuente et al., 2016c,a, 2017; Gulia-Nuss et al., 2016; Shaw et al., 2017; de la Fuente, 2018). Vaccinomics, a holistic perspective based on the use of omics technologies and bioinformatics in a systems biology approach for the characterization of tick-host-pathogen molecular interactions is our platform for the identification of candidate vaccine antigens (de la Fuente and Merino, 2013; de la Fuente et al., 2016a, 2018; Contreras et al., 2017). In this context, tick cell lines constitute a valuable resource because it is a proven model for the study of tick-pathogen and particularly tick-A. phagocytophilum interactions, easy manipulation without animal experimentation, and the fact that A. phagocytophilum infects mainly one cell type in vertebrates (neutrophils) but multiple cell types in ticks better resembled by these cell lines (Munderloh et al., 1994; Severo et al., 2015; Villar et al., 2015; Bell-Sakyi et al., 2018).

The regulome (transcription factors-target genes interactions) and interactome (protein-protein physical and functional interactions) play a critical role in cell response to different stimuli including pathogen infection. Both regulome and interactome are implicated in transcriptional regulation, which is one of the most fundamental mechanisms for controlling the amount of protein produced by cells under different environmental and physiological conditions and developmental stages (Gronostajski et al., 2011; Vaquerizas et al., 2012; Shih et al., 2016; Rioualen et al., 2017). Therefore, the application of regulomics and interactomics to host/tick-pathogen interactions would advance our understanding of these molecular interactions and contribute to the identification of new control targets for the prevention and control of tick infestations and tick-borne diseases (de la Fuente et al., 2018; Artigas-Jerónimo et al., 2018a,b; Estrada-Peña et al., 2018).

Few studies have addressed the role of the regulome or regulon (part of the regulome including a set of genes that share a common regulatory element binding site) in the interaction between tick-borne pathogens and vertebrate hosts (i.e., Bugrysheva et al., 2015; Boyle et al., 2019). However, limited information is available on the role of tick regulome in response to pathogen infection (Artigas-Jerónimo et al., 2018b).

In this study, we applied complementary in silico approaches to modeling how A. phagoctophilum infection modulates tick vector regulome, and the possibilities for the identification of new control target antigens. This proof-of-concept research provided new information on tick-pathogen interactions and potential targets for developing interventions for the control of tick infestations and pathogen infection.

#### MATERIALS AND METHODS

#### Datasets

The RNA sequencing (RNAseq) datasets of differential expression of I. scapularis transcription factors (TF) and target genes (TG) in response to A. phagocytophilum infection was obtained from previously published transcriptomics analyses in ISE6 cells, and fed adult female midguts and salivary glands (Ayllón et al., 2015; Villar et al., 2015). Gene ontology (GO) level-3 annotations for biological processes (BP) were conducted using Blast2GO software (version 3.0)1 (Villar et al., 2014; **Supplementary Dataset 1**). The RNAseq data is available at http://dx.doi.org/10. 5061/dryad.50kt0 and http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE68881.

#### Network Analysis of the Tick Regulome in Response to Infection

A network of interactions followed by a co-correspondence analysis (CoCA) was used for the integration of TF and TG interactions (regulome) of I. scapularis tick response to A. phagocytophilum infection. The methodology to build the network of interactions between proteins and functional metabolic processes has been previously described and validated (Estrada-Peña et al., 2018). This network consists of a set of nodes that are connected by edges where nodes are the interacting items, and links between nodes represent the strength with which they interact. In this development, a TF or TG is the source node and the cell metabolic process(es) in which it is involved are the target(s). The edge linking both nodes has a weight, which is the expression of either TF or TG. Networks were built separately for infected and uninfected I. scapularis ISE6 cells, salivary glands and midguts. Only TF and TG with GO functional annotations were included in the networks (**Supplementary Dataset 1**). Centrality is a fundamental property of a network because it refers to nodes that connect high score nodes (Opsahl et al., 2010; Estrada-Peña et al., 2018). In this context, "high score" applies to other nodes with high importance in the network. We calculated the importance of a node in the "traffic" between different nodes of the network using Betweenness Centrality (BNC), giving a higher score to a node that sits on many shortest paths of other node pairs (Barthelemy, 2004; Estrada-Peña et al., 2018). In our context, it is an indicator of the relative importance of a TF/TG in the links between two or more processes. Other calculated indexes included the PageRank (PR), a measure of the importance of the nodes linking with a given node, and the Weighted Degree (WD), which was calculated from the expression profile of each TF/TG linking to a cell process (Estrada-Peña et al., 2018).

The interactions between TF and TG were demonstrated using CoCA. Only TF/TG with values of BNC or PR higher than zero were included. We did the CoCA using the indexes of centrality obtained from the network explained above. The function "coca" of the package "cocorresp" was used for the R programming environment (Simpson, 2016). Data on BNC, PR and WD of each infected and uninfected datasets from ISE6 cells, salivary glands and midguts were entered into separate CoCAs. The analysis aimed to relate two different datasets from uninfected and infected samples to find patterns that are common to both and associating the TF and TG that are close in the reduced multivariate space and establishing correspondences. The plotting of the scores in the two first axes of the reduced space gives the interaction between TF and TG, i.e., the closer

1www.blast2go.com

they are in the space, the higher is the expected interaction. The method produces a cloud of interacting TF and TG. To improve the resolution of the charts, we plotted only TF/TG that were at a maximum of two score units of distance. We assumed that other TG separated by more than 2 score units from the values of TF were not interacting with these TF.

#### In silico Prediction of TF-TG Interactions

Putative DNA binding sites in TG for TF present only in infected tick ISE6 cells were predicted based on published information for TF-interacting sequences in other species, and the I. scapularis genomic scaffold whole genome shotgun sequence using cisTargetX2 and direct search for TF binding sequences in the predicted 50 gene regulatory regions of the I. scapularis genome (**Supplementary Dataset 2**; Potier et al., 2012; Rougemont and Naef, 2012; Vaquerizas et al., 2012). The BP with higher representation in the upregulated than in downregulated regulome (peptidase inhibitor and stress response) in response to infection were selected for further characterization of TF-TG interactions.

#### RNA Interference (RNAi) for Gene Knockdown in Tick ISE6 Cells

The I. scapularis ISE6 cells (provided by U.G. Munderloh, University of Minnesota, United States) was maintained in L-15B300 medium as described previously (Munderloh et al., 1994). Four different TF (HSF, Ap-2, Arx and Hox; **Supplementary Dataset 2**) were silenced using two siRNAs for each TF (HSF: 5 0 GCA CUC AGG GCC AGG AUU A 30 and 50 CCU CGG AAG CAG ACA GGA A 30 ; Ap-2: 50 AGA AAG AGG ACA CGA AGA A 30 and 50 CCA AGA AAG AGG ACA CGA A 30 ; Arx: 50 CCA AGA AAG AGG ACA CGA A 30 and 50 GAC CGA AGC CAG AGU GCA A 30 ; Hox: 50 CCU CCA GCU UCA ACA CAU A 30 and 50 ACG CCA CGG CCG AGC UUA A 30 ) provided by Dharmacon (GE Healthcare Dharmacon Inc., Lafayette, CO, United States). As control, two Rs86 siRNAs (50 CGG UAA AUG UCG AAG CAA A 30 and 50 GCG AAU AUG AAG UCG GUA A 30 ) were used. The siRNA experiments were conducted by incubating ISE6 tick cells with 100 nM of each siRNA diluted in 100 µl of serum-free medium in 24-well plates using four wells per treatment. To facilitate siRNA transfection, DharmaFECT (GE Healthcare Dharmacon Inc.) was used following manufacturer's recommendations. After 24 h, 0.5 ml/well of fresh medium was added. After 48 h of siRNA exposure, medium containing siRNA was removed and replaced with 1 ml fresh medium alone or containing cell free A. phagocytophilum NY18 obtained as previously reported (de la Fuente et al., 2005). Cells were incubated for a total of 72 h, and then collected for DNA and RNA extraction.

### Determination of Gene Knockdown and TG mRNA Levels by RT-qPCR

Total RNA was extracted from ISE6 cells using All Prep DNA/RNA/PROTEIN Mini Kit (Qiagen, Hilden, Germany) following manufacturer's recommendations. Gene knockdown levels after TF RNAi were assessed for TF and TG by RTqPCR on RNA samples using gene-specific oligonucleotide primers (**Supplementary Table 1**), the Kapa SYBR Fast One-Step RT-qPCR Kit (Kapa Biosystems, Roche Holding AG, Basel, Switzerland), and the QIAGEN Rotor-Gene Real-Time PCR Detection System (Qiagen). 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. The mRNA levels were normalized against tick rps4 using the genNorm method [Delta-Delta-Ct (ddCt) method] as described previously (Ayllón et al., 2013). Normalized Ct values were compared between test siRNAs-treated tick cells and controls treated with Rs86 siRNA by Chi2 -test (p = 0.05; n = 4 biological replicates).

### Determination of A. phagocytophilum DNA Levels by qPCR

Total DNA was extracted from infected cells using an All Prep DNA/RNA/Protein Mini Kit (Qiagen, Hilden, Germany). DNA samples were analyzed by qPCR using gene-specific primers for A. phagocytophilum msp4 as previously described (Ayllón et al., 2013). Normalized against tick rps4 Ct values were compared between test siRNAs-treated tick cells and controls treated with Rs86 siRNA by Chi2 -test (p = 0.001; n = 2–4 biological replicates).

#### RESULTS AND DISCUSSION

#### Rationale and Experimental Design

The tick regulome in response to A. phagocytophilum infection was characterized in the I. scapularis tick vector to provide insights into tissue-specific regulome profiles, and the identification of potential targets for the control of tick infestations and pathogen infection/transmission. The experimental design included two independent methods for the in silico characterization of the tick regulome in response to A. phagocytophilum using transcriptomics data previously obtained from infected I. scapularis ISE6 cells, and fed female midguts and salivary glands (**Supplementary Figure 1A**). The first approach was based on a network analysis in which the nodes were either TF or TG together with their corresponding GO BP annotations, and the link between two nodes represented the expression of the gene (**Supplementary Figure 1B**). The indexes of centrality were calculated separately for each network of uninfected and infected samples, and only nodes of TF and TG with indexes of centrality higher than zero were used for co-correspondence CoCA analysis (**Supplementary Figure 1B**). The second approach was used in parallel with network analysis, and consisted in the in silico prediction of TF-TG interactions based on described TF recognition sequences by searching in the I. scapularis genomic scaffold whole genome shotgun sequence (**Supplementary Figure 1C**). This analysis was focused on TF present only in infected ISE6 cells, and TG in BP overrepresented in the upregulated than in the downregulated regulome in response to infection as a proof-of-concept to facilitate the

2https://omictools.com/cistargetx-tool

identification of candidate target antigens for development of vaccines and other control measures. The results of the network analysis were then plotted with TF and TG together in the reduced space to demonstrate that the position of the TF correlates with the TG that are near to these TF after the CoCA (**Supplementary Figure 1D**). Finally, the results of both approaches were compared, and the TF-TG interactions predicted by both methods were functionally characterized by RNAi in A. phagocytophilum-infected and uninfected tick ISE6 cells (**Supplementary Figure 1E**).

#### The I. scapularis Regulome Shows Tissue-Specific Signatures in Response to A. phagocytophilum Infection

For the construction of networks, a total of 144, 86, and 93 TF (**Supplementary Dataset 3**), and 5225, 3919, and 4341 TG (**Supplementary Dataset 4**) were used derived from tick ISE6 cells, salivary glands and midguts, respectively. Tick midgut did not show detectable differences in the network indexes of TF BP between uninfected and infected samples (close to 100% BNC; **Figure 1A**). However, the TF multicell development and anatomical structure development processes increased to near 200% in infected versus uninfected ISE6 cells (**Figure 1A**). The network centrality of all the TF processes showed a clear increase in infected salivary glands when compared to unifected controls (**Figure 1A**). The multicell development process was represented in TF from ISE6 cells only (**Figure 1A**).

Other than minor variations in the network indexes of the TF, each test showed different TF that were present or absent in either uninfected or infected samples. Four TF were detected only in uninfected ISE6 cells, and other 4 were recorded only in infected ISE6 cells (**Figure 1B**). The most prominent TF in ISE6 cells (ISCW01819) was completely inhibited in infected ISE6 cells (**Figure 1B**). The pattern was more complex in the salivary glands showing up to 17 TF recorded only in infected, and 4 in uninfected samples (**Figure 1C**). Nine TF were recorded only in uninfected and 8 only in infected midgut (**Figure 1D**). As in ISE6 cells, the three most highly represented TF in uninfected midgut were not recorded in infected samples (**Figure 1D**). The number of TG detected only in uninfected or infected samples varied from 197 (uninfected) to 206 (infected) ISE6 cells, 159 (uninfected) to 585 (infected) salivary glands, and 360 (uninfected) to 129 (infected) midgut.

Every detected TF that was unique for uninfected or infected samples was included in CoCA. Based on the values of network indexes of TG, 70 TG in uninfected and 88 TG in infected ISE6 cells, 5 TG in uninfected and 58 TG in infected salivary glands, and 43 TG in uninfected and 25 TG in infected midgut were included in the analysis. The results from multivariate analyses showed a clear correspondence between TF and TG recorded only in uninfected or infected samples while the origin of the first correspondence axis (value = 0) separated completely TF and TG occurring only in either uninfected or infected samples (**Figures 2A–C** and **Supplementary Figures 2A–C**). These analyses suggested that the TF closer to individual TG are likely to regulate the expression of these genes (**Figures 2A–C** and **Supplementary Figures 2A–C**).

These results evidenced tissue-specific differences between infected and uninfected cells, thus supporting previous findings at the mRNA, protein and metabolic levels in I. scapularis ISE6 cells, a model for hemocytes, midgut and salivary glands, which are involved in A. phagocytophilum life cycle in the tick vector (Ayllón et al., 2015; Villar et al., 2015; reviewed by de la Fuente et al., 2017). These results evidenced that the regulome regulates various BP involved in tick-A. phagocytophilum interactions (**Figures 1A, 3A** and **Supplementary Dataset 2**), a finding previously reported in other organisms (Shih et al., 2016; Casella et al., 2017).

#### A. phagocytophilum Modulates the Tick Regulome to Upregulate Biological Processes That Facilitate Pathogen Infection

To complement the network analysis approach to tick regulome study, the putative DNA binding sites were characterized in silico for TF and TG in the upregulated regulome in response to infection in tick vector ISE6 cells (**Figures 3A,B** and **Supplementary Dataset 2**). In particular, the peptidase inhibitor and stress response BP with higher representation in the upregulated than in downregulated regulome (**Figure 3A**) were selected for further characterization of TF-TG interactions (**Figure 3B**).

The results showed a correlation between complementary in silico approaches (**Figure 4A**), therefore providing support for the network analysis of the regulome to predict at the transcriptomics level the most significant TF-TG interactions in response to stimuli such as pathogen infection.

To characterize the functional implications of selected TF-TG interactions predicted by both methodological approaches (**Figures 3B**, **4A**), RNAi was used in ISE6 cells to knockdown the expression of TF and characterize the effect on TG mRNA and A. phagocytophilum DNA levels (**Figures 4B–D**). The results showed that after 77–87% (average ± S.D., 81 ± 4%) TF silencing (**Figure 4B**), the levels of all predicted TG except for ISCW011771 and ISCW012363 decreased when compared to Rs86 siRNA-treated controls (**Figure 4C**). The two TG that were not downregulated after TF knockdown had the lowest mRNA levels (**Figure 4C**), which could affect the possibility of detecting differences between test and control cells. Alternatively, other TF or interacting proteins could be involved in the regulation of these genes. Nevertheless, except for ISCW011771 and ISCW012363 the results supported the prediction that these TF are implicated in the regulation of TG. Nevertheless, these TF-TG interactions should be corroborated in future experiments using different in vitro protein-DNA binding assays (Yang, 1998; Forde and McCutchen-Maloney, 2002; Deplancke and Gheldof, 2012; Ogawa and Biggin, 2012).

Regarding A. phagocytophilum infection, the results showed a 96–97% decrease in pathogen DNA levels after TF knockdown (**Figure 4D**). These results suggested that the TF

FIGURE 1 | Changes in the expression of TF in tick ISE6 cells, salivary gland and midgut in response to A. phagocytophilum infection. (A) The percentage change of the index Betweenness Centrality (BNC) infection among uninfected and infected target organs in the 13 BP GO annotations. (B) Values of BNC of the 8 TF that showed the highest changes between uninfected BNC (u) and infected BNC (i) ISE6 cells. (C) Values of BNC of the 21 TF that showed the highest changes between uninfected BNC (u) and infected BNC (i) salivary glands. (D) Values of BNC of the 17 TF that showed the highest changes between uninfected BNC (u) and infected BNC (i) midgut. Abbreviations: reg., regulation; (u), uninfected; (i), infected.

and corresponding TG are upregulated by A. phagocytophilum to facilitate pathogen infection.

#### Characterization of TF and Upregulated TG in Response to Infection as Putative Control Targets

The TF implicated in the regulation of selected TG included heat shock transcription factor (HSF), Ap-2, Aristaless-related homeobox gene (Arx) and Hox (**Figure 3B**). These TF has been described before to function in different transcriptionally regulated processes in other species. The mammalian Ap-2 TF has been shown to be involved in transcriptional activation and DNA binding/dimerization (Williams and Tjian, 1991). The HSF family has been implicated in the regulation of different physiological processes including cell response to stress and infection (Gomez-Pastor et al., 2018). Hox and Arx are members of a family of essential developmental regulators that bind to homeodomains target DNA sequences to regulate embryogenesis and neuronal processes in different organisms (Pellerin et al., 1994; Cho et al., 2012).

Few studies in other host-pathogen models have shown that some of these TF facilitate pathogen infection and therefore has been proposed as potential targets for control interventions. In primary peripheral blood monocytes, HSF1 is upregulated by human cytomegalovirus (HCMV) for pathogen survival, and has been suggested as a potential control target (Peppenelli et al., 2018). The tick-borne pathogen, Ehrlichia chaffeensis, upregulates the expression of certain Hox genes to facilitate infection through epigenetic mechanisms in human monocytic leukemia cells (THP-1) (Mitra et al., 2018). However, the role of these TF during pathogen infection in ticks has not been investigated before.

The BP modulated by the regulome of selected TF-TG interactions with a putative role in facilitating A. phagocytophilum infection included peptidase inhibitor

and stress response (**Figure 3A** and **Supplementary Dataset 2**). The stress response upregulated TG included genes coding for peroxinectin and uncharacterized protein with heme binding and peroxidase activity (n = 2), and glutathione peroxidase with glutathione peroxidase activity. The peptidase inhibitor TG encoded for a carboxypeptidase inhibitor precursor with carboxipeptidase and metalloendopeptidase inhibitor activity, uncharacterized protein, Kunitz-type proteinase inhibitor 5 II, serpin-2 precursors and secreted salivary gland peptides with serine-type endopeptidase inhibitor activity (n = 6), and cystatin and salivary cystatin-L with cysteine-type endopeptidase inhibitor activity (n = 2).

These proteins are involved in key physiological processes during tick life cycle such as heme/iron metabolism and detoxification and tick-host interactions (Chmelar et al., 2016; de la Fuente et al., 2016a). Glutathione peroxidase belongs to the glutathione antioxidant defense system and protects eukaryotic cells from oxidative damage (Espinosa-Diez et al., 2015). Glutathione peroxidase levels and activity are affected in different ways by Anaplasma marginale infection in both vertebrate and tick cells (Reddy et al., 1988; More et al., 1989; Kalil et al., 2017; Esmaeilnejad et al., 2018). While the activity and expression of glutathione peroxidase and other components of the antioxidant system was lower in A. marginaleinfected cattle and water buffaloes (Reddy et al., 1988; More et al., 1989; Esmaeilnejad et al., 2018), glutathione peroxidase coding gene was upregulated in embryonic Rhipicephalus microplus BME26 cells in response to A. marginale infection (Kalil et al., 2017). These results showed that A. marginale infection induces a differential response of the glutathione antioxidant defense system in the vertebrate and tick hosts. RNAi-mediated gene silencing of glutathione peroxidase and other antioxidant defense system genes increased A. marginale infection in BME26 cells, suggesting that the antioxidant response mediated by this molecule might play a role in the control of infection in ticks (Kalil et al., 2017). Another of the identified tick TG, peroxinectin, is a cell adhesion protein involved in melanization of pathogens in invertebrates (Sritunyalucksana et al., 2001; Cerenius and Söderhäll, 2004), was upregulated in crayfish resistant to white spot syndrome virus, and susceptible crayfish failed to upregulate this gene in response to viral infection (Yi et al., 2017). Strong cellular adhesion in response to the invading agent during crustacean encapsulation defense reaction was proposed as a

FIGURE 4 | Functional characterization of selected TF-TG components of the tick ISE6 cells regulome. (A) The predictive results of network analysis and in silico prediction of TF-TG interactions were compared in infected ISE6 cells. The TF-TG interactions predicted by both methods (squared in black letter for TF and red letter for TG) were then functionally characterized by RNAi in tick ISE6 cells. (B) Percentage of TF gene knockdown with respect to Rs86 siRNA control in ISE6 cells. Normalized against tick rps4 Ct values were compared between test siRNAs-treated tick cells and controls treated with Rs86 siRNA by Chi2 -test (p < 0.05; n = 4 biological replicates). (C) The TG mRNA levels were determined by qRT-PCR in ISE6 cells after TF gene knockdown or treatment with control Rs86 siRNA. Normalized against tick rps4 Ct values (average + S.E.) were compared between test siRNAs-treated tick cells and controls treated with Rs86 siRNA by Chi2 -test ( ∗p < 0.001; n = 4 biological replicates). (D) The A. phagocytophilum DNA levels were determined by qPCR in ISE6 cells after TF gene knockdown or treatment with control Rs86 siRNA. Normalized against tick rps4 Ct values (average + S.D.) were compared between test siRNAs-treated tick cells and controls treated with Rs86 siRNA by Chi2 -test (∗p < 0.001; n = 2–4 biological replicates).

protective mechanism mediated by peroxinectin in infected crayfish. The heme-binding lipoprotein (HELP), a transporter of heme in ticks, was previously found to be upregulated and downregulated in midguts and salivary glands, respectively, of A. phagocytophilum-infected ticks (Villar et al., 2016). HELP, together with Vitellogenin 1 and 2, was proposed to transport heme to other tick tissues such as salivary glands (Hajdusek et al., 2009). The uncharacterized protein with heme binding activity identified in this study may have a function similar to HELP, suggesting that A. phagocytophilum affects hemoglobin primary cleavage and heme transport in tick midguts and salivary glands, possibly to regulate the levels of heme in a tissue-specific manner with potential effects for pathogen and vector survival. The expression of Kunitz-type proteinase inhibitors have been found to be modified in several tick species in response to infection by tick-borne pathogens such as Bartonella henselae (Liu et al., 2014), flavivirus (McNally et al., 2012), Babesia bigemina (Antunes et al., 2012), and A. marginale (Zivkovic et al., 2010). Kunitz peptides are moonlighting proteins that perform multiple functions within the feeding lesion (Schwarz et al., 2014). Upregulation of Kunitz proteins in salivary glands of ticks infected with B. henselae (Liu et al., 2014) and flavivirus (McNally et al., 2012) may be associated with host immunity modulation at the feeding site. In contrast to upregulation, it is less clear why Kunitz peptides, including a Kunitz-type proteinase inhibitor 5, would be dowregulated in salivary glands following A. marginale infection (Zivkovic et al., 2010). Additional studies show that the expression of Kunitz peptides is complex and may be related to the tick and pathogen species (Rachinsky et al., 2007; Antunes et al., 2012).