Understanding the evolutionary structural variability and target specificity of tick salivary Kunitz peptides using next generation transcriptome data

Background Ticks are blood-sucking arthropods and a primary function of tick salivary proteins is to counteract the host’s immune response. Tick salivary Kunitz-domain proteins perform multiple functions within the feeding lesion and have been classified as venoms; thereby, constituting them as one of the important elements in the arms race with the host. The two main mechanisms advocated to explain the functional heterogeneity of tick salivary Kunitz-domain proteins are gene sharing and gene duplication. Both do not, however, elucidate the evolution of the Kunitz family in ticks from a structural dynamic point of view. The Red Queen hypothesis offers a fruitful theoretical framework to give a dynamic explanation for host-parasite interactions. Using the recent salivary gland Ixodes ricinus transcriptome we analyze, for the first time, single Kunitz-domain encoding transcripts by means of computational, structural bioinformatics and phylogenetic approaches to improve our understanding of the structural evolution of this important multigenic protein family. Results Organizing the I. ricinus single Kunitz-domain peptides based on their cysteine motif allowed us to specify a putative target and to relate this target specificity to Illumina transcript reads during tick feeding. We observe that several of these Kunitz peptide groups vary in their translated amino acid sequence, secondary structure, antigenicity, and intrinsic disorder, and that the majority of these groups are subject to a purifying (negative) selection. We finalize by describing the evolution and emergence of these Kunitz peptides. The overall interpretation of our analyses discloses a rapidly emerging Kunitz group with a distinct disulfide bond pattern from the I. ricinus salivary gland transcriptome. Conclusions We propose a model to explain the structural and functional evolution of tick salivary Kunitz peptides that we call target-oriented evolution. Our study reveals that combining analytical approaches (transcriptomes, computational, bioinformatics and phylogenetics) improves our understanding of the biological functions of important salivary gland mediators during tick feeding.


Background
Ticks have been successful as ectoparasites due to morphological and physiological adaptations [1]. As obligate hematophagous (blood-sucking) arthropods, ticks transmit various bacterial and viral diseases, e.g., babesiosis, theileriosis, anaplasmosis, Lyme disease and tick-borne encephalitis, and, thus greatly impact human and animal health [2]. To combat host defense mechanisms, ticks possess powerful pharmacological proteins in their saliva that are injected into the vertebrate host during blood feeding [3,4]. Our current understanding of host-parasite interactions has been greatly impacted by the revolution in sequencing technologies that has provided a massive amount of data previously unimaginable [5]. Among the various transcriptome studies to date, transcriptome projects from arthropod disease vectors are using Sanger or next generation sequencing techniques. These include several salivary gland (SG) transcriptome and proteome projects (also known as sialomes) from hard tick [6][7][8][9][10][11][12][13][14][15][16][17][18] and soft tick species [19][20][21][22]. Hundreds of novel SG transcripts that putatively encode proteins were discovered from these sialome projects, thus elucidating how ticks may complete a blood meal while providing insight towards their evolutionary expansion [23].
Kunitz SG peptides of I. scapularis were phylogenetically analyzed to uncover their emergence in ticks and the expression trends of these Kunitz peptides were also statistically analyzed [23]. These Kunitz peptides were categorized in three different groups (groups I, II and III) based on their Cys motif. The I. scapularis Kunitz peptides belonging to group I were suggested to represent the ancestor of all tick Kunitz-domain family (single and multiple domains). Several Kunitz peptides seemed to have lost their ability to function as serine protease inhibitors and instead to block and/or modulate ion channels, possibly related to the tick's necessity for prolonged feeding on the vertebrate host [23]. The authors are aware of only one study that functionally and structurally characterized a tick Kunitz peptide as an ion channel effector, the maxiK channel modulator Ra-KLP from R. appendiculatus [40]. This paucity must be remedied since Dai et al. [23] putatively characterize the majority of I. scapularis Kunitz peptides as ion channel blockers/modulators. Furthermore, Fry et al. [41] have argued that hematophagous secreted proteins, such as of the Kunitz family, should be classified as venomous. Most classified toxins are stabilized by their disulfide bridges and once these toxins become functionally essential as a venom, their adaptation is often reinforced by gene duplication [41]. Gene sharing and gene duplication are the main mechanisms advocated to explain the functional heterogeneity of tick salivary Kunitz family proteins [23,42].
In our study, we used computational, structural bioinformatics and phylogenetic methods to reevaluate tick salivary Kunitz peptides from a more in-depth structural point of view by analyzing the functional, antigenic, and evolutionary characteristics of Kunitz peptides from the recently annotated Ixodes ricinus SG transcriptome [18]; GenBank Bioproject PRJNA177622. Compared to classical biochemical analyses and classical Sanger sequencing techniques that revealed only a few thousands of sequences from tick transcriptome studies presented until today (with the exception of the A. maculatum transcriptome), the large amount of available data of the 454 pyrosequencing/Illumina I. ricinus SG transcriptome makes it feasible to thoroughly analyze multigenic protein families (i.e., the Kunitz-domain family). In the I. ricinus transcriptome, 4% of all Illumina reads were classified as transcripts encoding for Kunitz-domain proteins and of these, 1.4% accounted for single Kunitzdomain transcripts (that putatively encode archetypal Kunitz peptides possessing more or less 6 Cys residues) [18]. Our approach demonstrates how to interpret next generation transcriptome data to expand our understanding of the molecular, structural and functional evolution of hematophagy in ticks. Our overall analysis developed a novel concept that elucidates the divergence and adaptation of I. ricinus salivary Kunitz peptides that we call target-oriented evolution.

