We present an approach to the study of the biological process in the salivary glands of feeding ticks supported by the robust field of graph theory. The aim of this study is to present how the expression of proteins can be translated into biological processes, linking the partners into a construct known as a network or graph; a secondary focus is to outline the main differences of biological processes in feeding argasids and ixodids. Our purpose is to demonstrate the basis of a solid comparison of the changes of biological processes along the feeding period of a tick. The study is based on recent and high standard datasets, one of them referring to R. sanguineus s.l. (Tirloni et al., 2020) the other to O. rostratus (Araujo et al., 2019). After a preliminary filtering (see below) processes were obtained through queries to online repositories. To handle this data, we did use the graph theory framework, which establishes links between protein(s) and process(es). Therefore, significant indexes can be derived from the network. Since we obtained data for the two major families of ticks, we contrasted the differences in feeding types.

We conducted a literature search in public repositories, like Pubmed, Google scholar and Scopus following the PRISMA statements (Moher et al., 2009). In the compilation we selected those studies that were found under the search terms: (“tick(s)” OR Ixodidae OR Argasidae) AND (“salivary glands” OR “gut”) AND (“sialome” OR “sialotranscriptome” OR “proteome”) with special attention to the Supplementary Material of each study. The query produced a relatively low number of published papers; therefore, we also evaluated the bibliography of each paper to complete the dataset. Data was selected to contain some details, as follows: (I) the name of the protein had to be following standards in SWISS-Prot/Uniprot KB database for improving further queries (see below); (II) the protein reads had to be in TPM (transcripts per million); (III) the data should refer to “slices” of the feeding process (in ixodids); (IV) we aimed to include contemporaneous studies on argasid ticks. Aiming to comparison, it is important that the studies were carried out at similar periods of time, since fast pacing advances in these techniques may turn results matchless. Two studies fulfilled these standards, one of them dealing with Rhipicephalus sanguineus s.l., (Rs) and the other with Ornithodoros rostratus (Or). In the study dealing with Rs (Tirloni et al., 2020) data is partitioned according to the feeding moment using the weight of the tick, which is a better estimation than the elapsed time after the onset of feeding. The study determined seven feeding groups: unfed ticks; group 1 (G1) with an average weight of 1.8 mg, collected at day 2 of feeding; group 2 (G2), averaging 3.6 mg, collected at day 6 of feeding; group 3 (G3) averaging 7.0 mg, collected also at day 6 of feeding; group G4 (G4) averaging 10.9 mg, collected at day 8 of feeding; group G5 (G5) averaging 24 mg, collected at days 8 and 11 of feeding and group 6 (G6) averaging 36 mg, collected at days 6, 10 and 13 of feeding. The Or study (Araujo et al., 2019) is aimed to compare the proteins in both salivary glands and gut, considering the short feeding time of soft ticks. This selection did not want to exclude other datasets because its unreliability, but because the selected ones correctly matched the aims for developing our framework.

We performed a data cleaning of both datasets, including only proteins in the original spreadsheet about the complete transcriptome of the Supplementary Material by Tirloni et al., 2020, with a coverage equal or higher than 90% (column BM), an E value equal or less than 10^-5 (column BI) and a similarity equal or higher than 50% (column BL). A similar filtering was followed for Or, with the same parameters, using columns FT, FZ, and GA, respectively, of the original datasheet by Araujo et al. (2019). To note that a selection of proteins restricts their number, and therefore affects conclusions. Even more, both datasheets represent different experiments. The proteins fulfilling these criteria were queried to UniprotKB (UniProt Consortium, 2021). Uniprot is a collection of annotated proteins that is manually curated and constantly updated. We search the actual name of protein given by the original datasheets for Rs and Or, without performing a BLAST query, since the original datasheets contain these names. In short, the list of proteins to be queried was uploaded to Uniprot, using the “Retrieve/ID mapping” option, and marking the option “biological processes” in the selection of downloadable columns. This returned the set of biological process(es) in which each protein is involved, including the Gene Ontology value(s) associated to each protein (GO; Ashburner et al., 2000; Gene Ontology Consortium, 2021). GO is a database that assumes that if a gene has the same behavior as other in a specific category of gene ontology, the gene belongs to the same category. Not every protein had an annotation; a few proteins were therefore impossible to ascribe to any biological process; these proteins (less than 1%) were dropped from the study. All the biological processes have its own GO term associated number, from the most specific biological process to the most general one, which is also returned as a “cascade” of processes by Uniprot. This allows to correlate “child” and “parent” processes in the chain of events by GO (i.e., “DNA repair” is a part of “Metabolism”, but other child processes are also part of “Metabolism”; therefore a network approach makes sense, see point 2.3 below). We did a distinction between low-level (LL) processes (the “children”) and top-level (TL) processes (the “parents”). To obtain the complete chain of processes involved in the expression of every single protein we used the package “GO.db” (Carlson, 2019) for the R programming environment (R Core Team, 2020) that reads and handles the GO terms returned by Uniprot for each protein.

