Varroa destructor parasitism has a greater effect on proteome changes than the deformed wing virus and activates TGF-β signaling pathways

Honeybee workers undergo metamorphosis in capped cells for approximately 13 days before adult emergence. During the same period, Varroa mites prick the defenseless host many times. We sought to identify proteome differences between emerging Varroa-parasitized and parasite-free honeybees showing the presence or absence of clinical signs of deformed wing virus (DWV) in the capped cells. A label-free proteomic analysis utilizing nanoLC coupled with an Orbitrap Fusion Tribrid mass spectrometer provided a quantitative comparison of 2316 protein hits. Redundancy analysis (RDA) showed that the combination of Varroa parasitism and DWV clinical signs caused proteome changes that occurred in the same direction as those of Varroa alone and were approximately two-fold higher. Furthermore, proteome changes associated with DWV signs alone were positioned above Varroa in the RDA. Multiple markers indicate that Varroa activates TGF-β-induced pathways to suppress wound healing and the immune response and that the collective action of stressors intensifies these effects. Furthermore, we indicate JAK/STAT hyperactivation, p53-BCL-6 feedback loop disruption, Wnt pathway activation, Wnt/Hippo crosstalk disruption, and NF-κB and JAK/STAT signaling conflict in the Varroa–honeybee–DWV interaction. These results illustrate the higher effect of Varroa than of DWV at the time of emergence. Markers for future research are provided.

. The samples consisted of emerging worker honeybees collected just as they started opening the capped brood cells. The bees were or were not parasitized by Varroa mites in the capped cells and did or did not possess visible clinical signs of DWV. Thus, the following four variants were collected from heavily Varroainfested colonies: CON -control, nonparasitized bees without DWV clinical signs; VAR -Varroa-parasitized bees without visible DWV clinical signs; DWV_VAR -Varroa-parasitized bees with visible DWV clinical signs; DWV -Varroa nonparasitized bees with DWV clinical signs. Important protein groups. It is apparent that most of the various changes in the proteomes were indicative of upregulation by DWV signs and Varroa parasitism, and some of the changes were even amplified or, in contrast, reduced by the interaction between these two stressors. In a limited number of cases, we observed suppression of the protein abundance of some markers. The markers that specifically changed in DWV and those that changed similarly in DWV and DWV_VAR but not in VAR were considered highly specific to the effect of DWV. Analogously, the changes observed to be associated with VAR or both VAR and DWV_VAR but not DWV were considered specific to Varroa.
The influence of pathogens and parasites is connected to changes in the host immune system; thus, examining immune-related mechanisms is key for research on this subject. Therefore, we inspected honeybee immune The positioning of CON, VAR and DWV_VAR in the same direction, with the DWV variant located apart from them, indicates that Varroa had a higher effect on proteome changes than that with signs of the virus. Furthermore, it appears that the combination Varroa parasitism with DWV causing clinical signs is additive, as the shift of DWV_VAR away from CON is approximately two-fold greater than the shift of VAR alone. Each individual dot represents an nLC-MS/MS analysis result. RDA1 explains 65.83% of the variability, while RDA2 explains 23.83% of the variability. The results of this analysis somewhat correspond to the heatmap in Fig. 3.
genes. For this purpose, we first examined proteins listed in a review related to honeybee viral defense 21 ; however, of the genes listed (see Table 1 in Brutscher et al. 21 ), which we converted to protein equivalents (Table S1f), we identified 9 proteins in NCBI after a detailed investigation. The identification of only a small number of proteins is in part a consequence of the number of records removed from NCBI. Of these 9 proteins found in our study, 5 showed substantial changes (Table S1g). Although these markers were important for our story, it was necessary to investigate further. Therefore, we selected other important immune-related proteins from the list of identified proteins and found changes in a number of proteins involved in various immune pathways, i.e., JAK/STAT, MAPK, JNK, Toll, LLRs, NF-κB, RIG-I, and autacoid regulation. A recent study showed the upregulation of genes in the Toll, Imd, JAK-STAT, JNK and RNAi pathways in response to a virus 22 , and activation of the Toll pathway by DWV is believed to impair NF-κB signaling 8,9 . Another study showed that JAK-STAT, mTOR, MAPK and endocytosis genes were upregulated due to Israeli acute paralysis virus (IAPV) 23 . Furthermore, a metanalysis of transcriptomes showed that Toll and Imd pathway genes were differentially expressed after Varroa/virus infection 13 . Although we can find similarities to the markers expressed with DWV and Varroa among these studies, it is necessary to consider the honeybee stage at which analysis is performed as well as the experimental design 13 . Here, we analyzed the emerging bee, and therefore, the results should be unique to this stage. Heatmap that visualizes the proteome differences of the four variants at the time of honeybee emergence. The presentation clearly demonstrates that the proteome changes in DWV_VAR were more different than in DWV and VAR compared to CON. Thus, it appears that combination of Varroa parasitism and DWV clinical signs impact the proteome more substantially than Varroa and DWV clinical signs alone. Note that the difference between VAR and DWV variants is not visible in the hierarchical heatmap clustering in the columns; however, the difference between VAR and DWV characterizes the shift along the RDA2 axis (see Fig. 2). This heatmap was created in the Perseus environment and we also show a heatmap created in R computed from the RDA analysis (Table S1zl). Importantly, both heatmaps showed same hierarchical clustering in columns.

Continued
When we investigated the changes in more detail, we were able to categorize the proteins into different functional groups (see Tables S1g-zh), which facilitated assessment. Because many proteins showed dramatic changes among variants, we worked with a threshold of at least a 0.8 log2 fold change when individually inspecting hundreds of key proteins. Final marker selection (markers with notes shown in Table S1zi) was used to generate (gene names provided in Table S1zj) a functional association network by STRING analysis, including some interconnected markers that changed to a lesser extent than a 0.8 log2 fold change. STRING analysis helped us generate conclusions about how key marker changes are connected by showing primary and linked markers. In brief, STRING analysis identified 284 of 336 selected proteins, for which the network (for stats, see Table S2) revealed key protein groups (Table S3) associated with the following KEGG pathways: metabolic pathways (Fig. S2); oxidative phosphorylation (Fig. S3), such as altered electron transfer chain (ETC) and ATP production; valine, leucine, and isoleucine degradation (Fig. S4), such as alternative ATP recruitment; fatty acid metabolism, degradation, and elongation (Fig. S5), such as the induction of fatty acids associated with virus replication; glutathione metabolism (Fig. S6), such as stress protection and virus replication; the Wingless (Wnt) and Hippo signaling pathways (Fig. S7), such as the regulation of organ growth; RNA transport ( Fig. S8), such as virus replication; and others. However, some important markers, such as phospholipase A-2 (Pla2a), which is a key component in the arachidonic acid metabolism pathway (Fig. S9), could not be incorporated to interconnect their proteins into the current network version. Finally, the STRING analysis indicated the importance of the thioredoxin-like fold InterPro domain (Fig. S10). The selected markers that resulted from our analyses are listed in Table 1 and are further discussed.  Table 1. Selected and discussed markers. The log2 fold change differences between CON, DWV, VAR and DWV_VAR and gene names, including the indication of connections between markers in STRING, are shown. No means that the gene name was not found and/or connection in STRING was not indicated.