Results and discussion
Organizing the different I. ricinus SG Kunitz peptides based on their Cys motifs facilitates functional predictions From the 454/Illumina SG I. ricinus transcriptome [18], we found archetypal Kunitz peptides containing 6 Cys residues and variants containing an additional 7 th Cys residue adjacent to the conserved 5 th Cys residue ( Figure 1). Based on their Cys motif we were able to organize these tick salivary Kunitz peptides into 11 groups (G1-G11; Figure 1A) Kunitz groups based on Cys motif. Kunitz peptides from the annotated transcriptome of I. ricinus were organized into 11 groups (G1-G11) based on their Cys motif (C = Cys and X = intra-Cys residues). The Cys motif for each group is depicted in Panel A along with their predicted disulfide bond patterns (in brackets). The alignment of one representative for each group (B) shows the conserved Cys residues (red), the positions of aa variability (yellow), and the key residues that interact with serine proteases (blue) and ion channels (green). The box around G3, G6 and G11 indicates that they do not contain all of the standard residues for serine protease (blue) or ion channel (green) interactions. Two separate alignments in Panel C depict the key conserved residues for serine protease inhibition (blue) using BPTI (UniProt: P00974) and ion channel blockage (green) using beta-bungarotoxin (UniProt: P00989). Residues highlighted in black are conserved. and a group with ≤4 Cys in their amino acid (aa) sequenceswe will refer to these peptides as simple Kunitz (SK), hereafter. The Pfam server (http://pfam. sanger.ac.uk/) verified that all 11 groups plus the SKs possess the conserved single Kunitz-domain. Group G6 is the only I. ricinus Kunitz peptide group that possesses a different disulfide bond pattern since it lacks the archetypal Kunitz/BPTI Cys2 and Cys4 disulfide bridge ( Figure 1A). The 3D modeled structures of all group representatives from Figure 1 show that the overall archetypal tertiary Kunitz fold is highly conserved regardless of the deviated disulfide bridge pattern for G6 (Additional file 1A-B). Furthermore, the tertiary Kunitz fold for each group representative does not structurally deviate much from each other (Additional file 1C and Additional file 2). The deviation from the archetypal Cys motif (as in G6), however, has been shown to be functionally flexible [26,32,40].
Aligning a representative of each I. ricinus Kunitz peptide group ( Figure 1B) in our study shows that specific residues coincide with the putative functions described by Dai et al. [23], classified as serine protease inhibitors or ion channel blockers/modulators. Groups G3, G6 and G11, however, do not contain these conserved residues (boxed in Figure 1B). With the exception of the annotated I. scapularis sialome [6,9], both protein BLAST [43] and PSI-BLAST [44] searches do not reveal any known or functionally described single Kunitz-domain peptide with the Cys motif of G3, G6 and G11. The proteins that slightly resemble the Cys motif of G6, and lack the archetypal Cys2 and Cys4 disulfide bridge, are the I. scapularis salivary Kunitz proteins Ixolaris [38] and Penthalaris [39], both possessing multiple Kunitz domains, and some members of the Kunitz family found in R. microplus [45]. An alignment using the classical Kunitz serine protease inhibitor BPTI (UniProt: P00974) with G1 and G7 shows the conserved positively charged arginine (Arg or R) or lysine (Lys or K), known as the P1 site that interacts with the negatively charged active site of serine proteases (highlighted blue in Figure 1C). Furthermore, the key residues for the ion channel blocker beta-bungarotoxin (UniProt: P00989), from the venomous snake Bungarus multicinctus [46], are also conserved in G2, G4, G5 and G8-G10 (highlighted green in Figure 1C).
By organizing the Kunitz peptides from the I. ricinus transcriptome into specific groups, we were also able to visualize specific profiles occurring in the transcript reads among the four independent cDNA libraries, suggesting that there may be a correlation between the average transcript read value and putative function for each group. The three groups with the highest number of transcripts were G1, G2 and G6 (G6 being the largest group). Figures 2A-C depicts the average Illumina reads for all groups obtained from the I. ricinus transcriptome, complementing our functional predictions. We notice that three profiles are occurring among the organized transcripts encoding for Kunitz peptides. Figure 2A (G1 and G7) shows that the average transcript read level (Kunitz encoding transcripts that we classify as putative protease inhibitors) is low at the early stages of tick feeding for both nymphs (EN: 6 h to 24 h of feeding on the host) and adults (EA: 6 h to 2 days of feeding on the host) compared with the increased transcript reads during late feeding (LN: 2-4 days of feeding and LA: 3-7 days of feeding). Soft ticks have a shorter feeding cycle than hard ticks and this phenomenon has been thoroughly explained by comparing the evolutionary divergence and emergence of serine protease inhibitors in hard ticks and soft ticks (about 120 and 92 MYA) [42]. Giving that hard ticks (like I. ricinus) feed on blood for longer periods of time (females: ≥1 week) and that serine proteases are involved in blood coagulation and immune responses, we expect a higher transcript read encoding for serine protease inhibitors during late feeding stages as depicted in Figure 2A (for both nymphal and adult life cycles).
We putatively classify G2, G4, G5 and G8-G10 as ion channel blockers and/or modulators ( Figure 1C) and as Figure 2B shows they exhibit a separate profile than the serine protease inhibitors of Figure 2A (G1 and G7). These putative ion channel blockers/modulators show lower transcript reads during the nymphal stages of feeding (both early, EN, and late, LN), but the reads are increased in adult ticks during the early stage of blood feeding (EA). The profile then decreases during the late stages of feeding ( Figure 2B). A possibility, as noted in mayfly nymphs SGs, for the expression of ion channel blockers/modulators in I. ricinus nymphs may partially be due to maintaining osmoregulatory function during molting [47]. To resist desiccation, tick nymphs rely on both water vapor and water retention, while adult ticks only rely on enhanced water retention [48]; therefore, this may explain the increased transcript read encoding for ion channel blockers/ modulators during adult feeding stages (compared with nymphs - Figure 2B). In Figure 2C we show that G6 deviates from the other two profiles, however, G3 and G11 have similar profiles to that of ion channel blockers/modulators ( Figure 2B). The distinct profile and Cys motif G6 displays suggests that this group may belong to (a) a combination of serine protease inhibitors and ion channel blockers/modulators, (b) either serine protease inhibitors or ion channel blockers/modulators, or, (c) a completely new, undefined functioni.e., gene sharing (moonlighting).
The I. ricinus SG Kunitz groups possess intra-Cys residue variability, assemble into various secondary structure clusters, and vary in antigenicity and protein disorder The aa positions with high variability (highlighted yellow in Figure 1B) for the intra-Cys residues of each group were calculated by the Protein Variability Server [49,50] using the Shannon entropy equation (see Methods) [51] are represented in Figures 3A and D as the ratio of highly variable/highly conserved sites (V/C) and the average of aa variability (AVE) per group, respectively. On one hand, we consider G3, G5, G7, G10 and G11 are "frozen" groups having low aa variability and high proportion of highly conserved sites. On the other hand, G1, G2 and G6 are "melted" groups or highly variable. Due to their low aa variability, the former assemblage may be considered highly adapted in the host-parasite interaction from a functional and/or immunological point of view. The second assemblage may still be under evolutive pressures that force aa variability.
To demonstrate the structural outcome of these "frozen" and "melted" clusters, we predicted the secondary structure for each Kunitz peptide in all groups. We show that 22 models depict the variability of secondary structure among the 145 Kunitz peptides (excluding SK peptides) (Additional file 3). Five models in Additional file 3 (Representative models G1: Ir2:98; G1: Ir2:4878; G1: Ir2:19262; G6: Ir2:1982 and G7: Ir2: 4792) are combined β-strand and α-helical secondary structure models, differing in length, amount and order of α-helices and β-strand, and describe 81.6% of the secondary structure variation present in the I. ricinus Kunitz/BPTI family. All the groups contain different amounts of secondary structure models ( Figure 3B and Additional file 3). It seems, however, that the secondary structure presented in Ir2-4878 and Ir2-98 of G1 is the structural skeleton for the remaining secondary structures; these structures are also present in 10 of the 11 groups and explain 50% of the secondary structures. The fact that G1 represents the structural skeleton may be since its Cys motif (as displayed in Figure 1A) is considered the ancestor of Kunitz peptides in I. scapularis [20]. Structural convergence, however, has been reported for venoms, a general classification that the Kunitz protein family belongs to [Revised in 41].
When analyzing the aa variability (V/C or AVE) in relation to the secondary structure models per group (excluding G6) a high correlation is given (R 2 = 0.78; Figure 3E); however, by including G6 in this analysis the correlation decreases (R 2 = 0.18; Figure 3F). This dictates that the aa variability per group is partly justifiable for structural reasons (functionality), except in G6. Figure 3C shows that ratio V/C or AVE of aa variability are correlated, therefore, using V/C or AVE of aa variability is interchangeable for Figure 3E and F. The presence of different folds (i.e., secondary structures) suggests that gene sharing events may have occurred throughout the evolution of I. ricinus and may cause additional and unique functional properties. Such moonlighting proteins may explain the diverse inhibitory effect of tick salivary Kunitz peptides.
To elucidate additional molecular properties, we analyzed both the antigenic properties and protein disorder for all I. ricinus Kunitz peptide groups in our study. Our analysis reveals that the regional disorder was located at both termini for the majority of peptides (data not shown). It is worth noting that antigenic properties involve the mobility of the C-and N-termini [52]. In Figure 4A we show that G6 has the highest antigenicity mean (VaxiJen score). G6 has a "misplaced" Cys2 (archetypal) and lack of this Cys2 is responsible for the functional flexibility of the Kunitz-domain binding loop, that may permit target-specific interactions other than the catalytic domain of serine proteases [53]. We also notice that groups with only two sequences vary in antigenicity ( Figure 4A) that also coincides with differences in secondary structure (Additional file 3). Functional flexibility (or mobility) is also demonstrated in G6, since it has a higher variability in intrinsic protein disorder ( Figure 4B). Intrinsically disordered regions (or proteins) increase molecular recognition because of an ability to fold differently upon binding and possess large interacting surfaces [54]. Interestingly, G6 has one sequence that is completely disordered (Ir2-10518) that is also highly antigenic (>0.7).

The I. ricinus SG Kunitz transcripts are under different selective pressures
The Red Queen hypothesis was proposed by Van Valen in 1973 [57] to explain the molecular evolution through evolutive interactions of prey-predator or host-parasite. The hypothesis establishes "that a proportional amount of successful response by one species produces a total negative effect on other species[…] To maintain itself as before, the [other] species must increase its fitness [in a proportional amount]". The I. ricinus SG Kunitz encoding transcripts may adhere to the Red Queen hypothesis, where inverse selection (positive or negative) is a series of one species adapting (the tick) and the other counteradapting (the host). To elucidate or validate the type of selective pressure (positive or purifying) each I. ricinus Kunitz peptide group has undergone throughout its evolution we used the Synonymous Non-synonymous Analysis Program (SNAP) (http://www.hiv.lanl.gov) and Datamonkey servers [58] (see Methods for further details).
Assuming that a ratio >1 between the non-synonymous nucleotide substitutions per site (d N ) with the synonymous substitutions (d S ) is due to a positive selection and a ratio <1 is due to a purifying selection, we note in Table 1 that the majority of groups are under a purifying selection. As we depict in Figure 4C, only G4 had a strong Figure 3 Graphical representations of amino acid and secondary structure variability. The relation between highly variable and low variable sites (V/C) among the groups is shown (A) The number of secondary structure models per groups is also presented (B). Correlation between the average of aa variability (AVE) and the ratio V/C (C). Average of aa variability estimated by Shannon entropy (D) and the correlation between the ratio V/C and 2D models, without and with group 6 (E-F).
correlation between disorder and antigenicity; coincidentally, G4 is the only group that has undergone positive selection (Table 1). This correlation ( Figure 4C) suggests that the disorder in G4 may be under positive selection due to the immune system of the host, since protein mobility (disorder) may influence the antigenic properties of proteins [52]. Our SNAP analysis, however, only displays the average of synonymous and the non-synonymous nucleotide substitutions per site, therefore, is not to say that selection per site was not positive or negative for specific members of the remaining groups. We therefore used Datamonkey to verify our SNAP analyses for evidence of positive selection and to determine which site along the nucleotide sequence are undergoing positive (natural) or purifying (negative) selection.
Overall, the ratio d N /d S from Datamonkey slightly differ from that of SNAP. The Datamonkey analyses show that G2, G4 and G6 possess evidence for positive selection and note the type of selection for each site, using the SLAC algorithm (see Methods). Figure 5 graphically displays these results along the entire transcript, codonby-codon (x-axis). Both G5 and G10 do not possess any sites under any selective pressure, while the remaining groups possess positive, negative (purifying) or both (see Table 1). As depicted in Figure 5, the majority of selected sites are negative (purifying) and groups G2 and G6 possess a few positively selected sites; however, there are more negatively selected sites for G2 and G6 than positively selected sites. In a recently published book, Eugene Koonin [61] explains that purifying (negative) selection, in some phases of evolution, is more common (orders of magnitude more common) than positive selection. Koonin considers purifying selection is the default process of elimination of the unfit. With this regard we understand that Kunitz-domain proteinstaking into account the number of members in generalshould undergo a massive purifying selection in order to shape (constrain) the molecular diversity of this tick SG Figure 4 Antigenicity and intrinsic protein disorder. Variations in antigenicity (A) for each group as predicted by the VaxiJen server [55] that uses an alignment-independent algorithm based on the physicochemical properties of proteins [55]. Ratio of disordered compared to the total sequence (B) for each group as predicted by the MetaDisorder server that builds a consensus from 13 protein disorder predictors [56]. Panel C shows the correlation for G4 between antigenicity (x-axis) and intrinsic protein disorder (y-axis).
peptide family. Additionally, as noted by the reduced gene duplications events in platypus venom compared with the massive gene expansions reported in cone snails and spiders [62], selective pressures may fluctuate according to venom function.

The evolution of I. ricinus SG Kunitz peptide groups
The phylogeny of the I. ricinus Kunitz peptides from our study was reconstructed and divergence times were estimated. Maximum likelihood (ML) and Bayesian phylogenetic methods resulted in similar tree topologies ( Figure 6 and Additional file 4A). Four main clades containing 9 out of 11 I. ricinus Kunitz peptide groups were well supported in the ML (Additional file 4A: bootstrap values ≥72, marked in red) and Bayesian trees (Additional file 4B: posterior probabilities ≥ 0.7, marked in red). The largest clade out of these four clades was composed of only Kunitz peptides of prostriate ticks, namely from I. scapularis and I. ricinus groups G2, G4-5, G8-10 that we characterize as ion channel blockers/modulators. Group G11 also appeared in this clade and, regarding its transcript read profile ( Figure 2C), G11 may also function as an ion channel blocker/modulator like the other Kunitz members of that clade. Although the maxiK channel modulator Ra-KLP from R. appendiculatus did not group with the "ion channel blocker/modulator" clade, prostriate tick salivary proteins may completely function differently than metastriate tick proteins. Thus, our characterization of G2, G4-5, G8-11 as ion channel Figure 5 Positively and negatively selected sites of I. ricinus Kunitz transcripts. Codon-by-codon (x-axis) of the nucleotide substitution selected sites (positive or negative; y-axis) were calculated by the Datamonkey server [58] via the SLAC algorithm [59]. The graph depicts the groups that had evidence (p-value < 0.05) for selection as indicated in Table 1. Of note, G3, G7 and G11 (and simple Kunitz) were not included since they contained only two sequences. blockers/modulators may still be valid and future experimental verifications will clarify this uncertainty. The second largest clade also contained only prostriate Kunitz peptides composed mainly of I. ricinus group G6. As in the ion channel blocker/modulator clade, no functionally described Kunitz appeared in this clade, therefore any possible function remains unknown. Functionally described Kunitz peptides were present in the remaining two largest clades. One of these clades was made up of only argasid Kunitz peptides that included anti-platelet inhibitors. This group seem to evolved independently from hard ticks as previously postulated by Mans et al. [42], possibly due to different blood feeding behavior compared with hard ticks. The other clade mainly contained I. Figure 6 Molecular clock analysis of I. ricinus Kunitz peptides. Kunitz nucleotide sequences of I. ricinus groups (G1-G11 and SK) and of ticks and spiders (defined by GenBank accession numbers) were used to infer divergence times within ticks by a Bayesian uncorrelated relaxed lognormal molecular clock model. Four taxon sets were calibrated according to Jeyaprakash and Hoy [63] and the nodes with the recovered divergence times were labeled: Araneae/Scorpions/Pycnogonida/Acari 459±18 MYA (node S), Argasidae 214±28 MYA (node A), Prostriata 196±27 MYA (node P) and Metastriata 134±22 MYA (node M). The figure shows a maximum clade credibility tree including the mean divergence times (MYA) at each node. Well supported nodes and matching ML and Bayesian analyses (bootstrap values ≥70%, posterior probabilities ≥0.7, see Additional file 4) are indicated with black circles. Nodes with gray circles appeared in the ML and Bayesian trees, but they were not well supported (bootstrap values ≤70% and posterior probabilities ≤0.7). Black squares indicate well supported nodes by the Bayesian analysis (posterior probabilities ≥0.7), but were not resolved in the ML tree. Gray squares indicate well supported Bayesian tree nodes (posterior probabilities ≥0.7) that were less supported in the ML analysis (bootstrap values ≤70%). Non-labeled nodes were not well supported in both analyses (bootstrap values ≤70% and posterior probabilities ≤0.7) and they did not match in both phylogenetic reconstructions. The four best-supported clades (supported well by both analytical methods) are labeled with red circles at the nodes. The red star at the node highlights a clade containing G8 and G11 members at a divergence time of 71 MYA that is similar to the divergence time of G6. The scale bar in MYA is at the bottom center of the tree.
scapularis peptides together with I. ricinus group G3 and anti-clotting inhibitors of Ornithodoros spp.. Our aa sequence and structural analyses of G3 did not reveal any putative function for this I. ricinus group, although the transcript profile of G3 was similar to that of ion channel blockers/modulators (as G11 in Figure 2C). However, since G3 clusters with anti-clotting inhibitors and not with G11 may suggest a different function during tick feeding. I. ricinus groups G1 and G7 clustered together with hard and soft tick Kunitz peptides, but their relationships are uncertain since they were not well supported in the ML and Bayesian trees (Additional file 4). Nevertheless, in our Bayesian analysis, I. ricinus groups G1 and G7 appeared with functionally described serine protease inhibitors and G7 proteins were also grouped together with G1 proteins in our ML analysis. Both our phylogenetic and structural analyses of the I. ricinus groups G1 and G7, suggests that these two groups may act as serine protease inhibitors. SK peptides appeared throughout the phylogram. Some SKs were similar in their aa sequence (Additional file 5) with the other I. ricinus Kunitz peptides assuming a similar function during blood feeding. Others SKs may be more functionally distinct, such as Ir2-24967-SK, since it differs in its aa sequence from all other I. ricinus Kunitz peptides.
Our molecular clock analysis of tick and spider Kunitz peptides recovered similar divergence times for the main splits of 1) spiders from ticks (426 MYA, 95% HPD: 389-462 MYA, node S), 2) soft ticks from hard ticks (233 MYA, 95% HPD: 198-266 MYA, node A) and 3) prostriate ticks within Ixodidae (221 MYA, 95% HPD: 181-254 MYA, node P) as calibrated for these taxon sets; the tree root was estimated as 436 MYA (95% HPD: 395-479 MYA) ( Figure 6, Additional file 4B). Most of these splits were well supported in the chronogram. Metastriate Kunitz peptides (fourth taxon set), however, started to evolve earlier than calibrated. Nevertheless, this split is not well supported in the chronogram as previously stated (193 MYA, 85% HPD: 160-228 MYA, Figure 3, node M). Although I. scapularis homologs of I. ricinus group G1 Kunitz peptides are considered to be the ancestral form of all Kunitz-domain proteins [23], our phylogenetic reconstructions of I. ricinus Kunitz peptides cannot confirm these findings. Even though G1 is not well supported in the ML and Bayesian analyses, it can be assumed that this group of proteins evolved convergently in ticks as both phylogenetic trees indicate (Additional file 4). Expansion of protein families due to gene duplication events has been previously reported for ticks [23,64] and can explain the convergent evolution of G1. Overall, the I. ricinus Kunitz peptides in our molecular clock analysis underlie recent gene duplication events that explains how several members of the I. ricinus groups and other tick Kunitz peptides do not strictly follow the calibrated divergence times within the taxon sets in the phylogenetic trees.