A network is a set of nodes that interact through links. Networks can be used to represent i.e., food webs (Ings et al., 2009) associations of parasites with hosts (Krasnov et al., 2012) or gene activations (Smolen et al., 2000). Each link has a property called “weight” that is proportional to the strength of the interaction between any two members of the network. In our application, the weight is the TPM of each protein linked to a process; this is a directed network (i.e., protein-process), but a protein may take part of several processes. The framework is based on the background by Estrada-Peña et al. (2018) analyzing the metabolic processes of Anaplasma phagocytophilum growing in human cells.

The power of a network is its ability to obtain indexes that define the importance of each node, either protein(s) or process(es). There is wide variety of indexes that represent several features of the network allowing the evaluation of the main properties of either the network or the individual nodes. Our application measures the importance of each biological process throughout the feeding process or compares biological processes in different species of ticks. We selected the so-called PageRank (PR; Senanayake et al., 2015) that measures the connectivity of a node in reference to the nodes in the network to which it is attached. In other words, a process linked only to one highly expressed protein, will have a different PR than a process linked to many proteins slightly expressed. It is measurement of the “popularity” of the process. In other terms, PR indicates means how “cool” is the biological process as linked to proteins.

Figure 1 sketches the explanation above. There are proteins ( Figure 1A ) with different TPM’s linking to low level biological processes ( Figure 1B ). The width of the lines is proportional to the protein’s TPM ( Figure 1C ). The LL processes are the “children” of TL processes ( Figure 1D ). These LL or TL processes are provided by the hierarchy of GO, e.g., a LL process (like “activation of C-type lectin signaling”) is part of a TL process (in the example, “response to stress”). It may happen that some proteins are missing ( Figure 1F ) or their TPM’s are higher or lower; therefore, the PR of both LL and TL will change ( Figures 1G, H , schematized according to the width of lines and the size of the circles). The PR index is thus a way to measure the expression of a biological process. Figure 2 includes an actual network of the LL in Rs G1 feeding ticks, as well as some sections augmented to provide with better context to the readers. The purpose of the Figure 2 is merely demonstrative and is not intended for the visual observation of each individual node or link.

We generated separate networks for either the LL or TL in the case of the salivary glands of Rs for each period of engorgement, as well as for salivary glands and gut of Or. The PR was calculated for each network. All the networks were built using NodeXL (https://nodexl.com; last accessed, June 2020) and Gephi v0.92 (Bastian et al., 2009). The purpose was to compare the different PR of LL and TL in the feeding timeline of Rs, and to pinpoint the main differences with Or. We will refer to “unfed” or “G1” to “G6” the periods for Rs (according to the original definitions by Tirloni et al., 2020), and “OrSG” or “OrGut” to the corresponding salivary glands or gut of Or.

We wanted to estimate the variability of the PR of each LL or TL along the feeding period of Rs. We employed the percent variance of the PR as a reliable statistical analysis: the value of PR in each group gives the percent variance with respect to the past value along consecutive feeding groups. This procedure produced information for literally thousands of LL, which resulted inadequate for a graphical representation (see Supplementary Table 1 for the complete list of LL used in this study). Therefore, after the calculation of the percent variance between any two consecutive groups, we selected the LL that have a variance of change above the percentile 98 (i.e., the top 2% of higher variability) providing a total of 71 LL that display the maximum rates of change throughout the feeding period. As explained before, our aim is not to analytically examine every single process, but to provide with examples about the network application and its possibilities. We however included the complete set of TL and charted the variability of PR. Charts plotting values of PR were generated using the R environmental programming (R Core Team, 2020). The number of shared proteins or LL among feeding groups and species was compared using Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/, last accessed August 2021). The topology of Venn diagrams does not allow to use more than 5 different groups in one single chart; we therefore produced these charts separately into meaningful groups.