VAR -CON
www.nature.com/scientificreports www.nature.com/scientificreports/ Excessive activation of the TGF-β pathway in the Varroa-honeybee-DWV interaction. In the latter investigation, we considered that Varroa saliva secretions impair host hemocyte-mediated wound-healing and plugging responses 24 . Among possible candidates, transforming growth factor beta 1 (TGF)-β1 or its analogs have been previously identified in tick saliva 25 . To ascertain whether there is a connection between the Varroa-induced protein changes observed in this study and TGF-β1, we compared our proteome-based changes with expression changes related to TGF-β1 treatment previously observed in a mouse mammary gland epithelial cell line 26 . Interestingly, we identified 51 unique homologous proteins in the VAR variant (compared via BLASTP) relative to those from a previous study 26 in which 50 markers were analogously changed (for details, see Table S4). In DWV_VAR, 49 proteins changed in a similar manner, and the one protein that differed was Pla2a. However, in the DWV variant, only 43 markers were analogously changed. This result indicated the presence of TGF-β1 or its influencer in Varroa saliva. However, notably, viruses have also been shown to influence TGF-β1 by either promoting or suppressing the TGF-β signaling pathway 27,28 . The upregulated DWV-specific marker, which likely promotes TGF-β, is chloride intracellular channel exc-4 (exc-4/clic4) 29 . Furthermore, increases in AP-1-β1 and AP-1-mu-1, which were highest in DWV and DWV_VAR, can represent a marker associated with TGF-β activation of vascular endothelial growth factor expression 30 . Therefore, we suggest that DWV positively regulates TGF-β signaling and that Varroa and DWV can have a synergic effect on TGF-β pathway activation.
The TGF-β pathway is impacted through modulation of eicosanoid metabolism 31 . Interestingly, the key player in the eicosanoid pathway, Pla2a 32 , was specifically downregulated in DWV_VAR, implying the suppression of eicosanoid generation by arachidonic acid associated with infection and injury 32 when DWV infection and Varroa parasitism of the host act together. Moreover, we suggest that the effect of Pla2 suppression can augment 15-hydroxyprostaglandin dehydrogenase [NAD + ] (Hpgd), which is known to inactivate prostaglandins and/or lipoxins 33 . Hpgd upregulation was not specific to Varroa or DWV, but importantly, its upregulation in DWV_VAR amplified the Pla2-downregulation-induced suppression of inflammation specific to DWV_VAR. Furthermore, the effect of TGF-β1 on proliferation is enhanced by eicosanoid inhibition 31 . Thus, we suggest that, in addition to the activated TGF-β pathway, the suppression of the eicosanoid proinflammatory component of the innate immune response is a previously unknown key factor in the Varroa and DWV interaction.
TGF-β induction has been reported in some cases to stimulate 34 or reduce 35 synthesis of the basement membrane-specific heparan sulfate proteoglycan core protein, commonly known as perlecan. Perlecan, as a multifunctional component of the extracellular matrix, shows different effects on distinct cell types 36 . Because we observed DWV-specific upregulation and opposing Varroa-specific downregulation (almost qualitative) of perlecan, including in DWV_VAR, specific, context-dependent cell migration and proliferation regulation is suggested; for example, perlecan blockade inhibits proliferation of vascular endothelial cells, and exogenous perlecan promotes the proliferation of keratinocytes (see Nakamura et al. 36 ). Perlecan is a known key factor in vascularization, sequestering growth factors such as fibroblast growth factor 2 37,38 , and its suppression by Varroa should thereby be connected to the wound-healing delay after injury 38 . In all treatments, the Varroa-induced perlecan effect can synergistically act with upregulated CD109 antigen, which inhibits the TGF-β1-induced-antiproliferative effect and wound closure 39 . Although CD109 attenuates TGF-β1 signaling in certain cells, it also promotes epidermal growth factor 40 . p38 modulation of p53 and excessive JAK/STAT signaling. At least one other important feature is prominently connected to the aberrant TGF-β signaling pathway. The aminoacyl tRNA synthase complex-interacting multifunctional protein 2 (Aimp2/p38) p53 proapoptotic factor 41 is a mediator of TGF-β signaling known to suppress the proto-oncogene c-myc 42 . Upregulation of Aimp2/p38 in DWV and DWV_VAR indicates the specificity of the effect of this marker due to DWV, but its effect was more significant in the stressor interaction. p53 is known as a general repressor of RNA polymerase III (Pol III) 43 , but in our study, general upregulation of the DNA-directed RNA polymerase III subunit RPC1 (Polr3a) indicated the opposite as well as induction of the RIG-I pathway connected with the induction of type I interferon (IFN-I) and NF-κB 44 . Polr3a activity can be linked to upregulated JAK/STAT signaling as an antiviral response, so the virus needs to modulate the JAK/STAT pathway 45 . Therefore, we suggest that DWV modulation of JAK/STAT occurs through p38/p53. Furthermore, upregulated casein kinase II subunit alpha (CkIIalpha) likely plays a role in p53-mediated apoptosis, and its activity should be activated by p38 MAPK 46 . However, it is important to note that CkII is multifunctional and plays a crucial role in PI3K/AKT/mTOR, NF-κB, JAK/STAT and Wnt 47,48 . Additionally, several other kinases, i.e., protein kinase C (Pkc) and stress-activated protein kinase JNK (Sapk/bsk), are able to phosphorylate p53 46 but are also of key importance in other pathways such as Wnt, which is discussed further. Notably, a number of other markers were indicative of p53 modulation, e.g., RANBP2-like and GRIP domain-containing protein 5/6 (RanBP2), SUMO-conjugating enzyme UBC9-B (Ubc9-b/lesswright), proliferating cell nuclear antigen (Pcna) and small ubiquitin-related modifier 3 (Smt3). While Aimp2/p38 upregulation was likely a response to the virus and seems to be disadvantageous for DWV, the upregulation of another marker, ubiquitin carboxyl-terminal hydrolase 5 (Usp5), may be favorable for the virus. We suggest that upregulated Usp5 is associated with p53 deactivation 49 .
p53 activity is connected to antiapoptotic Bcl-2 50 , which should mainly be downregulated in DWV and DWV_ VAR as a result of increased levels of PRKC apoptosis WT1 regulator (Pawr) 51 . Furthermore, p53 is known to induce the expression of the proto-oncogene BCL-6, which in turn suppresses p53 52 , and the Drosophila homolog of BCL-6, known as ken and barbie (ken), modulates the JAK/STAT pathway 53,54 . In all investigated treatments, we observed excessive JAK/STAT signaling through high upregulation of the signal transducer and activator of transcription 5B (Stat5B, also known as stat92e). JAK/STAT signaling upregulation is the key feature of the antiviral response in honeybees [21][22][23] ; however, the expression pattern of the Stat5B marker observed here suggests high JAK/STAT activation by Varroa as well and even attenuation of Stat5B in the stressor interaction. Furthermore, www.nature.com/scientificreports www.nature.com/scientificreports/ we mostly observed upregulated tetra-peptide repeat homeobox protein 1-like (tprx1) in DWV_VAR, which, according to a BLAST search, can also be annotated as BCL-6 corepressor-like protein 1-like (Bcl-6c), zinc finger protein 512B-like, or with a lower probability, cyclin-dependent kinase inhibitor 1 C (Cdkn1c/p57). We suggest the interesting idea that the marker tprx1/BCL-6c represses BCL-6 55 , thereby facilitating the activity of Stat5B, but the function of the putative BCL-6c must be investigated in future analyses. Importantly, BCL-6 exerts much stronger DNA binding than STAT and is able to repress transcription via STAT binding sites 56 . Opposing effects on target genes and inverse proportional expression by STAT5 and BLC-6 were observed 57 , further supporting our suggestion. Taken together, JAK/STAT was hyperactivated, and the p53-BCL-6 feedback loop was disrupted, so interestingly, p53 and JAK/STAT simultaneously had large impacts in DWV_VAR. Finally, we suggest that the p53 pathway impacts the JAK/STAT pathway via the STAT-masking mechanism 58 . Non-Smad pathways. Although TGF-β signaling induces the Smad2/3 pathway, it can also promote non-Smad pathways, such as MAPK/ERK (Ras-Raf-MEK-ERK), PI3K/AKT/mTOR and Rho-like GTPase 59,60 , and it is also linked to the markers described in the above section dedicated to p53. Importantly, non-Smad pathways modulate p53 activity; therefore, Smad and non-Smad pathways are integrated 61,62 . Furthermore, important MAPK/ERK markers, i.e., dual specificity mitogen-activated protein kinase kinase 3 and dual specificity mitogen-activated protein kinase kinase dSOR1 63 , seemed to be upregulated differently by DWV and Varroa. In this context, prosaposin, a relevant corresponding marker that was upregulated in DWV_VAR, activates MAPK pathways 64,65 through G-protein signaling 64 ; moreover, as part of the PI3K/AKT pathway, prosaposin induces Schwann cell survival 66 , which is an important feature of the response to viral and mite stress. The elevated G-protein-signaling in DWV_VAR nicely supports the upregulation of G-protein-signaling modulator 2. Furthermore, the activation of MAPK/ERK and PI3K/AKT/mTOR 67 corresponded to the mostly specific or DWV_VAR-specific upregulation of Ras. The Ras proteins include ras-like protein 2, ras GTPase-activating protein-binding protein 2, and ras-related proteins Rab-7a, Ral-a, Rab-1A, Rab-14; nonupregulation of ras suppressor protein 1 is also important in this context. The Rab proteins are essential in vesicle intracellular transport, and upregulation of these markers supports increased endo/exocytosis 68 . We stress that Rabankyrin-5, a FYVE-finger effector of Rab5 69 , and Rab-7a were both upregulated in all treatments and should mostly be associated with increased endocytosis, an important part of the virus life cycle 70 . The next set of important markers represents the ras-related p21 GTPases Rac1 and Cdc42, for which rho GTPase-activating protein 44 (Gap44) can act as a Gap. The Gaps switch off signal transduction 71 , and therefore the downregulation of Gap44 in DWV_VAR led to the promotion of the Rho-like GTPase signaling pathway. Moreover, serine/threonine-protein kinase Genghis Khan, which functions as a downstream effector, Cdc42, and regulator of actin polymerization 72 were upregulated in DWV_VAR. The upregulation of the regulator complex protein LAMTOR1 (p18/Lamtor1; late endosomal/ lysosomal adaptor, MAPK and mTOR activator 1) indicated mTORC1 activation, mainly in the stressor interaction 73 . Importantly, mTORC1 activation is amino acid dependent 74 , and leucine is crucial for mTOR activation 75 , which we found to increase metabolism; this event is discussed further in relation to the manipulation of energy requirements. Lamtor1 upregulation is likely a response to ROS accumulation, which is connected to p53-dependent apoptosis 73 .
Specific upregulation of RanBP3 and Src42A in the stressor interaction. Importantly, we observed DWV_VAR-specific upregulation of ran-binding protein 3 (RanBP3), which negatively regulates TGF-β signaling by interacting with Smad proteins 76 . However, RanBP3-mediated Smad2/3 nuclear export requires dephosphorylation by protein phosphatase 1 (Ppm1a) 76 , which should be obstructed by the upregulation of exc-4/clic4 in DWV_VAR 29 . Moreover, activation of the components of the non-Smad Ras and PI3K pathways signals RanBP3 phosphorylation and modulation of Ran-dependent nuclear transport 77 , and this event, including RanBP3 upregulation, can increase ribonucleoprotein nuclear transport, especially in the early and late phases of viral infection 78 . Next, we propose that a possible function of upregulated RanBP3 is to downregulate excessive JAK/STAT signaling in DWV_VAR through nuclear transport of the highly upregulated Stat5B 79 . RanBP3 is also known to function in the nuclear export of active β-catenin (armadillo in Drosophila), thereby leading to suppression of the Wnt pathway 80 . Finally, we indicate the disruption of the p53-BCL-6 feedback loop described above 52 , and we also suggest modulation of the p53-β-catenin feedback loop 81 , in which the RanBP3 marker may play an important role.
The next key marker that was specifically upregulated in DWV_VAR is tyrosine-protein kinase Src42A (Src42A), which is linked to many partners according to the STRING analysis (see Table S1zk); its activity can be associated with Ras, MAPK, JNK, Hippo including the PDZ domain, LLRs, cytokine receptors, serine/ threonine-protein phosphatase 2 A, DE-cadherin, acetyl-CoA carboxylase, EF-hand calcium binding protein, and exc-4/clic4. We identified Src42A as a specific theoretical partner for Stat5B (see analysis in Table S1zk), and this connection suggested the modulation of JAK/STAT signaling in DWV_VAR due to the consequent upregulation of Stat5B and Src42A as a link between Src kinase activation and the function of Stat92E, which has also been indicated in Drosophila 82 . Importantly, Src42A has been described to regulate the JNK signaling pathway and epidermal closure in Drosophila 83,84 , and upregulation of Src42A should be specifically associated with the activation of inflammatory cell signaling pathways in wounds 85 . It is also likely that Src42A upregulation correlates with upregulation of proteins connected to the elimination of reactive oxygen species (ROS; see Table S1h), particularly H 2 O 2 , which is detected by Src42A 85 . Finally, Src42A has been described to act as a negative regulator of RTK/Ras/Raf/MAPK signaling 86 . Thus, it is likely that upregulated Src42A functions to inhibit wound-induced transcription in epidermal cells 87 . www.nature.com/scientificreports www.nature.com/scientificreports/ Wnt and Hippo signaling. Upregulation of Src42A correlates with Sapk/bsk activation 84 , and the STRING analysis (Fig. S7) indicated that Sapk/bsk is a shared protein in Wnt and Hippo signaling. DWV_VAR-specific upregulation of MOB kinase activator-like 1 (mats) functions as an activator of large tumor suppressor (Lats)/ Warts (Wts) kinases in Hippo signaling, but it can be hypothesized that the additional effect of Src42A in DWV_ VAR is to abolish the effect of mats 88 . Hippo and JNK signaling regulate wound-induced polyploidization (WIP), and JNK through AP-1 limits rather than stimulating Yki activation and polyploidization in the Drosophila epidermis 89 . Considering that the activated AP-1/JNK pathway can promote expression of proapoptotic genes 90 , the upregulation of AP-1-β1 and AP-1-mu-1 is an important result. In the case of Wnt, it is important to consider that the canonical Wnt is β-catenin-dependent and that the noncanonical planar cell polarity Wnt (Wnt/PCP) and Wnt/Ca 2+ transduction cascades function independently of β-catenin 91 . Activation of Wnt/PCP is characterized by upregulation of Rac and JNK, Rho, and profilin, whereas the activation of Wnt/Ca 2+ is characterized by the upregulation of Pkc, calcium/calmodulin-dependent protein kinase II (Camkii) and calcineurin 91 . Moreover, extrinsic Wnt is connected to the polarity component of atypical protein kinase C (aPkc) 92 , which was mainly upregulated in DWV_VAR. A gain in aPkc function has been shown to transform growth by disrupting Hippo/ Yap (Yki) signaling 93 . This aberrant aPkc activity disrupts cell polarity, and its full activity requires partitioning of defective protein 6 (Par-6) containing the Cdc42-and Rac-interactive binding (CRIB)-PDZ domain 94 . Here, the PDZ-domain-containing uncharacterized protein LOC409904, which showed a similar expression pattern to that of Par-6, may also participate. Further, aPkc activity should be affected through its interaction with RanBP2/ Nup358 95 , which was mainly upregulated in DWV_VAR.
The DWV-upregulated Camkii is also a Wnt signaling member. This multifunctional kinase is best known to affect learning and memory but is also linked to Wnt activation 96 and can be connected to viral infection in the brain 97 , where DWV is localized 98 . Upregulation of calcyclin-binding protein in DWV_VAR enhanced the ubiquitin-mediated degradation of nonphosphorylated β-catenin 99 , and importantly, this event was connected to p53 activation observed here 100 . Rac1 is an important member of the Wnt signaling pathway that mediates β-catenin phosphorylation and stimulates the β-catenin-dependent transcription of Wnt target genes 101 ; activation of Rac signal transduction was described above. Furthermore, coronin-1C upregulation supports Rac1 activation 102 and is a marker of myosin II disassembly 103 . Additional markers include the upregulated Toll-like receptor (TLR) activator leucine-rich repeat flightless-interacting protein 2 (LRRfip2) 104 , which activates Wnt upstream of β-catenin 105 . Furthermore, slit homolog 1 protein-like has been described as a Wnt/β-catenin target 106 .
Considering the expression patterns of the markers, both canonical and noncanonical Wnt pathways were activated in the stressor interaction. Moreover, disruption of the crosstalk between Wnt and Hippo was indicated but may be cell specific. We note the recent finding that activated Hippo can paradoxically promote JNK-dependent cell migration 107 , but this event opposes macropinocytosis, which should result from the upregulation of aPkc and Src42A observed here 108 .
Lethal(2). The activation of Wnt/β-catenin signaling has been shown to be necessary for TGF-β-mediated fibrosis, but this activation is Smad independent and occurs through p38 MAPK 109 . Among the proteins showing the greatest changes in response to stressors was the upregulated protein lethal(2)essential for life (lethal(2)), which has an α-crystallin domain and is a member of the group of small heat shock proteins (SHSPs) that are induced by stress and are ubiquitous and analogous among vertebrates and invertebrates 110,111 . Recently, it has been suggested that lethal(2) has a function in antiviral defense 22 ; however, we suggest that the strong increase in DWV and even its amplified expression in DWV_VAR are consequences of the stress-responsive MAPK cascade because the protein is activated by the p38 MAPK cascade and MAPKAP kinases. Additionally, and importantly, the expression pattern of lethal(2) was similar to that of Aimp2/p38. Furthermore, many studies have shown a connection between α-crystallin/SHSPs and ischemic injury [112][113][114] . The overexpression of lethal(2) in honeybees is likely associated with DWV symptoms and is further amplified by Varroa, which sucks the host and inhibits wound healing, thereby restricting the blood supply to tissues. Overall, one alternative non-Smad pathway targeting p38 results in the excessive production of lethal(2) as a corresponding marker.
Varroa NF-κB downregulation overrides DWV activation. TGF-β1 has been reported to be an apoptosis activator that downregulates NF-κB 115 . However, NF-κB should be activated in response to infection by the virus, which hijacks NF-κB through modulation of various pathways 116 . Therefore, we inspected proteins connected to NF-κB regulation.
An important marker of the host inhibitory effect on DWV replication is upregulated hsp70-binding protein 1 (HspBP1), which inhibits viral gene expression mediated by NF-κB 117 . Thus, upregulation of HspBP1 in DWV and DWV_VAR may be considered a marker indicating virus suppression 117 and this effect was more significant in the stressor interaction. Furthermore, suppression of NF-κB illustrated significant upregulation of the key marker NF-kappa-B inhibitor cactus (Cact2) 118 in DWV_VAR. The association of Toll-mediated Hippo with downregulation of Cact2 119 indicates the possible conflict between the effects of Varroa and DWV in the stressor interaction. We hypothesize that the abovementioned Src42A-abolishing effect of upregulated mats participates here as a Hippo component 88 . Furthermore, mitogen-activated protein kinase-binding protein 1 ( = JNK-binding protein; Mapkbp1 = Jnkbp1) should downregulate NF-κB 120 , but it can also activate JNK as a result of the TGF-β activation 121 in VAR and DWV_VAR. In this connection, the inhibition of NF-κB in dendritic cells has been shown to induce strongly augmented JNK/AP-1 activity because of elevated levels of ROS 122 , illustrating the number of increased markers (see Table S1h). The elimination of dendritic cells is important in that it leads to immunologic indifference 122 . The expression patterns of AP-1s suggest that DWV drives immune indifference, but accounting for other relevant markers, the effect is greatest in the interaction with Varroa. (2019) 9:9400 | https://doi.org/10.1038/s41598-019-45764-1 www.nature.com/scientificreports www.nature.com/scientificreports/ Although we first mentioned strong indicators of NF-κB suppression, other markers indicate the opposite; the aforementioned Bcl-2 suppression by Pawr should be connected to NF-κB activation 123 . The TLR-related upregulated proteins connected to NF-κB activation are LRRfip2 104 , dorsal protein isoform (dorsal-1A) 124 , and an evolutionarily conserved signaling intermediate in the Toll pathway (Ecsit) 125 . Additionally, Ecsit is required for bone morphogenetic proteins (Bmps), which are members of the TGF-β family 126 . Another important marker change indicating activation of NF-κB is Polr3a, which has been shown to induce IFN-I 44 . Furthermore, the function of Pol III is promoted by activated Ras signaling 127 . Upregulation of these markers was observed in all treatments, but similar to the markers of NF-κB suppression, it was most significant in DWV_VAR.
Aberrant NF-κB also has connections to the proteins responsible for chromatin remodeling, mainly by acetylation/deacetylation 128 . Among these, we encountered dynamic gene regulation following DWV infection and Varroa parasitization. According to the three established differentiating enhancer-binding transcription factor models 128,129 , we recognized two principally affected regulatory classes: (i) upregulated JNK/AP-1 contributed to the above discussed histone acetylases/deacetylases in local chromatin modification and remodeling, and (ii) chief inducible factors suppressed NF-κB and upregulated Stat5B. At this point, we linked the above discussed BCL-6c and Stat5B to NF-κB such that BCL-6 and NF-κB cistromes mediate the opposing regulation of the innate immune response 130 . Note that if BCL-6c targets BLC-6, its overexpression can represent an attempt to boost Smad signaling 131 . Moreover, the similarity to Cdkn1c/p57 suggests a link to TGF-β, i.e., the cyclin-dependent kinase inhibitor specifically induced by TGF-β should be required for cell cycle arrest 132 . Furthermore, it is imperative to note that in the DWV and Varroa interaction we observed upregulation of RanBP3, which should be able to transport Stat5B and downregulate excessive JAK/STAT signaling 79 . Finally, NF-κB is the common TLR mediator 130 . Therefore, according to our results, the TLR-NF-κB subnetwork in DWV_VAR is likely specifically influenced by the impaired BCL-6 130 .
Hyperactivated STAT in DWV_VAR should induce supercompetitor cells that kill the loser cells that function independently of Wnt, Yorkie and Myc 133 . However, highly elevated JAK/STAT in DWV_VAR may be the result of Varroa attempting to downregulate NF-κB through TGF-β1 or an analog 115 . Thus, it appears that Varroa and DWV are conflicting in their regulation of the NF-κB pathway via suppression or promotion, respectively, during their interaction. Although a negative effect of DWV on NF-κB was recently described using a few markers 8,9 , our high-throughput experimental approach using emerging bees is different from the approach utilized in that study.
We conclude that Varroa-induced NF-κB inhibition prevails on DWV infection-induced NF-κB activation in the DWV_VAR variant. The competition between NF-κB activation or suppression is connected to excessive JAK/STAT signaling, which can be diminished by the RanBP3 marker.
Histone deacetylase Rpd3 as a conjunction marker. Based on the protein changes observed and the STRING analysis, we identified histone deacetylase Rpd3 (Rpd3) as a key marker. According to the STRING analysis, the function of HDAC1 ortholog Rpd3 is associated with upregulated Stat5B, Pcna, Ubc9-b/lesswright, Cact2, Polr3a, CkIIalpha, Rabankyrin-5, actin-related protein 1, enolase-phosphatase E1, downregulated Sumo3, and rather constant Arm/β-catenin. Considering the previous discussion, these members are associated with JAK/STAT, p53, NF-κB, RIG-I, PI3K/Akt/mTOR, Wnt, endocytosis, and amino acid metabolism. Similar to Rpd3/HDAC1, one of the above key markers, Src42A, is directly linked to Stat5B, Cact2, and Rabankyrin-5 according to the STRING network and thus has an additional impact on the corresponding pathways. Another important marker that is likely connected to Rpd3/HDAC1 is tprx1/Bcl-6c 55 . Although Rpd3 upregulation was not Varroa-specific, its effect was primarily highlighted by associated markers that were specifically or mostly changed as a result of the interaction between DWV and Varroa, in which the highest transcriptional activation by STAT was generated by chromatin remodeling via acetylation/deacetylation 134 . Similar quantitative changes to Rpd3 were observed for host cell factor 1, which has been described to activate viral gene transcription, scaffolding histone-modifying proteins and regulation of various cell cycle stages 135 . Thus, it is likely that hcf-1 influences transcription collectively with Rpd3. Overall, our results imply that histone posttranslational modifications drive the altered pathways. Relevant additional markers are, for instance, histone acetyltransferase type B catalytic subunit and protein arginine N-methyltransferase 1, which were mostly upregulated in DWV and DWV_VAR.