Target-oriented evolution, a model for I. ricinus SG Kunitz peptide evolution
The Red Queen hypothesis [57] explains that hostparasites interactions must be characterized by constant mutations in both systems in order for parasites to efficiently infect the host and for the host to survive or avoid the parasitic attack. We understand the evolution of I. ricinus salivary Kunitz family in this dynamic framework. In the following paragraphs we explain our basis to propose one model of target-oriented evolution for this family of proteins and briefly discuss our findings regarding the evolution of venomous proteins.
Our model considers G6, compared with all other I. ricinus Kunitz peptide groups, as a cornerstone in understanding the Kunitz family evolution of I. ricinus. Firstly, members of G6 show unique molecular, structural and evolutive properties by possessing the highest amount of intra-Cys residues between Cys2 and Cys4 ( Figure 1), a distinct Illumina read profile ( Figure 2C), a high aa variability ( Figure 3A and D), a higher number of different secondary structures ( Figure 3B), and an increased number of negatively selected sites ( Figure 5 and Table 1). Secondly, in the phylogenetic analysis, the clade formed by G6 is mainly composed of I. ricinus Kunitz peptides (Figure 6 and Additional file 4), also representing the largest monophyletic clade suggesting a recent and specific evolution for this group in I. ricinus. Thirdly, based on the molecular clock analysis G6 is the fastest evolving group, diverging about 70 MYA ago and having the highest number of members. Taking into account the time of divergence and the number of members, we should also consider G6 as an example of possessing an accelerated rate of evolution that has been reported for several other families of venomous proteins [5]. At this point we consider it important to make a distinction. As depicted in the phylogenetic tree ( Figure 6) members of other groups appear to evolve around the same time point as G6. It is crucial to note here that even when some genes show the same time of evolutive divergence, as they belong to a distinct group, they do have different functional states at the moment of divergence. For example, our tree depicts a clade containing members of G8 and G11 that diverged 71 MYA ago ( Figure 6, marked with a red star), but, as we previously showed (Figures 1, 2, 3B, 4A-B), these groups have different properties compared with G6. In this way we consider that molecular mechanisms working for the observed Kunitz family expansion (e.g., gene duplication) must be operating on the whole genome background, thus, functionally and structurally different gene(s) may show genetic amplification at similar evolutive times. But we consider that the genetic expansion must be occurring at different rates for different groups of genes, depending on the function of the gene and its importance in the interaction with the hostas previously reported for venomous protein families [62].
We consider that G6 is an important proof that tick salivary Kunitz peptides possess a target-orientated evolution. By target-orientated evolution we mean the molecular properties that arise in the parasitic proteins through its interactions with the host targets and the host immune system that eventually determine its structural state. Our model (Figure 7) takes into consideration that different groups are undergoing different evolutive pressures (Table  1). This may explain, from a structural point of view, the evolution of the various functions found in the Kunitz family. In the model, we propose three main evolutionary categories: "working", "short-term" and "long-term" (Figure 7). The "working" group (i.e., G6) contains an increased number of highly variable members and will eventually reduce its number of members and gain in molecular specialization due to immunological and functional constrictions. Thus, becoming a variable but less flexible, "short-term" evolving group (i.e., G1, G2, G8 and G9). The high variability found in G1 and G2 compared with G6 may be interpreted differently. Non-monophyletic Kunitz groups, like G1 and G2 may be remnants of variable-monophyletic groups, like G6, that are evolutionarily static in expressed traits while accumulating genetic changes [65]. Finally, evolution will lead to a "long-term" evolved group (i.e., G3, G4, G5, G7, G10 and G11) that is less variable and highly specific in their molecular functions. Immune and functional pressures will remain during the whole process of both differentiation and maintenance of a gained specificity. As we showed in Table 1, I. ricinus SG Kunitz groups are evolving under different selective pressures and we predict they may have different functions (Figure 1 and Figure 2). Genes with different functions in venomous proteins families have been shown to evolve under different selective pressures [5] and thus at a different rate of evolution. In this sense we do not claim the three categories of "working", "short-term" and "long-term" as relatively exclusive to time, but we also refer them to a functional state (see above the example of G6 in relation to G8 and G11).