We aimed to demonstrate the usefulness of the method selecting one of the TL that did show changes along the feeding of Rs, namely “response to stress”. We chose every LL involved in “response to stress”, according to the hierarchy in GO, and produced a heatmap of the PR values of the top LL showing the maximum rates of change (about 30% of the total). We also obtained subnetworks involving only the selected LL’s. To note, this is not a purposely selected chain of processes aimed to demonstrate a particular event. It is intended only as a proof-of-concept of the reliability of PR to selectively track the changes of biological processes in the metabolism of the tick and to know the processes of the salivary glands via the relative composition of the sialome. It must be noted, too, that because we did a selection of proteins, driving to a selection of the processes under change, the results obtained cannot be interpreted as the complete sequence of changes through feeding, but a demonstration of the study flow.

The discussion of this study includes possible gaps that this approach could produce on the results obtained, mainly regarding the filtering proteins and the hypothesized problems derived from the original BLAST hit. With the actual knowledge and limitations in the field, we consider studies of this kind will move forward with the improvement of annotations of proteins in ticks.

Although this study is not devoted to the proteins identified in the sialome of both species of ticks, we think a background is necessary for further presentation of results. Some of these results have been already published in the original papers and are herein provided aiming to completeness. We describe first the general results about the number of proteins and LL in Rs, then we focus on the LL and TL of Rs. We will briefly describe the results for Or and compare between Rs; last, we demonstrate the framework with an example made on data about “response to stress” in Rs. As mentioned later, some caution should be exerted reading the Results, mainly regarding to the preliminary selection of proteins of the study.

After cleaning of the datasets of proteins, we obtained a set of 1166 proteins in Rs (all groups) and 2150 in Or (SG and Gut). Some of these proteins are shared among the different groups of Rs and with Or, but in different proportions regarding the species and the feeding groups ( Figure 3 ). In the case of Rs, there is a set of core proteins (767-771, Figures 3A, B ) that are always present in all the groups. It is important to remark some interesting associations: (i) the high number of shared proteins between G1-G2 (406, Figure 1A ) and G4-G5 (410, Figure 1B ) but not G6 (36, Figure 1B ) (ii) the fewer shared proteins between Or and the groups of Rs (299, Figure 1C ) (iii) the high similarity of the proteins detected in both SG and guts of Or (2083, Figure 1C ); only 31 proteins have been detected solely in SG of Or, and 40 in its gut. In summary ( Figure 1D ) unfed Rs lack 464 proteins already found (according to our filtering rules) in other feeding pools; groups G3 and G6 of Rs lack 430 and 432 proteins, respectively, that were detected in the other groups. The maximum number of non-shared proteins among the groups of Rs is six (in a total of 1246). Collectively, Rs and Or do not share the 91% of the total proteins recorded in both species (but read “Discussion” about the possible gaps in the design of the calculations).