Surprising upregulation of hexamerin 70c and 110 due to Varroa and the stressor interaction.
Metamorphosis, which is characterized by total body reconstruction from a larva to an adult, is connected to extensive protein depletion until the first feeding after emergence 17 . The principal proteins used in metamorphosis are hexamerins 17 . Depletion of the storage proteins due to Varroa mite feeding on the host 17,136 or increased demands for protein synthesis due to virus infection has been suggested; however, recent transcriptome studies did not report important changes in hexamerin expression in relation to virus infection and/or Varroa parasitism 4,22 . Thus, the strong upregulation of hexamerin 70c and 110 associated with Varroa parasitism and the upregulation associated with the combined effects of DWV and Varroa observed in this study are very interesting. The opposite situation was found for hexamerin 70b, which was downregulated by DWV and Varroa and to a lesser extent by their interaction. Moreover, downregulation was observed for hexamerin 110 in the DWV variant. These results suggest contradictory mechanisms of action by the investigated stressors affecting the abundance of these proteins, which serve as crucial protein sources during metamorphosis. The increased abundance of hexamerin c and hexamerin 110 may be an attempt of the host to compensate for the loss of the key proteins needed for successful transition to the adult stage. The finding that hexamerin 70a levels were almost constant in all investigated variants supports our previous suggestion that this protein is not a storage protein used during metamorphosis 17 . Overall, the highly contradictory changes in the storage of hexamerins 70c, 110, and 70b due to DWV infection and Varroa parasitization resulted in amplification of the effects of Varroa parasitization, with this interaction likely leading to a fatal retarding effect on development. (2019) 9:9400 | https://doi.org/10.1038/s41598-019-45764-1 www.nature.com/scientificreports www.nature.com/scientificreports/ MRJPs, vitellogenin, and JH-related proteins. Additional nutrition-related proteins are the major royal jelly proteins (MRJPs) 137 and the reproductive protein vitellogenin 138,139 . However, and importantly, these storage proteins are also linked to immunity 137,138,140 . The results of this study showed upregulation of MRJPs in DWV, VAR, and DWV_VAR and downregulation of the vitellogenin precursor, mainly in DWV and DWV_VAR. The higher content of MRJPs, especially MRJP1 and MRJP3, is likely related to immune defense 140 . Interestingly, MRJP1 and MRPJ3 are also present in the honeybee brain 141 .
Importantly, low vitellogenin levels are connected to higher levels of juvenile hormones (JHs) 139 . The strong decrease in vitellogenin levels in DWV and DWV_VAR correlated to increased levels of farnesol dehydrogenase-like, which are connected to JH synthesis 142 , but upregulation of JH epoxide hydrolase 1, which inactivates JHs 143 , seems to be a counteracting mechanism and an attempt to decrease JH levels. Importantly, a decrease in vitellogenin has been found to be associated with apoptosis of hemocytes; this mechanism to decrease immunity has been described for worker bees, which dramatically downregulate their defense machinery when they switch from being hive workers to foragers 144 . Thus, the results indicate that vitellogenin suppression is driven by DWV to prevent viral destruction by hemolymph-based immune protection. This mechanism is even amplified by the interaction of DWV with Varroa, but it does not seem to be initiated by the parasite alone.