Conclusions
The Red Queen hypothesis is an arms race proposed as a model for co-evolution for host-parasite interactions [57]. Today's increasing knowledge about the interface between host-tick interactions has provided insight that ticks successfully feed to repletion by counteracting the host immune response, platelet aggregation, inflammation and vasoconstriction [4,64,66,67]. We therefore infer that two main evolutionary pressures may be driving the evolution in I. ricinus Kunitz peptidesi.e., immune response and the necessity of functional diversity against the diverse host molecular targets. The former is in agreement with specific immune response against tick salivary proteins [68], and, specifically against Kunitz members [69]. The pressure for functional diversity on tick Kunitz peptides may be that several targets have been associated with their inhibitory activity such as trypsin, elastase, kallikrein [70], chymotrypsin [71], tissue factor complex inhibitor [38] and ion channels [40]. We consider our model of target-oriented evolution (Figure 7) to explain the evolutionary dynamics found in the I. ricinus SG Kunitz peptides in the framework of the Red Queen hypothesis. Figure 7 Proposed mechanism for target-oriented evolution. We infer from this model that aa variability, structural changes and protein disorder may have different purposes in what we label as working, short-term and long-term evolved groups of tick salivary Kunitz peptides. For working evolved groups we propose that they possess high functional plasticity, but are highly susceptible to immune recognition. In short-term evolved groups they possess low functional plasticity and less immune recognition susceptibility. Long-term evolved groups have extremely low functional plasticity and possess low immune recognition susceptibility and may possess immune escape mechanisms other than antigenic variability.