The queries to Uniprot and GO produced a variable number of LL ( Figure 4 ). There is a total of 3121 LL in Rs, considering every feeding group; Or has 2535 LL in the SG and 2541 in the gut (2531 of them are shared). The core group of LL processes shared by Rs is 2172-2176. However, when compared against Or, the number of shared LL drops to 1980 comparing the unfed group with fed Or. Also, this number continued to fall to 535 when comparing a fed group of Rs and feeding Or. This is a striking result, since the lack of 91% of shared proteins between Rs and Or drives the lack of only the 10% of LL biological processes. These results should be observed with caution. Up to 902 LL (collectively detected in the feeding groups) are missing in unfed Rs. Of interest, both G3 and G6 feeding Rs lack 812 and 849 LL.

The LL represent the most fundamental bricks of the biological processes of the feeding tick. A chart of the changes of PR of each LL through the feeding sequence of Rs would be hardly readable. We thus chose the top 2% of most variable LL’s charting the log10 of their PR values. The heatmap in Figure 5 includes a dendrogram grouping the LL that tend to have similar values of PR along the tick feeding. The analysis of the LL processes in the top 2% reveals important differences in the functional biology of groups of feeding ticks. These LL can be grouped into three categories: those that have a relatively low PR (but in G2), those that keep a high PR throughout the feeding process, like “translation” (but in G2) and those that show intermediate to high values of PR, with slight variations through the whole feeding (but in G2). To note, the heatmap shows only the top 2% of LL with highest rates of change, obtained from an already filtered subset of proteins of the total proteome.

Without questions, G2 represents a special breakthrough point in the feeding process of Rs at the LL observational level, as it could be expected for ticks beginning the feeding process. The analysis shows that some LL are suddenly activated at G2 (breakpoint “A” in Figure 5 ), while some others decrease its PR, even if they keep the highest values in other groups (like “translation”, breakpoint “B” in Figure 4 ); others are “arrested” specifically at G2, even if they are expressed at medium (breakpoint “C” in Figure 5 ) or higher values (breakpoint “D” in Figure 5 ) in the rest of feeding groups. As mentioned, the figure is intended to show the analytical abilities of a network approach to extract conclusions using only the TPM of the detected proteins.

The analysis of TL processes narrates a different history of the feeding Rs. We identified 47 TL processes in unfed and feeding Rs. The networks produced for the TL for all the groups are provided as part of the Supplementary Material Figures (S1–S7) . Interestingly, some TL processes remained undetected at some time slices of the feeding of Rs. The analysis of TL clearly states that major changes in the sialome of Rs are obviously produced at the beginning of feeding (G1) in which many TL increase its activity. This finding is not contradictory with the previous comments about LL. We consider that the analysis of TL is a “panoramic view” of the biological processes of the SG of the tick, while LL are the bricks of the chain of events of the tick feeding. The number of LL per TL is very variable. Some TL are represented by only a few LL (or even only one). This immediately translates into the very variable number of proteins that are involved in the set of LL.

The Figure 6 charts the log10 of PR of every TL in the consecutive feeding groups of Rs. This figure re-arranges them to produce a dendrogram of the TL that change its PR in similar patterns along the feeding. The beginning of feeding (G1) marks the expected breakpoint in comparison with the unfed ticks. The pool of G1 have the most notorious changes of PR in more than 80% of TL processes. Some TL experience a decrease of PR from unfed to G1, and then continue to be highly represented in ongoing groups (breakpoint A in Figure 6 ). Exceptions are “reproductive process”, “biogenesis”, and “developmental processes”, all of which have a low PR that increases in G1 and continue at high values until the termination of feeding (breakpoint B in Figure 6 ). To note that the TL associated to reproduction show a clear increase as soon as G1, and not at the end of the of feeding: the morphological and physiological changes associated with the reproduction start as soon as the tick begins to feed. Some other TL experience a relative increase in G1 to decrease again in the subsequent groups (breakpoint C in Figure 6 ).

In terms of absolute rates of change of PR between unfed and G1, the PR of “cell proliferation” increased more than 65,000 times, “localization” (a TL involving the cellular organelles) increased more than 13,000 times and the “immune system processes” increased almost 2,000 times. All these TL returned to levels comparable with unfed ticks after G2. However, ticks in G2 group experienced dramatic changes, mainly those related to “metabolic processes” (with a PR 54,000 times higher). Most other TL processes have lower changes through the feeding. The changes occurring in the G6 group that includes ticks detached after feeding are of interest. In example, the “metabolic process” that remained at relatively low levels through the previous phases of the feeding period (but for G1) increase again more than 9,000 times, overlapping with the final period of intense blood intake.