Manipulating energy requirements. The increase in energy requirements resulting from viral infection
is met by mitochondria. The upregulation of various proteins of the ETC (Table S1i), including the electron transfer flavoprotein (ETF), was consistent with the observation that viruses direct mitochondria to produce more ATP to supply energy 145 . Based on that connection, we further identified the strong upregulation of an isovaleryl-CoA dehydrogenase, mitochondrial (IVD), which is a mitochondrial matrix enzyme that catalyzes the third step in leucine catabolism, providing electrons to ETF 146 . Thus, high IVD upregulation is likely a mechanism connected to the increased need for ATP during nutrient starvation 147 . Moreover, IVD could hypothetically regulate mTORC1 activation by affecting leucine metabolism 75 . Finally, IVD upregulation may also be important for eliminating toxic isovaleryl-CoA 148 . Overall, we suggest that the elevated expression of IVD in response to DWV and Varroa plays a role in cell-mediated immunity against pathogens and consequently is an alternative generator of ATP in a starved, metamorphosing honeybee. One of the markers that should affect mitochondrial function is the abovementioned Ecsit 149 ; we suggest a possible link between this TLR/BMP pathway-related marker and upregulated oxidative phosphorylation, the inflammatory response, and developmental changes 149 .

ROS elimination and detoxification -distracted metabolism responses. Pathogen infection
is connected to increased production of ROS, which kills pathogens directly or through signaling pathways that activate the immune response, including the NF-κB pathway 150 . In our study, various proteins associated with ROS (Table S1h) reduction were upregulated, mainly thioredoxin-related proteins and peroxiredoxin-6. Therefore, ROS elimination due to DWV, Varroa, and their interaction seems to be primarily connected to these pathways.
Among the classes of proteins involved in detoxification (Table S1j), similar increases in the coexpressed cytochrome P450 4G11 (Cyp4g11) and NADPH-cytochrome P450 reductase in DWV and DWV_VAR are likely connected to the altered expression of cuticular hydrocarbons 151 . This observation is consistent with the effects on many proteins related to the cuticle (Table S1k). Next, among the glutathione S-transferases and UDP-glucuronosyltransferases that showed changes in expression, the greatest changes were observed for glutathione S-transferase 1 and UDP-glucuronosyltransferase 1-3-like, which were both increased similarly in DWV and VAR and were even amplified in DWV_VAR. This finding suggests that these proteins play important roles in the detoxification of endogenous compounds related to virus infection and Varroa parasitism.
An interesting finding related to detoxification was the upregulation of S-formylglutathione hydrolase (Fgh), which is an esterase D that is involved in the detoxification of formaldehyde 152 . Similar upregulation was observed for alcohol dehydrogenase [NADP( + )]-like (Adh). Since adenosylhomocysteinase controls methylation through regulation of intracellular adenosylhomocysteine (SAH) 153 , an increase in its levels should correlate with methanol and, consequently, formaldehyde production 154 . An increase in Adh may correlate with NADPH-dependent aldehyde elimination of alcohols 155 . According to the STRING analysis, Adhn activity is associated with other highly changed carbohydrate metabolizing enzymes, i.e., sorbitol dehydrogenase-like (Sdh) and trans-1 ,2-dihydrobenzene-1,2-diol dehydrogenase (Dhdh), and through sorbitol dehydrogenase, there is a further connection with L-lactate dehydrogenase (L-ldh) and Hpgd. The activated formaldehyde elimination pathway suggests an adaptation to formaldehyde resistance 156 ; this adaptation may be attributed to aldehyde accumulation due to injury 157 . Upregulation of Sdh indicates activation of the polyol pathway, which is connected to insufficient insulin, ROS increases, and activation of MAPK and NF-κB 158 . Additionally, Sdh is an important marker affecting injury 159 . We emphasize that the activated polyol pathway, which has a reduced ability to regenerate glutathione due to NADPH consumption, leads to a decrease in the defense against oxidative stress 160 . Dhdh is an important detoxification enzyme that can regulate ROS generation and the apoptotic pathway response 161 , which, based on its expression, was specific to DWV. The elevated qualitative/quantitative levels of L-ldh in DWV_VAR were consistent with the finding that this enzyme is a known marker of ischemic damage 162 . Thus, these enzyme markers merit further investigation because they illustrate the distracted metabolism connected to immune system disruption and injury. Apidermin 3. Important changes representing possible conflicts between the effects of DWV and Varroa are indicated by apidermin 3-like protein, which was downregulated in DWV but upregulated in VAR and even in DWV_VAR. Unfortunately, little is known about apidermin 3, a cuticular protein that is expressed in the (2019) 9:9400 | https://doi.org/10.1038/s41598-019-45764-1 www.nature.com/scientificreports www.nature.com/scientificreports/ exoskeletal epidermis and is uniquely associated with nonpigmented cuticles such as the eye cover and external cuticle of white pupae 163 . Due to the limited information available, apidermin 3 stands alone in our story.