Identifying single Kunitz-domain peptides from the I. ricinus transcriptome and peptide sequence alignments
We chose all full-length, secreted and non-truncated single Kunitz-domain peptide sequences from the SG I. ricinus transcriptome [18]. We then submitted all sequences to the Pfam server [72] to verify and identify them as single Kunitz-domain peptides. To further verify that these sequences are classical Kunitz-domain peptides we performed a protein BLAST alignment using BPTI (GenBank accession no.: P00974) as the query. As a final filtering point we submitted the remaining sequences to SignalP 4.0 server [73], since we were only interested in secreted peptides. All signal peptides were removed from all sequences and this resulted in a total of about 200 mature, single Kunitz-domain peptide sequences.
Only sequences that possessed a minimum of 6 Cys residues and a maximum of 7 Cys residues were used for our analysis. The majority of functionally described tick Kunitz peptides possess 6-7 Cys and usually >7 Cys are proteins with more than a single Kunitz-domainwith a few documented exceptions [26,40]. We specifically wanted to investigate the archetypal Kunitz-domain peptides of I. ricinus (hence, those possessing more or less 6 Cys residues). All peptides were divided into groups based on their Cys motif (i.e., the number intra-Cys residues) and a group was defined as whether there were two or more representative sequences for each motif. We chose a minimum of two sequences as a criterion since G3 has been previously reported by Dai et al. [23] and we only found two sequences for this group in the I. ricinus salivary transcriptome. This grouping resulted in a total of 145 sequences that formed eleven groups (G1-G11). We additionally found several Kunitz-domain peptides (n = 7) in the I. ricinus salivary transcriptome that possessed less than 6 Cys residues that were included in our phylogenetic analysis. These small peptides were defined as 'simple Kunitz' (SK) due their atypical Cys motif compared to archetypal Kunitz-domain peptides. Due to limitations of dealing with extremely small peptides and variations on their Cys motif, these SKs were used in only a few of our analyses. Thus, a total of 152 I. ricinus salivary Kunitz peptide sequences were used in our subsequent analyses.

