Differential gene expression analysis
The mechanisms underlying resistance and susceptibility against ticks of phenotypically divergent Braford hosts were addressed using RNA-Seq in cattle skin. With this approach, the extension of genetic modulation induced by ticks could be better examined in terms of differential gene expression and their collective participation in resistance (and susceptibility) at system biology levels. The gene expression profiles found before challenge (pre) could be partially explained by the immunity previously acquired by hosts, since they were not naive. Additionally, 24 h after the fourth tick infestation (post), we expected to see both innate and acquired immune responses, as already demonstrated for re-infested animals in skin-level studies of gene expression [14, 22, 23, 28].
Through the intergroup comparisons, we investigated the gene expression profiles of R and S hosts, looking for candidate genes that would explain phenotypic differences, according to the infestation challenge: Rpre vs. Spre (1) and Rpost vs. Spost (2). Both differential expression analyses revealed highly redundant DEG lists, suggesting that mechanisms of resistance in Braford cattle were present before challenge. Natural infestations before the first skin sampling, the time of skin sampling after infestation (24 h) or even the absence of the tick attachment site at the biopsy could contribute to these findings. Almost all annotated DEGs were classified as TopDEGs with biological functions correlated to cattle resistance against ticks. Overall, they were involved in hemostatic and immune mechanisms, many of them already associated with tick infestation, such as coagulation, iron metabolism, and inflammation [42,43,44]. It is conceivable that the modulation of these functions could help hosts achieve a better initial cellular response against tick antigens when tick feeding starts.
AWAT1, which had the highest variation in the intergroup comparisons, encodes an enzyme that acts in lipid metabolism and sebum composition in the skin [45]. Lipid metabolism has a role in inflammation control [46] and semiochemistry [47], and probably also in anti-tick resistance [24]. The tachykinin family (represented by TACR2 gene) was widely associated with inflammatory response, playing a role in wound healing and immune cell differentiation [48, 49]. Seric amyloid (SAA) family members were associated with anti-tick resistant response [50], and SAA3 expression was up-regulated in R hosts in addition to the intragroup comparison (Additional file 2: Table S2). Amine oxidase copper containing 3 (SAO) gene product is committed to leukocyte trafficking [51, 52], a function already associated with anti-tick responses [53], as well as in vascularization and tissue organization. Genes, such as bolA family member (BOLA) and tetratricopeptidase repeat domain 36 (TTC36), are related to immune mechanisms, including antigen processing and presentation. Moreover, genes involved in skin organization as collagen type XI alpha 1 (COL11A1), laminin alpha 1 (LAMA1) and agouti signaling protein (ASIP) may also contribute physically to tick bite healing, unsuccessful attachment, and establishment of the feeding site. Finally, the ectonucleotide pyrophosphatase/phosphodiesterase 3 (ENPP3) gene encodes an enzyme committed to the hydrolysis of extracellular nucleotides, also known as CD203c. This enzyme is better described as an activation marker of basophils, a cell type associated with acquired anti-tick resistance [25, 54, 55].
Genes modulated only before infestation (Rpre) are probably involved in primary cellular functions that may prepare R hosts to promptly react to infestation. DPYSL4 gene is described as a target of p53 in the apoptotic response to DNA damage [56]; once up-regulated, it could contribute to the apoptosis in the tick bite site. MTHFS is associated with folate turnover rate and depletion [57], and folate is essential for healthy skin [58]. Finally, over-expression of UTS2R, even in the absence of active infestation, suggests better control of vascular dynamics and osmoregulation in R hosts.
On the other hand, genes modulated only in Rpost may indicate putative mechanisms of protective response induced by the presence of ticks. Although the proteins encoded by genes such as UPK3A, TET1 and RBP1 do not play a direct role in immune cell functions, they indirectly help immune response. The function of UPK3A gene product has been linked to epithelial differentiation [59]. TET1, inhibited in Rpost, encodes a transcription repressor of a subset of genes, including interleukin 1 beta (IL1b), IL8, IL23, intercellular adhesion molecule 1 (ICAM1) and chemokine C-X-C motif ligand 1 (CXCL1) [60]. Most of these genes, modulated by TET1 expression, were up-regulated in Rpost (Additional file 2: Table S2). RBP1, also up-regulated in R hosts, encodes the carrier protein of retinol (vitamin A alcohol) transport, making it accessible in peripheral tissues. Besides the role in epithelial tissues, the hormone-like properties of vitamin A binding to nuclear hormone receptors retinoic acid receptors (RAR) and peroxisome proliferator-activated receptor (PPAR) affect immunity [61, 62]. This gene was up-regulated in Rpost compared to Spost and down-regulated in Spost compared to Spre, suggesting its involvement in tick resistance. Although the MYADML family is still incompletely annotated, some products have been described as modulators of cell spread and migration to the skin [63]. Three members of this family were modulated in Rpost with contrasting modulation (log2 FC -5.88, − 2.4 and 6.33; Table 1), which could indicate antagonistic effects regarding response to tick infestation.
In the intragroup comparisons, Rpost vs. Rpre (3) and Spost vs. Spre (4), we were looking for genes elicited by tick challenge, according to the resistance group studied. They revealed extensive gene lists, with thousands of DEGs broadly distributed over the genome, showing both redundancy and high correlation between hosts DEGs as well as non-redundant DEGs for R or S hosts (Additional file 2: Table S2). Critical functions against infestation, such as immune responses, coagulation, vascularization, and ion transport, were shared between hosts. This indicates that tick infestation caused a comprehensive response in Braford skin, regardless of being an R or an S host. On the other hand, non-redundant DEGs with potential anti-tick functions may represent mechanisms of response that contribute to differentiating hosts regarding resistance. This hypothesis was supported by the activation of critical functions in R hosts with massive participation of all immune cells (Table 2), which had the potential to protect them. Among S, cellular activation was not predicted for granulocytes (Table 3), an important cell type for tick resistance [19, 23, 54]. Overall, the data suggest the relevance of cytokines in R skin reactions. The inhibition of grooming predicted for R hosts is probably due to the time of skin sampling (24 h) when larvae had already been rejected [64]. Although very little is known about this, grooming behavior is considered to have an impact on tick load [15] and probably to ectoparasites other than ticks as well.
More than a hundred immune-related genes encoding cytokines, chemokines, CD markers, acute phase proteins, complement proteins, integrins, and TFs were found (Additional file 2: Table S2). Our results revealed novel DEGs associated with resistance/susceptibility, DEGs that support previous results and predictions, and DEGs that diverge with authors comparing gene expression between taurine and zebuine breeds [14, 21,22,23,24,25], within other cattle breeds [28,29,30,31] or murine models [37, 65]. According to Piper and colleagues [21, 23], CXCL2 and CCL2 were modulated only in S hosts; however, we found them also modulated in R. CCR1 and IL2R were described as modulated in R and CD14 in S [26], in contrast with our results where CCR1 and CD14 were modulated in both hosts, and IL2RG up-regulated only in S. IL8 was up-regulated in both R and S, but was previously described as down-regulated in R [33] and up-regulated in S [23]. IL13RA1 was described as up-regulated in S and down-regulated in R [28], whereas we found IL13RA1 to be up-regulated in S and IL13RA2 to be up-regulated in R. Regarding T cell participation in resistance [14, 22, 23, 25, 27, 30, 31], CD4 was up-regulated only in R hosts, in agreement with previous results. These inconsistencies may be the result of a myriad of factors, including sensibility of different techniques for gene expression analysis, experimental design, epistatic interactions with the genetic background and environmental effects, and indicate the need for extensive research on specific tick-host relationships.
Enrichment analysis (EA) based on functional ontologies
In the intergroup comparisons, similar protein classes were enriched among TopDEGs, although the number of TopDEGs was insufficient for network analysis. In the intragroup comparisons, S hosts showed more enriched protein classes than R. Functional enrichment was applied to TopDEGs (Additional file 3: Table S3) as this analysis is a powerful tool to determine the relevance of a gene in the context of a pathway. Some processes, which can be crucial to managing the interference of tick in the skin immune biology were over-represented in both hosts, including inflammation, chemotaxis, and immune response. Connective tissue degradation, matrix remodeling and immune process networks found only in R hosts likely indicate mechanisms they use to respond to ticks which are not elicited in S.
Even though both hosts had enriched leukocyte chemotaxis, differences between R and S in the modulation of this network may contribute to susceptibility (Additional file 3: Table S3A). For example, CCR3 was a TopDEG up-regulated only in R, although the corresponding ligand CCL13 was up-regulated in both phenotypes, thus suggesting differences in CCL13 effectiveness. CCL13 is involved in chemotaxis of eosinophils, basophils, and Th1 cells, all of which were previously associated with resistance [19, 23, 54]. The modulation of cytoskeletal proteins may also contribute to more efficient chemotaxis. On the other hand, the down-regulation of CCL20 in R hosts could indicate inhibition of the recruitment of lymphocyte subpopulation to the resistant skin [66]. The over-expression of MIP-1 alpha and I-TAC in S hosts are other examples: the former encodes a chemokine involved in inflammation, especially in attracting neutrophils and macrophages, which are known to have their functions impaired by tick saliva [43, 67, 68]; I-TAC, a chemoattractant of activated lymphocytes produced mainly by basal keratinocytes [69], did not find its receptor coding gene C-X-C motif chemokine receptor 3 (CXCR3) modulated, which was expected in activated T cells responding to infestation. Although tick infestation also elicits a cellular response in S hosts, it may attract less efficient cells to their skin.
Additionally, Th17-derived cytokines network was also over-represented among TopDEGs from both hosts (Additional file 3: Table S3A). Th17 cells belong to a distinct subset of CD4+ Th lymphocytes characterized as preferential producers of IL17A, IL17F, IL22, and, in humans, IL26. The receptors for IL17 and IL22 are broadly expressed in various epithelial tissues. The effector cytokines released from Th17 cells affect different cellular populations in inflammatory sites, mediating host defense through the activation of several signaling pathways. The participation of the IL17 family in cattle anti-tick response has already been proposed [25].
Finally, Jak-STAT pathway was also enriched for both hosts (Additional file 3: Table S3A). In this pathway, membrane-associated tyrosine kinases lead to the phosphorylation of signaling proteins and transcription factors such as STATs, regulating the expression of genes involved in the inflammatory response, cell survival, and cell cycle in response to numerous cytokines and growth factors. S hosts showed down-regulation of CNTF whose product, in the nervous system, appears to protect tissue during inflammation, and SOCS2, a Jak-STAT regulator coding gene whose product interacts with major molecules of signaling complexes to block further signal transduction.
Connective tissue degradation and ECM-remodeling networks were over-represented only among TopDEGs from R. These are two particularly important processes which require the over-expression and activation of proteinases able to destroy the connective tissue, allowing non-resident inflammatory cells arrival to the skin [23]. Most proteolytic enzymes are MMPs, and proteins including disintegrin and metallopeptidase domain (ADAMs), KLKs, elastases, trypsin, and others also participate. The MMPs comprise a family of enzymes that can collectively degrade all components of the ECM. They are the main proteases involved in ECM remodeling and are negatively regulated by TIMPs. The ability to modulate the skin immune-biology to impair tick success may represent a great difference between R and S hosts.
Inflammation via IL6-signaling was also over-represented only in TopDEGs from R hosts, in contrast with previous reports in which both R and S hosts show up-regulation of IL6 [25]; it is also reported in mice-infested skin, a model considered susceptible to tick [65]. This cytokine is secreted during inflammation by macrophages, endothelial cells, and fibroblasts. Its primary effect is to induce the secretion of acute-phase proteins, whose role as mediators of inflammatory responses has already been associated with resistance against ticks [24, 50]. IL6 also participates in the differentiation of some cell types, particularly in inducing B lymphocytes to differentiate into plasma cells; it can also lead to anti-apoptotic effects [70].
Given this, the up-regulation of IL6 together with many acute-phase protein coding genes and IL22 may create a favorable scenario for the maintenance of healthy skin barrier, since IL22 plays a critical role in skin immunity, inflammation, and repair [71]. Th22 participation in tick resistance was supported by the pathway analysis using TopDEGs from R hosts. Previous works have demonstrated the massive infiltration of CD4+ cells in the skin of infested animals as well as an increased number of these cells in the bloodstream and lymph nodes [14, 22, 26, 27, 30, 31]. IL22 is a member of the IL10 family of cytokines [72] that acts via a heterodimeric receptor complex (IL22 receptor) consisting of the subunits IL22RA1 and IL10RB (Fig. 2). IL22 binding to IL22 receptor activates the Jak-STAT pathway, inducing the phosphorylation of JAK1, Tyk2, and of the transcription factors STAT1, STAT3, and STAT5. IL22 also activates three major mitogen-activated protein (MAP) kinase pathways, c-Jun NH2-terminal kinase (JNK), p38 MAPK and extracellular signal-regulated kinase (ERK1/2), which are possibly required for activation of c-Jun/c-Fos to elicit the expression of pro-inflammatory cytokine genes [73, 74]. The primary sources of IL22 are immune cells, whereas its receptor is found in non-hematopoietic cells, especially at outer body barriers such as skin and the digestive and respiratory tracts [71, 75, 76]. Many genes in the Th22 pathway were DEGs in R hosts after infestation, including IL22 itself and other genes involved in response to this cytokine in tissue cells, as the tyrosine kinases and c-Fos (Fig. 2). Additional file 2: Table S2 illustrates substantial evidence for the differential activation of IL22-mediated immune response in R hosts, with the up-regulation of IL20, IL22, HP, JUN, FOSL1, IL9R, MMP13, CCR3, colony-stimulating factor 2 (CSF2) and IL6 in R but not in S, suggesting a putative role of this immune mechanism in resistance against ticks in Braford cattle. S hosts also showed some DEGs in this pathway, though neither IL22 was up-regulated, nor IL6, whose product acts in the Th22 cell differentiation [71]. Robbertse and colleagues (2018) [19] suggested the modulation of cytokines, such as IL6, IL13, IL5, chemokines, and of receptors as CCL2 and CCR1, acute phase proteins, cell adhesion molecules, CD4 and CD14 molecules, with the resulting attraction of T and B lymphocytes and granulocytes to the skin, to be associated with resistance. Our results partially agree with the predictions of these authors, since the modulation of some DEGs was redundant between R and S Braford hosts.
Network analysis
This analysis was carried out to identify the central mechanisms and to better explain contrasting host phenotypes based on their differential gene expression. It showed partially overlapping size-limited sub-networks, expanded from the TopDEGs lists, giving preference to objects with more connectivity considering that highly connected genes tend to be involved in similar biological functions.
WNT-signaling, the most relevant pathway predicted for R hosts (Additional file 4: Table S4 and Fig. 3a), is an evolutionarily conserved pathway that regulates the key process of organogenesis during embryonic development, cell fate determination, cell migration and polarity [77]. Recently, this pathway has also been associated with functions in immune response [78, 79] and skin biology as well [80,81,82]. Its effects on inflammation can be both pro- and anti-inflammatory, depending on the stimulus, cell type, activation context, and its crosstalk with other signaling pathways [78, 83]. It becomes inflammatory in infections with pathogenic bacteria and in inflammatory bowel diseases [78, 84]. However, it has also been associated with the induction of tolerogenic dendritic cells (DCs), which participate in immune response regulation mainly through the inhibition of Th1/Th17 cytokines [79]. It also participates in skin wound healing, through its ability to induce self-renewal and proliferation of skin stem cells [80, 82]. Finally, this pathway plays a critical role in many physiological processes of the skin from the earliest stages of development to postnatal control of hair cycling, inter-follicular epidermis maintenance as well as hair and skin pigmentation [81]. Furthermore, the deregulation of this pathway has already been associated with skin diseases [81] and host-tick interface [65]. Hair coat morphology and pigmentation were also associated with tick counts [85, 86]. Additionally, the WNT-signaling pathway was predicted in our analysis of genes mined from TagSNPs associated with resistance to ticks in the same breed [36] (Additional file 5: Table S5 and Fig. 3b). The importance of wound healing and structural proteins for tick resistance has been suggested [21, 37, 65]. Our results suggest the participation of WNT pathway in anti-tick resistance through an inflammatory scenario comprising the regulation of cytokines, receptors, matrix proteinases, and transcription factors and, at the same time, through its role in skin healing and remodeling, which hinders tick success.
On the other hand, the most relevant pathway for TopDEGs from S hosts includes MMP1, OSM, stromelysin-1, CCL2, and IL3 (Additional file 6: Table S6 and Fig. 4), which act on the regulation of immune response. This pathway included molecules related to tissue remodeling, such as OSM, MMPs, and TIMPs, potentially activating Jak-STAT and MAPK signaling cascades [87]. Most of the modulated genes in this pathway were DEGs (not TopDEGs), although Piper and colleagues [23] suggested ticks can more easily modulate the skin of S hosts.
Considering the complexity of host anti-tick resistance, which involves multiple mechanisms of activation and regulation of gene expression to drive the cellular and immune responses, we investigated the role of TFs as key molecules to orchestrate this modulation. Despite the general lack of knowledge concerning the regulation of gene expression in cattle, the TF lists were filtered according to the recently published manually-curated bovine TF compendium [41] to show only bovine TF candidates. For R hosts, Fos family members (c-Fos, c-Jun, Fra1, and FosB) were predicted as the most relevant TFs, with the latest two being TopDEGs (Additional file 7: Table S7), reinforcing the complexity of the regulation of gene expression in R host response against ticks. For S hosts, a nuclear factor kappa B (NF-kB) family member was the main TF based on its g-score (Additional file 8: Table S8), although Fos family members were also identified with lower relevance. Even though we extensively analyzed the data integrating different approaches, it was impossible to cover all aspects of the complexity of the host anti-tick response in this study. Experiments designed to compare different skin biopsies, using co-expression analysis and proteomic approaches could certainly drive the next steps of the investigation of this complex trait.