DWV in the investigated samples -(un)certainty of detection.
On average, the results of the label-free nLC-MS/MS proteomics analysis indicated the highest virus load in the DWV variant (see protein hits in the Table S1zm), while in the remaining sample types (including control) it was similarly low. We stress that it is necessary to carefully consider the proteomic detection of pathogen proteins in the host to prevent the pitfalls in data evaluation 5,164,165 . Therefore, to verify the accuracy of DWV identification in the bee samples, we extracted the peptides related to the DWV polypeptide identification. Overall, we recognized 36 different peptides identifying the DWV-polypeptide (see Table S1zn). Analysis of these peptide sequences and alignment to the divided polypeptide 5,166 revealed that the peptides identified were both structural and nonstructural of the DWV polypeptide (see Table S1zn), indicating the DWV replication state 5 . The results show that the highest diversity of peptides identifying DWV was in the DWV variant; however, in the remaining sample types CON, VAR and DWV_VAR, quantitative identification of DWV was not reliable. Thus, in the future, it will be necessary to perform detailed quantitative analysis including qPCR on the presence of DWV in the emerging bees. Furthermore, it will be necessary to consider the influence of Varroa on DWV strain diversity 2,4,167 , which was not investigated in this study.
Summarization. In summary, using label-free proteomics, we demonstrated that the proteome changes in the interaction of Varroa and DWV were directed by Varroa (Fig. 2). It appears that DWV and Varroa influence biochemical pathways by analogous and opposing transitions, which is somewhat consistent with the concept ascertained in the tick-host-pathogen interaction 19 . We stress one important objection to the comparison to ticks in our study, which is that multiple Varroa parasitize a single infirm bee in a capped cell. Therefore, the bee is pricked many times by Varroa, and the changes caused by this stressor become more obvious. One can imagine that similar effects can be achieved in humans parasitized continuously by hundreds of ticks. Moreover, the specific impact of Varroa is that parasitism occurs during metamorphosis and the bee cannot eat, all underlined by that the mite is tightly associated with DWV transmission. The analogy to ticks is also supported by our finding that Varroa saliva contains tissue growth influencers 25 . We proved the impact of Varroa on TGF-β signaling, which impairs host immunity.
Based on proteomic changes, which were thoroughly inspected individually, we provide a simplified schematic illustration of the impacted biochemical pathways in the Varroa-honeybee-DWV interaction (Fig. 4). We indicate that Varroa and DWV compete to manipulate the immune system. Importantly, we note that Varroa and DWV have conflicting effects on activation and suppression of NF-κB, respectively. During their interaction, eicosanoid synthesis is suppressed, which augments the TGF-β pathway. We indicate that JAK/STAT is hyperactivated and that the p53-BCL-6 feedback loop is disrupted, so p53-induced apoptosis and JAK/STAT act simultaneously. Excessive JAK/STAT signaling likely occurs through the collective effect of stressors on JNK, i.e., TGF-β induces JNK and JNK activation due the host response to DWV. We hypothesize p53-dependent masking of STAT 58 and suggest that the prevailing effect of RanBP3 in the interaction is not to diminish the Smad-dependent TGF-β pathway but to participate in viral replication and diminish excessive JAK/STAT and/or Wnt signaling.
Although Varroa and DWV cooperate in some respects, they have opposing effects on NF-κB and JAK/STAT. In the stressor interaction, the markers indicate facilitated entry of viruses into cells, which, compared to simple transmission, is a case of cooperation.
This study provides a set of important markers that should be further investigated in particular situations. Because our study was conducted at the time of honeybee emergence, the time at which the mite effect wins over DWV remains undetermined.