Tertiary protein modeling and disulfide bridge predictions
To predict the disulfide bond patterns we first modeled a representative of each group (those indicated in Figure 1) using the evolutionary protein model building web-server Phyre2 [74]. Using the Schrödinger's Maestro software [75], each predicted 3D model was minimized, prepared, refined and re-minimizedminimization was used to remove steric clashes. Minimization uses a Truncated Newton algorithm [76], an OPLS-AA force field and a surface generalized Born (SGB) implicit solvent [77]. The total CPU time for minimizing was 2-3 min. After minimization, Cys resides were connected, hydrogen atoms were added, and the overall hydrogen-bond network was optimized using the Schrodinger's Protein Preparation Wizard [78], which optimizes the entire hydrogen bond network by means of side chain sampling. The Protein Preparation Wizard analyzes the structure then builds a cluster of hydrogen bonds and 100000 Monte Carlo orientations for each cluster. The algorithm determines the quality of the hydrogen-bond network to produce an optimized structure. The modeled structures were then reminimized. The disulfide bond patterns were verified using the web-server BRIDGED [79,80].
Amino acid (aa) variability (Shannon entropy) and secondary structure analysis Shannon entropy analysis [51] was used to calculate the aa variability. The Shannon entropy (H) was calculated for every position with the following equation: i¼1 P i log 2 P i P i is the amount of each aa type (i), and M is the number of aa types (maximum 20 aa). When H ≥2.0 the site is considered variable, while H ≤1.0 count for highly conserved positions [81]. Secondary structure was predicted using the position-specific scoring matrices method [82] from the PSIPRED server [82,83]. The presence (or not) of β-strands and α-helices, their length, amount, and position in the sequences were used as criteria for 2D model differentiation.