The changes observed in LL have not a direct translation to TL; they are not directly correlated, since i.e., most prominent changes of LL were noticed in G2. We hypothesize that there are two levels of observation: the LL allow to capture a finely tuned description of the actions of tick sialome, while TL only allow a wide (but reliable) description of changes. Anyway, even if TL alone may be a rough way to define the physiology of the feeding ticks, regarding the fine detail obtained examining LL, the changes in tick sialome through feeding are obvious and have a biological support.

This study is mostly focused to capture the dynamics of the changes of sialome in the feeding R. sanguineus. However, we aimed to compare data for two species of the major families of ticks, using data of Or as an example. To note, the dataset of Or only reported one slice of time, and is not extended into the process of digestion. Also, the dataset of Rs lacks data on gut, so the comparison is limited. As before, we did chart the LL processes that have significantly different log10 PR in SG and gut ( Figure 7A ) and the complete set of TL processes ( Figure 7B ). Of interest, most variable LL have higher PR in the gut than in SG. Variability in the remaining LL (more than 2000 LL) is lower. Figure 7B displays the log10 of PR of the complete set of TL in both organs. Differences are minimal and there is an almost pair-to-pair correspondence of the PR values of TL in the examined organs of Or. The networks produced for the TL for both organs are provided as part of the Supplementary material Figures (S8, S9). Of interest, metabolism related TL have far smaller PR in Or than in Rs, while the cellular processes and a wide array of TL involved in several types of regulation are more abundant.

However, these results, together with those comparing the number of shared proteins in both families of ticks are striking. These results could be derived from an unreliable selection of proteins, since they were extracted and queried from different proteomes. We will further discuss about these findings and the possible sources of error.

We focused on the analysis of the chain of events involved in the TL “response to stress” in unfed and feeding Rs. The TL includes a total of 226 proteins and 149 LL. As before, it is not possible to produce a readable chart for every LL; we therefore opted for focusing the analysis on the LL with highest rate of change across the feeding, resulting in the 30% of LL. The Figure 8 includes a heatmap charting the log10 transformed values of PR throughout the feeding period of Rs including only the LL for “response to stress”. The chart includes a dendrogram grouping the LL that tend to change the PR values in a similar pattern. A major increase of PR is observed for tick pools collected at G2, G4 and G5 (breakpoint A in Figure 8 ). They share the large increase of PR, in LL, like the C-type lectin receptor signaling pathway, the negative regulation of inflammatory response and the blood coagulation. There are LL that increase its PR but at are almost absent at some periods of feeding like i.e., the negative regulation of oxidative stress, the acclimation of cells to heat, and the defense response to viruses (breakpoints B, C in Figure 8 ). Other LL ( Figure 8D ) have relatively medium to low values of PR throughout the complete feeding period, showing some increase mainly at G2, G4 and G5 (i.e., response to cold, or hyperosmotic response) or G1 (processes involved in double strand DNA repair). The rest of the processes ( Figure 8E ) are predominantly common in unfed, G1, G3 and G6, with a small increase in the rest of the feeding phases. Other LL involved in “response to stress” are well below the percentile 70 of rate of change, and they exhibit very small changes through the complete feeding.

The subnetworks obtained for that specific chain of biological process are also included in Figure 8 , only for some slices of the feeding. The topology of these networks, consisting in a relatively dense nucleus with many dots in the periphery explains the timeline of the process. Unfed ticks have a relaxed structure of proteins, to move into a cohesively linked network of LL in ticks belonging to G1. However, G3 and G6 show a large group of dots tightly located in the center of the network, denoting large interactions among them, while some of the proteins and processes do not belong to that nucleus of the network and are located in the periphery, playing secondary actions.