Materials and Methods
Study design. For the study presented in this paper, it was pivotal to show whether mites and/or signs of DWV occurred in the capped cells. The sample design (see Fig. 1) was selected to determine the interaction between the Varroa mite and DWV. Thus, the collected samples comprised the following four bee variants collected at the time when they were just emerging, i.e., when they started opening the capped cells: (i) control bees that were not parasitized by Varroa in the capped cells and had no recognizable DWV impact (marked as "CON"); (ii) bees parasitized by the Varroa mite in the capped cells but with no visible DWV impact (marked as "VAR"); (iii) bees parasitized by the Varroa mite in the capped cells and showing visible signs of DWV impact (marked as "DWV_VAR"); and (iv) bees impacted by DWV but not parasitized by Varroa in the capped cells (marked as "DWV"). Note that in the study design, sample (iv) served as the second control, along with the healthy bees (i). Furthermore, note that the DWV variant was difficult to find relative to DWV_VAR, suggesting that the mite could parasitize bee larvae and that nurse bees could remove it before capping; therefore, the mite was not present in the capped cells.
Bee samples. The honeybees originated from the Bee Research Institute at Dol, Czechia. The samples, including controls, were collected from colonies infested by Varroa and with manifestations of DWV clinical signs. All samples were from naturally Varroa heavily infested colonies and had the same bee genetics, and no manipulative experiments were performed in this study. The samples were collected from 3 colonies and randomly selected for analysis. The honeybee workers were collected just as they started opening the capped cells. Four people attended the collection of samples from the comb: one person collected the bees with sterile tweezers and placed them into 0.5 mL Eppendorf ® LoBind tubes (Cat No. Z666491, Sigma-Aldrich, St. Louis, MO, USA), two others verified the number of mites present in the cells, and a fourth person recorded the data. The honeybee brood comb was taken from the Varroa-infested colony, and the collection time was limited to ca. 20 min to minimize the impact of the www.nature.com/scientificreports www.nature.com/scientificreports/ colony environment absence; then, the comb was placed back in the hive. The emerging bees were collected, and the mites, if present in the cells, were counted to obtain sample No./number of mites: 1/0, 2/0, 3/9, 4/7, 5/7, 6/12, 7/5, 8/4, 9/4, 10/7, 11/0, 12/0, 13/0, 14/0, 15/0, 16/0, 17/0, 18/0, and 19/5. Upon subsequent checking (Multiple Comparison p-value adjusted by Hochberg's method 20 , all p-values > 0.5), there were no significant differences among samples with different mite numbers 168 .The collected samples were transported to the laboratory in a rack on dry ice and further stored in a deep freezer at −80 °C until use. Proteomic analysis. The emerged individual bees were homogenized using a 2 mL Potter-Elvehjem tissue grinder (Kartell Labware division, Noviglio, Italy) and drilling machine (PSB 650RE, Bosch, Stuttgart, Germany). In brief, each bee was initially homogenized in 0.5 mL of 50 mM Tris-HCl, pH = 7.4 with 1% Triton X-100 and an EDTA-free protease inhibitor cocktail tablet (Cat No. 05 056 489 001, Roche, Indianapolis, IN, USA) dissolved in 30 mL of buffer. The sample was initially homogenized for 1 min, left on ice for 10 min, homogenized again for 1 min and then a final time after incubation for 15 min on ice in 2 mL of nonpure water (ddH 2 O; Thermo Fisher Scientific, Waltham, MA, USA). The samples were centrifuged at 20,000 × g for 20 min at 4 °C in an MR 23i centrifuge (Jouan Industries, France). The samples were then processed analogously to FASP, which is available online 169 .
All samples digested with trypsin were analyzed by nLC-MS/MS in two analytical replicates. All 19 biological samples were analyzed in one analytical series; thus, 38 nLC-MS/MS runs were performed consecutively. First, one test nLC-MS/MS analysis marked as 1_1 was performed; thus, the analyses of sample 1 included in the dataset are marked 1_2 and 1_3, and the rest of the analyses are marked X_1 and X_2 (where X is sample No. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Note that we originally selected 20 biological samples, with 5 biological replicates for each treatment, but the last two analyses in the series of 40 analyses lacked spectra and therefore were not evaluated. Nano reversed-phase columns (EASY-Spray column, 50 cm × 75 µm ID, PepMap C18, 2 µm particles, 100 Å pore size) were used. Mobile phase buffer A was composed of water, 2% acetonitrile and 0.1% formic acid. Mobile phase B was composed of 80% acetonitrile and 0.1% formic acid. The samples were loaded onto a trap column (Acclaim PepMap300, C18, 5 µm, 300 Å Wide Pore, 300 µm × 5 mm) for 4 min at 15 μL/min, and the loading buffer was composed of water, 2% acetonitrile and 0.1% trifluoroacetic acid, with analysis in a Thermo Orbitrap Fusion Tribrid mass spectrometer (q-OT-IT) instrument performed as previously described 170,171 . . Simplified schematic representation of the affected pathways in an emerging honeybee worker that impact Varroa and DWV symptoms via their interactions. Multiple markers indicated that Varroa mite salivary compounds promote TGF-β signaling, which strengthens Smad2/3 pathway signaling, but alternative signaling pathways are propagated. Connection of the activated Wnt and Hippo pathways, which undergo crosstalk disruption, leads to aberrant tissue growth, which is the likely cause of the increased probability of wing deformity. Considering that Varroa overrides DWV-induced proteome changes, as shown in Fig. 2, this suggests that the mite suppresses mechanisms that normally facilitate DWV infection. Thus, there are conflicts in JAK/ STAT and NF-κB signaling. We indicate the collective effect of Varroa and DWV on JNK, the domination of p53-induced apoptosis and that JAK/STAT is hyperactivated and the p53-BCL-6 feedback loop is disrupted. The cooperating effect is the facilitation of viral endocytosis by the mite. The activated Akt/mTOR (through upregulated p18/Lamtor1) pathway inhibits autophagy and recycling of cells damaged by ROS, which are activated by p53. The high number of damaged cells decreases the efficacy of ATP acquisition, resulting in loser cell status.
www.nature.com/scientificreports www.nature.com/scientificreports/ Proteomic data evaluation. The data were analyzed and quantified with label-free algorithms using MaxQuant software (version 1.5.3.8) 172 . The FDR was set to 1% for both proteins and peptides, and we specified a minimum length of seven amino acids. The Andromeda search engine 173 was used for the MS/MS spectra search against the NCBInr database downloaded on January 31, 2017. The actual database consisted of 150,263 records and contained not only the Apis mellifera genome but also pathogens (importantly viruses) and symbionts. The enzyme specificity was set as C-terminal to Arg and Lys, also allowing for cleavage at proline bonds and a maximum of two missed cleavages. MethylThio was selected as a fixed modification, and N-terminal protein acetylation and methionine oxidation were selected as variable modifications.
Statistical data evaluation. The proteomic data evaluated via MaxQuant 172 were further processed in Perseus (version 1.5.2.4) 174 . Contaminants, reverse, and only identified by site hits were removed, and the matrix was reduced to at least two positives in at least one group. The protein abundance data of intensities were logarithmically transformed as suggested by Anderson et al. 175 . We preferred log2 transformation, which permits better visualization of less-abundant species. The data were further analyzed using RDA and multiple comparison statistics using the "vegan" R package 168,176 . We evaluated whether it was beneficial to present our data including the technical replicates of the nLC-MS/MS analyses, because it allows us to show all variances influencing the results. Furthermore, this means of data presentation does not influence the statistical results from a mathematical perspective 168 . Distribution of the protein abundance to create the heatmap was computed from the abovementioned RDA analysis using the "predict" function in the "vegan" package 176 . The corresponding heatmap depicting the protein abundance and clustering in dendrograms was produced using the "gplots" R package 177 . In addition, we created a heatmap in the Perseus 174 environment. The missing values were replaced based on a normal distribution (width: 0.3; down shift: 1.7; total matrix) and the data were normalized through the Z-score. Hierarchical clustering was performed using the average Euclidian distance (Preprocess with k-means; 300 clusters) for both row and column trees. The resulting heatmap was adjusted to grayscale. While the first heatmap is provided in the supplementary material, the latter is shown in the main manuscript; however, both heatmaps showed the same hierarchical clustering in columns.

Individual data evaluation.
We individually inspected the proteins that showed the greatest changes. We worked with the threshold, which was characterized by at least a 0.8 log2 fold change, but we enriched the final list of markers with some interconnected markers that changed to a reduced extent or did not change. Note that some markers are important just because they do not change, while the other do. We further grouped the proteins based on their function. The selected proteins were searched using UniProt and NCBI, and the description, selected biological and/or molecular function or protein region or each protein was linked to the description to facilitate characterizations.
Data evaluation using STRING. To ascertain the functional association network, we identified genes for the selected proteins and entered them into STRING 10.5 178 . From the STRING network, we further selected the key KEGG pathways. It is important to note that not all proteins could be connected in STRING, but the analysis facilitated data evaluation.
Search for proteins that are influenced by TGF-β. Varroa parasitization is likely to affect host wound healing. However, whether this phenomenon occurs during the Varroa-honeybee interaction remains unknown. Therefore, we searched previously described tick salivary compounds, which have been reviewed 25 . Among the possible candidates, we selected tissue growth factor influencers 179 . This selection was facilitated by the observation that tissue growth factor influencers have been reported to interact with our key selected markers. Among the factors inhibiting healing, we identified TGF-β 25 . We further searched the literature related to protein changes caused by TGF-β and found a relevant publication by Xie et al. 26 , which we used to compare markers and changes. Each sequence of the protein markers presented in Tables 1 and 2 in Xie et al. 26 was processed via BLASTP 180 , and the list of proteins obtained using this process was searched for A. mellifera homologs. When the identified honeybee homologs displayed a homology higher than 50% and the same function, we determined whether the changes in the honeybee proteome corresponded to those reported by Xie et al. 26 .