Nucleotide substitution rate analysis (SNAP and Datamonkey server)
We first performed a codon-based alignment using the GUIDANCE server [84] (MAFFT, Localpair, maxiterate 1000 with 100 bootstraps) to analyze the nucleotide substitution of the Kunitz encoding transcripts. We then submitted this codon-based alignment to the SNAP server that employs the Nei and Gojobori method [85]. Nei and Gojobori used the following equations to denote the type of substitution based on the assumption that a mutation in the second nucleotide position of any codon will cause a non-synonymous substitution, but mutations at the first or third position may cause either a non-synonymous or a synonymous substitution: Where s i and n i are the respective sites for synonymous (S) and non-synonymous (N) substitution for the ith codon with the total number of codons for the sequences studied (r). Nei and Gojobori [85] also used the following formulas to estimate the synonymous (d S ) and the non-synonymous (d N ) substitutions per site by applying the Jukes-Cantor formula [86]: We then used d S and d N to identify whether a group has undergone positive (>1) or purifying selection (<1); where p S and p N are the synonymous differences per synonymous site and non-synonymous differences per nonsynonymous site, respectively. Standard deviation (σ) and variance were also calculated by the SNAP server using the methods described by Nei and Gojobori [85].
Under the Datamonkey server [58] we used the algorithm PARRIS [60] to detect positive selection and the SLAC algorithm [59] was used to calculate the ratio d N /d S . A two-tailed extended binomial distribution set at p-value < 0.05 was used to assess significance of both algorithms. The SLAC and PARRIS algorithms use a neighbor-joining tree with a maximum likelihood for branch lengths and substitution rates. A few of our codon alignments contained multiple segments (i.e., recombination) as analyzed by GARD [87], therefore the base frequencies and substitution rates obtained from GARD were inferred together from the entire alignment, while branch lengths were fitted to each segment separately. Based on these analyses a global d N /d S ratio was obtained.
The SLAC algorithm was also used to detect which nucleotide substitution site were positively or negatively selected. For each synonymous and non-synonymous nucleotide substitution site, four measurements were made: normalized expected (ES and EN) and observed numbers (NS and NN). The SLAC algorithm then calculated: If dN < dS a codon was negatively selected and if dN > dS a codon was positively selected.

Phylogenetic analysis of Kunitz-domain peptides from hematophagous arthropods
All 152 I. ricinus Kunitz peptide sequences were used for a phylogenetic reconstruction of the Kunitz-domain family within ticks. Accordingly, we searched the NCBI database for Kunitz sequences from hard and soft ticks, and additionally spider sequences were included as outgroup. Similar to the I. ricinus Monolaris, only secreted and full-length Kunitz-domain proteins with up to 7 Cys residues in their aa sequence were chosen. All candidate proteins were re-confirmed as single Kunitz-domain peptides using the Pfam server [72]. Additionally, sequences of single Kunitz-domain proteins that have been functionally described were also included into the phylogenetic analysis.
A codon alignment was constructed using TranslatorX [88] and then all mature 249 single Kunitz-domain protein sequences were aligned using the program MAFFT v7 [89]. For the multiple sequence alignment we applied an iterative refinement method (L-INS-I) with the BLOSUM 62 matrix (New scoring scheme, gap opening penalty: 2.0, Offset value: 0.5). The protein alignment was then used in TranslatorX to align the corresponding Kunitz nucleotide sequences of ticks and spiders to create the final codon alignment. All unaligned, flanking sequence regions of the codon alignment were trimmed and all gaps were closed apart from gaps that were introduced by up to two sequences only (final sequence lengths: 126 bp, Additional file 5). The final alignment was analyzed using jmodeltest v2.1.2 in order to select the best-fit nucleotide substitution model for the phylogenetic reconstruction [90]. A phylogenetic tree was constructed under the maximum likelihood (ML) optimality criterion and the GTR model with a gamma distribution of substitution rate and estimated base frequencies using the program RAxML (HPC v7.2.6) [91]. Tree searching and bootstrapping were performed simultaneously (−f a, 1000 bootstrap replicates).

Divergence time estimation
The codon alignment of coding single Kunitz-domain proteins was also used to estimate divergence times among ticks using BEAST v1.7.5 [92]. For this Bayesian analysis, the following four taxon sets and divergence times (unit: million years ago, MYA) were calibrated according to Jeyaprakash and Hoy [63] using normal prior settings: Araneae/Scorpions/Pycnogonida/Acari 459 ± 18 MYA (included taxa: all tick sequences), Argasidae 214 ± 28 MYA (included taxa: all argasid sequences), Prostriata 196 ± 27 MYA (included taxa: all prostriate sequences) and Metastriata (included taxa: all metastriate sequences) 134 ± 22 MYA. All tick Kunitz-domain sequences in the Araneae/Scorpions/Pycnogonida/Acari taxon set were enforced to be kept monophyletic. Default settings were used for all other prior and operator settings in all analyses. Similar to the ML analysis, a GTR substitution model with a gamma rate distribution across sites and estimated base frequencies were applied. Two independent Beast runs with four independent Markov-chain Monte Carlo (MCMC) chains and 100,000,000 generations were performed. The starting trees were randomly generated and all trees were sampled every 10,000 generations. Different molecular clock (uncorrelated exponential and lognormal relaxed clock) [93] and speciation models [94] (Yule process and birth-death process) were tested. A uniform prior distribution for the calculations of the mean branch rates under the uncorrelated lognormal or exponential relaxed molecular clock was set. All BEAST analyses were assessed if each parameter converges on the same posterior distribution in the MCMC runs and if an effective sample size was reached using Tracer v1.5 (http:// tree.bio.ed.ac.uk/software/tracer). Additionally, the marginal likelihoods were estimated using path sampling and step stone sampling methods of each analysis in order to select the appropriate model for our analysis [95]. According to our evaluation of the different BEAST runs, the uncorrelated relaxed lognormal molecular clock model was chosen with the Yule tree prior for speciation. The two tree files from the latter BEAST analyses were combined using LogCombiner v1.7.5 from the BEAST software package using a burn-in of 10% of all sampled trees was set. Finally, TreeAnnotator v1.7.5 was used to summarize all trees in order to obtain a single maximum clade credibility tree including mean node heights.

Additional files
Additional file 1: Tertiary structural representation of G6 and protein structural alignment of each group representative. The tertiary structure of the archetypal Kunitz BPTI (A; PDB: 1BPI) and the modeled G6 representative (B; Ir2-1983) show the disulfide bridges (indicated by roman numerals), loops (L1 and L2), the beta-sheets (β1-β2) that form the β-hairpin, and the alpha-helices (α0 and/or α1). Both BPTI and G6 are colored from the N-terminus (blue) to the C-terminus (red). The tertiary structural alignment in Panel C depicts that the Cα protein backbone for each group representative do not drastically deviate from BPTI (color codes for each structure is presented at the far right).
Additional file 2: Tertiary structural deviations. The root mean square deviation (rmsd) for each group representative and BPTI were calculated using the Protein structural alignment tool, from the Maestro molecular modeling platform [75].
Additional file 3: Changes in putative secondary structure of I. ricinus Kunitz peptides. The PSIPRED server [82,83] was used to predict the secondary structure. The different Kunitz members were grouped into 22 secondary structural models. Kunitz members shown represent prototypes for the corresponding Kunitz secondary structures. The shaded grey box is the SKs (<6 Cys residues) and the graphical legend is at the bottom righthand corner. sequences and from other tick and spider sequences was constructed using TranslatorX [88]. Firstly, all mature 249 single Kunitz sequences were aligned using the program MAFFT v7 [89]. Secondly, the protein alignment was used in TranslatorX to align the corresponding Kunitz nucleotide sequences of ticks and spiders to create the final codon alignment. All unaligned, flanking sequence regions of the codon alignment were trimmed and all gaps were closed apart from gaps that were introduced by up to two sequences only (final sequence lengths: 126 bp).

Competing interests
The authors declare that they have no competing interests.
Authors' contributions AS drafted the manuscript and performed and interpreted the phylogenetic analysis. ACC drafted the manuscript and performed and interpreted aa variability and secondary structure analysis. JK assisted in drafting the manuscript and in interpreting the results. JJV developed the overall concept, assisted in drafting the manuscript, organized the I. ricinus Kunitz sequences and Illumina data and, performed and interpreted the substitution rate, structural analysis, and bioinformatics. All authors read and approved the final manuscript.
Authors' information Alexandra Schwarz and Alejandro Cabezas-Cruz joint first authorship. Jan Kopecký and James J Valdés joint senior authorship.