Skip to main content

Transcriptional changes in the peripheral blood leukocytes from Brangus cattle before and after tick challenge with Rhipicephalus australis

Abstract

Background

Disease emergence and production loss caused by cattle tick infestations have focused attention on genetic selection strategies to breed beef cattle with increased tick resistance. However, the mechanisms behind host responses to tick infestation have not been fully characterised. Hence, this study examined gene expression profiles of peripheral blood leukocytes from tick-naive Brangus steers (Bos taurus x Bos indicus) at 0, 3, and 12 weeks following artificial tick challenge experiments with Rhipicephalus australis larvae. The aim of the study was to investigate the effect of tick infestation on host leukocyte response to explore genes associated with the expression of high and low host resistance to ticks.

Results

Animals with high (HR, n = 5) and low (LR, n = 5) host resistance were identified after repeated tick challenge. A total of 3644 unique differentially expressed genes (FDR < 0.05) were identified in the comparison of tick-exposed (both HR and LR) and tick-naive steers for the 3-week and 12-week infestation period. Enrichment analyses showed genes were involved in leukocyte chemotaxis, coagulation, and inflammatory response. The IL-17 signalling, and cytokine-cytokine interactions pathways appeared to be relevant in protection and immunopathology to tick challenge. Comparison of HR and LR phenotypes at timepoints of weeks 0, 3, and 12 showed there were 69, 8, and 4 differentially expressed genes, respectively. Most of these genes were related to immune, tissue remodelling, and angiogenesis functions, suggesting this is relevant in the development of resistance or susceptibility to tick challenge.

Conclusions

This study showed the effect of tick infestation on Brangus cattle with variable phenotypes of host resistance to R. australis ticks. Steers responded to infestation by expressing leukocyte genes related to chemotaxis, cytokine secretion, and inflammatory response. The altered expression of genes from the bovine MHC complex in highly resistant animals at pre- and post- infestation stages also supports the relevance of this genomic region for disease resilience. Overall, this study offers a resource of leukocyte gene expression data on matched tick-naive and tick-infested steers relevant for the improvement of tick resistance in composite cattle.

Peer Review reports

Background

Rhipicephalus microplus is a tick species complex of hematophagous arthropods that parasitise cattle in tropical and subtropical countries, of which the representative species found in Australia is Rhipicephalus australis [1, 2]. Heavy tick infestation and the associated risk of transmission of tick-borne pathogens have detrimental consequences for animal health, welfare, and production in the cattle industry. Global economic losses due to cattle ticks are estimated to be US$22–30 billion annually [3,4,5]. Strategies to control the spread of cattle ticks include treatment with acaricides, pasture spelling, the use of tick-resistant cattle, or a combination of all of these [6]. However, treating animals with synthetic acaricides raises concerns about potential residual effects in meat and milk for human consumption, as well as the challenge of the development of chemical resistance [7, 8]. Genetic selection to improve host resistance to ectoparasites remains the most sustainable approach to minimise disease burden and enhance animal welfare in the cattle industry [5, 9, 10].

Female cattle ticks have an average parasitic phase of ~ 21 days on the host. During this period, the ticks secrete a diverse range of bioactive molecules to modulate host responses in order to achieve feeding success and survival [11]. This complex range of secreted protein and non-protein tick salivary molecules has been reviewed elsewhere [12]. Some of the host processes that are known to be disrupted by the hematophagy of cattle ticks include blood coagulation, cytokine secretion, cell adhesion, as well as leukocyte recruitment to the feeding site. Conversely, host resistance in cattle manifests by preventing attachment of tick larvae, which results in early larval death, as well as limiting hematophagy, thus reducing the reproductive success of the female adult ticks [11].

Host resistance is a moderately heritable trait which can range from high (tick-resistant) to low (tick-susceptible) as determined by the number of engorging adult ticks (4.5-8 mm) following a natural or artificial larval infestation [13, 14]. Bos indicus cattle generally have very high tick resistance, whereas Bos taurus cattle are mostly susceptible [14]. However, stable composite breeds (B. indicus x B. taurus) range widely from low to high host resistance as reported in Santa-Gertrudis [15], Bradford [16] and Brangus [17]. Targeting tick resistance in breeding programs is desirable, but measuring the phenotype requires collection of tick counts or scores on an individual animal basis, which is costly, laborious, and requires standardisation of exposure [5, 9]. Hence, knowledge of the genes controlling this trait will contribute towards the strategic development of biomarker-assisted selection methods for tick resistant cattle.

Gene expression and quantitative trait loci (QTL) studies have previously shown that host resistance to the cattle tick is controlled by multiple genes [16, 18,19,20,21,22,23,24,25]. Additionally, the complex interplay of innate, adaptive and mal-adaptive responses elicited during tick infestation have been shown to contribute to the expression of divergent phenotypes [15, 26,27,28,29,30,31,32,33,34,35]. As a result, identification of genes that can be used as universal biomarkers of tick resistance is particularly challenging. Transcriptomic profiling of resistant versus susceptible cattle has been well characterised at the skin level [16, 18, 19, 27, 36, 37], the primary site of tick attachment and feeding, whereas immunologically relevant tissues are still far less explored [21, 26, 31, 38]. Since blood collection is a less invasive procedure than skin biopsy collection for biomarker research, this study hypothesized that it is possible to use peripheral blood leukocytes to elucidate differentially expressed genes that contribute to divergence in host resistance to cattle ticks. Therefore, the aim of this study was to apply RNA-sequencing in the leukocytes of Brangus cattle (B. indicus x B taurus) collected before and after repeated artificial challenge with R. australis tick larvae to identify genes associated with 1) host response to tick infestation; and 2) divergent phenotypes of host resistance.

Results

Host resistance phenotyping

Thirty tick-naive Brangus steers were artificially infested with R. australis larvae for 12 consecutive weeks and host resistance phenotype was determined by weekly measurement of tick burden with scores (Fig. 1A). The scores were less variable between week 8 and week 15, and this period was considered to represent stabilisation of the phenotype, hence, the mean tick score (MTS) was used to rank animals. A summary of the mean tick scores for 10 animals selected as divergent in host resistance phenotype is given in Fig. 1B. The highly resistant steers (HR, n = 5) and least resistant steers (LR, n = 5) had a group MTS of 1.35 ± 0.49 and 4.08 ± 1.08, respectively, and the difference was significant (p < 0.001) (Fig. 1C). The scoring data table for all 30 Brangus steers is provided in Additional File 1.

Fig. 1
figure 1

Tick scoring for host resistance phenotyping of Brangus steers. A Heatmap representation of the tick scoring data collected from 30 Brangus steers across twelve timepoints showing clustering of the most resistant (bottom) and least resistant (top). Tick scores represent the estimated number of adult ticks that develop 21 days after a single infestation with R.australis larvae, where 1 = 0–50 ticks, 2 = 50–100 ticks, 3 = 100–200 ticks, 4 = 200–300 ticks, 5 = more than 300 ticks, blank = not scored. The vertical colour bar indicates the value that each colour represents in the heatmap from low (green) to high (red). B The mean tick score (MTS) ± standard deviation (SD) calculated across week 8–15 for 10 animals selected as divergent in host resistant phenotype. C Boxplots showing the comparison of MTS (week 8–15) values between the high (HR, n = 5) and low (LR, n = 5) host resistance phenotypes

Genomic estimates of Bos indicus content and relationship with mean tick score

Genomic estimates of B. indicus content (BIC) were used to investigate if there was an association between low tick burdens and high indicine content in this Brangus population. The BIC values obtained for 29 steers ranged from 25 to 49% with a median value of 40% (Additional File 2). The inverse correlation between MTS and BIC values was not significant (r = − 0.17, p = 0.37), as shown in the figure provided in Additional File 2.

RNA-sequencing statistics summary

Leukocyte samples from the HR and LR group at three timepoints (T0, T3, and T12) of the infestation trial were included in this RNA-Seq study (see Section “Methods” for further information). Following the removal of two RNA samples from the HR-T0 group, there were 28 sequenced libraries which generated, on average, 36,170,898 raw reads (100 bp single-end) per sample. Of these, 0.3% of the reads were discarded after adapter trimming and read quality control and 93.4% of the remaining sequences were uniquely mapped to the reference B. taurus genome ARS-UCD1.2. The average library size was 23,460,531 reads assigned as counts to 28,786 annotated genes in the reference genome. The RNA-seq mapping statistics per sample are shown in Additional File 3.

Differentially expressed genes associated with host response to tick infestation

To study the effect of tick infestation on leukocyte gene expression, comparisons were performed between timepoint sample sets T3 (n = 10) vs. T0 (n = 8), and T12 (n = 10) vs. T0.

The total number of expressed genes tested by edgeR was 13,984 in T3-vs-T0, and 13,865 in T12-vs-T0. Exploratory analysis with multidimensional scaling (MDS) plot showed sample clustering in the first dimension according to timepoint as expected, whereas the secondary dimension separated samples according to host resistance phenotype. This was more clearly evident in T0 than in the T3 and T12 sample sets (see Additional File 4).

Differential expression analysis in T3-vs-T0 identified 3065 significant DEGs (FDR < 0.05), of which 1674 were upregulated and 1391 were downregulated genes (Fig. 2A; Additional File 5). In the ranked gene list according to expression fold change, the top upregulated (logFC> 3) DEGs were PROCR, LOC104972252, ESYT3, PDLIM1, SPMD3, ALOX15, whereas the top downregulated (logFC< − 3) DEGs were MARCO, KIAA1324L, FOS, LOC515150, GRO1, MFSD4A (Table 1).

Fig. 2
figure 2

Differentially expressed genes (DEGs) in the leukocytes of tick-infested compared to tick-naive Brangus steers. A Volcano plot of DEGs from T3-vs-T0 comparison (3-week tick-exposed vs. tick-naive). B Volcano plot of DEGs from T12-vs-T0 comparison (12-week tick-exposed vs. tick-naive). The x-axis represents fold change and y-axis represents statistical significance. Symbols are shown for genes with fold change |logFC| > 2.5 at a significance threshold FDR < 0.05. C Venn diagram of upregulated DEGs between comparisons. D Venn diagram of downregulated between comparisons

Table 1 Top significant DEGs from the comparison between 3-week tick-exposed and tick-naive steers

The T12-vs-T0 analysis identified 1959 significant DEGs (FDR < 0.05), of which 989 were upregulated and 970 were downregulated genes (Fig. 2B; Additional File 6). In the ranked gene list according to expression fold change, the most upregulated DEGs (logFC> 2.5) were SPMD3, ALOX15, GATA, CCR3, LOC100297044, GPAT2, LOC100848006, HRH4, LOC104969122, ALOX5, PDLIM1, HBB, whereas the most downregulated DEGs (logFC<− 2.5) were FOS, GRO1 (Table 2).

Table 2 Top significant DEGs from the comparison between 12-week tick-exposed and tick-naive steers

There were a total of 1380 DEGs captured in both comparisons, which represents 45% of all DEGs from T3-vs-T0 and 70.4% from T12-vs-T0. Among these, 714 genes were upregulated and 666 were downregulated (Fig. 2C-D). The most differentially expressed genes (|logFC| > 3) including SMPD3, ALOX15, GATA1, CCR3, GPAT2, GRO1 had higher expression in T12-vs-T0 than in T3-vs-T0, whereas PDILIM1, FOS, LOC107133088, LOC514978, FCRL6 had higher expression in T3-vs-T0 than in T12-vs-T0. The vitamin K epoxide reductase complex subunit like 1 (VKORC1L1) was the most significant DEG in both comparisons (FDR < 5.4E-06). Interestingly, a relevant positional candidate gene for tick resistance RIPK2 [24] was found downregulated (logFC T12vsT0 = − 0.27, FDRT12vsT0 = 0.0026) in T12-vs-T0. Related to this, the interacting serine/threonine kinase 1 (RIPK1; logFC = − 0.3, FDR < 0.02) and gasdermin D (GSDM; logFC = − 0.4, FDR < 0.02) genes [39] were also found to be significant and downregulated in both timepoint comparisons, and gasdermin E (GSDE, logFC T3vsT0 = − 1.2, FDRT3vsT0 = 0.0024) downregulated only in T3-vs-T0 comparison.

Over-represented terms and biological pathways in response to tick infestation

Over-representation analysis of Gene Ontology (GO) terms and Kyoto Enrichment of Gene and Genomes (KEGG) pathways was performed to identify the specific biological functions associated to the DEGs detected in the comparisons between tick-exposed (3-week and 12-week) and tick-naïve steers. The GO biological process (BP) functions that were significantly enriched (p-adj < 0.05) in response to tick infestation were related to cell chemotaxis, blood coagulation, angiogenesis, inflammatory response, cell adhesion and cytokine secretion (Fig. 3A). Additionally, the KEGG enrichment analyses (p-adj < 0.05) showed that pathways including osteoclast differentiation (bta04380), cytokine-cytokine interaction (bta04060), and the IL-17 signalling pathway (bta04657) were significantly changed in host leukocytes upon tick infestation (Fig. 3B). For instance, the complement and coagulation cascade (bta04610), and arachidonic acid metabolism (bta00590) pathways were more significantly enriched in 3-week than in the 12-week exposure period. Conversely, the pathways Staphylococcus aureus infection (bta05150) and cholesterol metabolism (bta04979) were more enriched in the 12-week period than in the 3-week exposure period.

Fig. 3
figure 3

Functional enrichment analysis of the differential expression analysis between tick-exposed and tick-naïve steers. A Enriched GO Biological Process terms and B Enriched KEGG pathways in the differentially expressed genes (DEGs) from the comparison T3-vs-T0 (3-week tick-exposed vs. tick-naive) and T12-vs-T0 (12-week tick-exposed vs. tick-naive) performed by clusterProfiler. The dot colour represents significance (p-adjusted < 0.05) of the term, and dot size (GeneRatio) represents the ratio of input genes that are annotated in a term

The DEGs annotated in the top 5 over-represented KEGG pathways are shown in Fig. 4. The highly downregulated genes FOS and GRO1, FOSB and JUN genes were common to various pathways including the IL-17 signalling pathway, osteoclast differentiation and cytokine-cytokine receptor interaction pathway. Highly upregulated genes such as PROCR and ALOX15 were annotated in the complement and coagulation cascades and arachidonic acid metabolism pathway, respectively. A cluster of upregulated genes in the cytokine-cytokine receptor interaction pathway included CXCL10, PPBP, PF4, CCR3, IL5RA in short-term tick exposure, and CCR3, CCR4, IL9R, BMPR1B and LOC100297044 in long-term tick exposure. The category network plots showing the fold changes of the featured genes in the enriched pathways are provided in Additional File 7.

Fig. 4
figure 4

Top enriched biological pathways in tick-exposed steers. Network plot of the enriched KEGG pathways (p-adjusted < 0.05) in the differentially expressed genes (DEGs) from the comparison T3-vs-T0 (3-week tick-exposed vs. tick-naive) and T12-vs-T0 (12-week tick-exposed vs. tick-naive). Light blue = enriched in T3-vs-T0; dark blue = enriched in T12-vs-T0)

Visualization of pathway enrichment results with the Pathview R package demonstrated that several components of IL-17 signalling pathway were primarily downregulated by tick infestation, including IL17RC, IL17RA, IKBKE, TRAF6, MAPK1, FOS resulting in the repression of chemokine genes GRO1, cytokine genes PTGS2, antimicrobial genes DEFB10, S100A7, S100A8, S100A9, except for a few genes including CASP8, TBK1 and CXCL10 which were activated at either or both stages of the infestation (Fig. 5). Additional KEGG maps featuring DEGs for the pathways represented in Fig. 4 are provided in Additional File 8.

Fig. 5
figure 5

Gene expression changes in the IL-17 signalling pathway. Representation of the differentially expressed genes that were enriched in the IL-17 signalling pathway (KEGG: bta04657). Red and green boxes indicate the up- and downregulated genes in response to tick infestation, respectively. Left-side coloured box represents a gene with differential expression in the T3-vs-T0 comparison (3-week tick-exposed vs. tick-naive). Right-side coloured box represents a gene with differential expression in the T12-vs-T0 comparison (12-week tick-exposed vs. tick-naive). Fully coloured box indicates a gene that was changed in both comparisons. Blank box indicates no changes in expression. Pathway data was sourced from the KEGG database [40,41,42] and rendered with the Pathview R package

Differentially expressed genes associated with divergent phenotypes of host resistance to tick infestation

Leukocyte gene expression changes between low (LR) and high (HR) host resistance steers were investigated at timepoints T0 (tick-naïve), T3 (3 weeks post-initial infestation), and T12 (12 weeks post-initial infestation). The three comparisons performed were: LRT0 (n = 3) vs. HRT0 (n = 5), LRT3 (n = 5) vs. HRT3 (n = 5), and LRT12 (n = 5) vs. HRT12 (n = 5).

The total number of expressed genes tested by edgeR was 14,326 in T0, 13,757 genes in T3, and 14,060 genes in T12. Based on a significance threshold of FDR < 0.05, there were 69 DEGs (27 upregulated, 43 downregulated) in T0, 8 DEGs (4 upregulated, 4 downregulated) in T3, and 4 DEGs (1 upregulated, 3 downregulated) in T12 (Fig. 6; Table 3). A total of 75 unique DEGs were identified in the comparison between LR and HR steers across the three timepoints, with only 3 common genes between T0 and T3 (LOC524810, FBN1 and C26H6orf52). Furthermore, the gene expression fold changes were higher at pre-infestation compared to the other two post-infestation timepoints. Genes of the bovine leukocyte antigen (BoLA) complex, such as BOLA-DQA2, BOLA-DQA5 were downregulated in LR steers prior to tick exposure, whereas BLA-DQB and BOLA-DQRB2 appeared upregulated in LR compared to HR steers after tick exposure at T3 and T12, respectively (Table 3; Additional File 9). A gene of interest, kallikrein-1 (KLK1) was found to be upregulated in LR steers (logFC =1.3; FDR =5.3E-03) before tick infestation.

Fig. 6
figure 6

Differentially expressed genes (DEGs) in the leukocytes of low (LR) compared to high (HR) host resistance Brangus steers. A Volcano plot of DEGs at T0 (tick-naive) timepoint. B Volcano plot of DEGs at T3 (3 weeks post-initial infestation) timepoint. C Volcano plot of DEGs at T12 (12 weeks post-initial infestation) timepoint. The x-axis represents fold change and y-axis represents statistical significance. Symbols are shown DEGs with fold change |logFC| > 1 and FDR < 0.05 in T0, and for the DEGs in T3 and T12 (FDR < 0.05)

Table 3 Significant DEGs from the comparison between low and high host resistance steers at different timepoints

The functional enrichment analysis did not report any significant GO terms or KEGG pathways in the lists of LR-vs-HR DEGs.

Overlap between genes associated with host response and host resistance to infestation

Using sample sets from high and low resistant phenotypes at three timepoints of the infestation trial, this study performed overall 5 types of gene expression comparisons. As described above, these comparisons were implemented to produce lists of DEGs between tick-exposed and naïve steers (T3-vs-T0 and T12-vs-T0) and between host resistant phenotypes (LRT0-vs-HRT0, LRT3-vs-HRT3, LRT12-vs-HRT12). Differentially expressed genes that could be relevant according to the effect of tick infestation and phenotype of interest were investigated among these lists and presented in Table 4.

Table 4 Differentially expressed genes by tick infestation and host resistant phenotype

ITGB3 and CXCR5 were leukocyte-expressed genes initially downregulated in tick-naïve steers with low host resistance, which then became upregulated in all steers following 3 and 12 weeks post-initial infestation. BOLA-DRB2 was the only gene that was upregulated by tick infestation, but that remained downregulated in low host resistance steers after 12 weeks. Other common genes such as C18H19orf18 and LTBP1 showed increased downregulation after 3 weeks of tick exposure, whereas the overall trend for the rest of the genes was towards upregulation in response to tick exposure.

Discussion

Selection for host resistance to ectoparasites is a desirable approach in the beef cattle industry as part of the overall genetic improvement strategy for animal productivity and health [4, 5]. Host resistance reduces the tick burden and therefore is a measurable phenotype, but in general, tick counts or scores are difficult to obtain, which constrains the ability to obtain large numbers of animal phenotypes for genetic improvement [9]. Although there is a heritable component for differences in host resistance, the dynamics of the host immune response also have an indispensable contribution, but these remain to be well defined [43]. The advances in high-throughput transcriptomics and bovine genomic resources offer the opportunity to put forward candidate genes for genomic selection programs.

Host resistance in Brangus cattle

This study showed that a group of naïve Brangus steers infested repeatedly with the same dose of ticks larvae exhibited divergent phenotypes of high (HR) and low (LR) host resistance after the sixth infestation (resistance at week 5, scored on week 8) which remained steady until the end of the trial. This result is broadly consistent with the infestation study in naïve Santa-Gertrudis cattle [15] which found that significant tick burden differences between resistant and susceptible hosts appeared at week 5 of their trial, although this study used tick counts instead of scores [15]. Importantly, there is no consensus regarding the timing for host resistance onset in the literature, as it will be mostly influenced by the underlying variation in immune response among individuals. Moreover, in natural infestation conditions with fewer timepoint observations spread over a long period of time, it may be more difficult to determine when an animal has acquired resistance. As such, it is highly relevant for the cattle industry to find a more objective and unbiased measure of tick resistant phenotype.

Tick resistance is moderately heritable and the proportion of B. indicus genetics in a breed can influence the variation in host resistance to ticks [5, 43]. A study comparing the resistance levels of Brazilian Brangus against purebred B. indicus (Nellore) cattle showed that under the same natural infestation conditions Brangus animals had a higher tick burden than the Nellore [44], thus attributing the resistance of the latter to the high level of indicine genetics. The Brangus breed is a composite with typically 3/8 of indicine (37.5%) and 5/8 of taurine (62.5%) composition, but a recent study of the genomic architecture of the American Brangus showed taurine composition values deviating from the expected towards 70.4 ± 0.6% [45]. The present study found that the genomic estimates of B. indicus content in this Australian Brangus steer population were, in general, consistent with the theoretical expectation showing a median of 40%, but there were also some individual extremes where the lowest had 25% and the highest 49% indicine composition, the latter being one of the animals classified in the highly resistant group. However, the results also showed the inverse correlation between the indicine content and the mean tick scores was not significant, which could be due to small sample size (n = 29, out of 30). Hence, based on this data, it cannot be concluded with certainty that no association exists between the variables tested, thus warranting further investigation with a larger animal population. At the same time, one may question whether the demonstration of an inverse correlation between average tick burden and indicine content is sufficient to explain host resistance, as it is still of great interest to determine which indicine alleles have been selected for host resistance in Brangus [46], and also whether there are functional candidate genes associated with the immunological response of the host to tick infestation.

To the best of our knowledge, this is the first study to collect timepoint transcriptome data from peripheral blood leukocytes of Brangus cattle to explore gene expression changes related to a) the host response to repeated artificial R. australis challenge, and b) the host resistance differences at pre- and post-infestation. The results showed that naïve steers (both HR and LR) with exposure to a short infestation period (3 weeks) produced 1.5 times more differentially expressed genes than exposure to a long infestation period (week 12). However, among the group differences within each timepoint, it was found that the number of differentially expressed genes was the highest at pre-infestation, and then decreased significantly at early and late post-infestation. Immune responses are highly complex and dynamic requiring the interplay of innate and acquired immunity, and therefore a stronger immune response is not necessarily the most protective given that it can result in immunopathology [47, 48]. It is also context-dependent, whereby the host needs to balance off the cost/benefit of eliciting immunity and removing or maintaining pathogen challenge [49]. With this into consideration, the significance of the gene expression results is discussed next.

Effect of tick infestation on peripheral blood leukocytes’ expression of immune and inflammatory pathways

The results showed a high number of differentially expressed genes across all animals (both HR and LR) in response to the short-term and long-term tick exposure which appeared enriched in GO terms including leukocyte chemotaxis, haemostasis, coagulation, and inflammatory response. These varied biological processes are consistent with a response to tissue damage inflicted by heavy infestation and host pathology in the attempt to eliminate the burden and have been also reported in gene expression studies in blood [31, 38] and skin [16, 18]. Moreover, tick infestation altered the expression of genes that participate in the immune system (IL-17 signalling, complement and coagulation cascades), cell signalling (cytokine-cytokine interaction), development and regeneration (osteoclast differentiation). Interestingly, these pathways appeared mostly downregulated in both short and long infestation periods. Whether these expression changes are mediated by tick salivary proteins needs further research. However, another possible explanation is that pathway downregulation in the leukocytes may have resulted from a reduction of the cell populations expressing these genes when moving out of the peripheral circulation into secondary lymphoid organs and skin [30, 31].

IL-17-mediated-immunity has emerged as an important host defence mechanism against pathogens and ectoparasites [16, 18]. In this study, several components of the downstream IL-17 signalling were downregulated by tick infestation including IL17RC, IL17RA, TRAF6, IKBKE, MAPK1, IKBKB, IKBKE, NKFBIA, FOS, the chemokine gene GRO1, cytokine gene PTGS2, and antimicrobial genes DEFB10, S100A7, S100A8, S100A9. Similar findings in mice studies have been reported for the downregulation of IL-17 receptors (IL-17R) during pathogen infection [50], and the host Th17 immunity during infestation with Ixodes scapularis ticks [51]. Sun et al. [52] found TRAF6 (tumor necrosis factor receptor-associated 6) to be a target of the tick salivary protein Ds Cystatin, effectively resulting in downregulation of TRAF6 in the Toll-like receptor signalling pathway. Gene expression changes in TRAF6 in cattle skin have also been previously reported in various cattle breeds following artificial tick challenge [19, 21].

Proinflammatory cytokines and chemokine contribute to host defence and inflammation by recruiting neutrophils and T-cells to sites of pathogen invasion [53]. The migratory capacity of activated T-cells depends on their antigenic experience and type of polarization (Th1 or Th2) which stimulates the selective upregulation of their chemokine receptors [54]. High expression of the chemokine receptor CXCR3 and ligand CXCL10 have been implicated with Th1 cell-mediated inflammation and chronic inflammatory conditions [55], and this study found the pair CXCR3/CXCL10 to be upregulated in response to long-term infestation. CXCL10 expression has also been reported in previous blood gene expression studies in tick-infested cattle although with discrepant patterns among breeds and resistance status [31, 38, 56], suggesting a relevant role for this chemokine gene in the adaptive immune response of cattle. Other studies have also reported expression changes in genes from the cytokine-cytokine receptor pathway in skin [16] and lymph node [26] transcriptomes, as well serum proteomes [57] of cattle exposed to ticks. Although establishing similarities in the cytokine/chemokine and receptor profiles across tick infestation studies is challenging due to differences in experimental design and timepoints of observation, it is clear that cytokine responses could contribute both to protection and immunopathology to tick feeding.

The immunological relevance of RIPK2 gene (receptor-interacting serine-threonine kinase 2) was first proposed by Porto Neto et al. [24]. The authors suggested this to be a positional candidate gene for tick resistance, demonstrating in a knockout model that RIPK2-deficient mice had reduced antibody production against tick salivary gland extract. Furthermore, Wang et al. [27] showed that RIPK2 expression was downregulated in the skin of Hereford-Shorthorn cattle of high and low tick resistance when infested with R. microplus. Similarly, the present study found that RIPK2, but also RIPK1, were downregulated in the leukocytes of Brangus steers following long-term tick exposure, thus this further supports the relevance of this gene, also at the transcriptional level, in host responses to tick infestation. More insights about the mechanism of host defence action involving the RIPK1 gene in a murine model of pathogen infection have been recently reported in the study by Chen et al. [39]. Here it is understood that RIPK1 activates gasdermin D (GSDMD) in macrophages and the related gasdermin E (GSDME) in neutrophils to promote host defence. In Brangus steers, both GSDMD and GSDME genes were also found to be downregulated by tick challenge, which could represent a reduction in the host innate immune signalling and therefore explain the overall enhanced suppression of host defence which benefits tick survival. Whether there is a role of tick salivary proteins in the downregulation of RIPK1 and gasdermin genes should be verified by future research.

Before tick infestation HR and LR steers reveal differences in expression of immune, tissue remodelling, and angiogenesis genes

This study found that leukocyte gene expression differences between HR and LR steers were the highest in the absence of tick challenge, which supports the hypothesis of genetic variations influencing the outcome of host resistance phenotypes among animals of the same breed. The surprising finding that gene expression profiles remain largely unchanged in the presence of tick burdens is in stark contrast to the previous findings of altered host defence responses to the cumulative effect of tick infestation. Several different factors may have contributed to the reduction of the leukocyte expression signals between phenotypes at post-infestation: rapid immune response signal (< 24 h), immunopathology interference with resistance, removal of circulating inflammatory cells by apoptosis or transmigration to the skin, immunosuppressive effects by tick feeding. The lack of temporal expression studies looking at several timepoints in the progression of bovine hosts from naïve to resistant states, as well as expression profiling of isolated immune cells, currently limits the ability to assess all of the factors mentioned previously. However, a similarity was found in a gene expression study in the lymph node of Bonsamara cattle, where there was a significant reduction in the number of DEGs in tissue infested with R. microplus adult ticks compared to the larval stage, concluding that lymphocyte maturation signals may be suppressed [26]. The number of animals in the group comparisons of both studies was small and therefore any implications of the commonalities found here require further investigation, where possible, in larger animal cohorts with consideration to the analysis of multiple host tissues and additional timepoints.

A key finding of this study was that observation of genes of the major histocompatibility complex (MHC), also known as bovine leukocyte antigen region (BoLA) had altered expression between HR and LR animals at pre- and post-infestation. The MHC of cattle is located in the bovine chromosome 23 and is highly polygenic and polymorphic, which allows for a wide expression repertoire of cell-surface glycoproteins that are essential for antigen processing and presentation to T cells, thus linking the innate and adaptive immune responses [58]. Several studies have reported the association of allelic diversity in BoLA genes and tick resistance/susceptibility, particularly in the class II DRB3 locus (BOLA-DRB3) [59,60,61,62]. In this study, the genes BOLA-DQA2 and BOLA-DQA5 presented as highly upregulated in LR compared to HR steers at pre-infestation, whereas after 3 weeks of repeated infestation, BOLA-DQB was found to be upregulated in HR animals. Additionally, BOLA-DRB2 is one of the genes that was found to increase in expression following repeated tick exposure in both phenotypes, but after 12 weeks of tick exposure, its expression remained higher in HR than in LR steers. In this study, there was no specific association of BOLA-DRB3 expression with resistance phenotypes, but it was observed to be upregulated in both phenotypes in response to short-term tick exposure. The contribution of these genes to tick resistance in Brangus cattle can be further explained by the findings of Goszczynski et al [46], which revealed signatures of positive selection towards indicine (Brahman) variants on the BoLA region of a Brangus population, including the location of DQA2, DQA5, DQB, DRB2, DRB3 genes. Therefore these genes should be considered as relevant candidates in host resistance to ticks.

In regards to the most relevant DEGs detected at pre-infestation, HR steers displayed increased expression of the chemokine receptor CXCR5 (C-X-C motif chemokine receptor 5) and ligand CXCL12 (C-X-C motif chemokine ligand 12), tissue remodelling genes such as FBN1 (fibrillin 1), COL26A1 (collagen type XXVI alpha 1 chain), and cell adhesion genes including ITGB3 (integrin subunit beta 3) and LAMA4 (lama subunit 4). CXCR5, which also was upregulated in response to tick infestation, is typically expressed by mature B cells and participates in the recruitment of naïve B cells into the lymph node, whereas CXCL12 promotes leukocyte chemotaxis [63]. FBN1 encodes an extracellular matrix glycoprotein and mutations in this gene are associated with deficient protein synthesis and weakening of connective tissue by perturbing microfibril assembly and function [64]. COL26A1 encodes a collagen type XXVI, although it has not yet been fully characterised in cattle, mice studies suggest it differentiates from other fibrillar collagen family subgroups and functions as an extracellular matrix component [65]. ITGB3, also upregulated following long-term tick exposure and enriched in the osteoclast differentiation pathway, encodes the beta3 subunit of the integrin heterodimer which is a cell adhesion receptor that mediates the attachment of cells to the extracellular matrix and signalling in osteoclast adherence to bone [66]. Integrin can stimulate changes in intracellular calcium levels, and an increase in the expression of calcium genes has been previously reported in skin of tick-infested cattle and linked to mechanisms of host resistance [36, 67]. LAMA4 expresses alpha4 subunit of laminins which are extracellular glycoproteins with an important role in the regulation of endothelial cell adhesion and this subunit is a major structural component of the basement membrane of developing blood vessels via complex interaction with integrins [68, 69]. In this way, transcriptional changes around cytokine signalling, tissue homeostasis and angiogenesis could explain how animals may be able to better resist tick infestation and these results correlate well with previous findings in skin transcriptome studies [27].

Conversely, susceptibility to tick infestation in Brangus steers could be related to a maladaptive response due to decreased transcriptional levels of the previously discussed DEGs, while also reporting high expression of immune BoLA genes. Before tick exposure, the LR Brangus also reported increased expression of KLK1 (kallikrein 1) and various immunoglobulin-like genes such as LOC618463 (sialic acid-binding Ig-like lectin 13 / SIGLEC12), LOC112441499 (Vlambda1) and LOC100847119 (immunoglobulin lambda-1 light chain-like), as well as LOC524810 (IgM). KLK1 encodes a serine protease involved in the kallikrein-kinin system and has been implicated in proinflammatory activities but also assigned a protective role in skin wound healing [70]. While More et al. [16] reported KLK1 to be enriched in the skin of tick resistant Braford cattle, this study found downregulated expression of KLK1 in the leukocytes of resistant Brangus. The discrepancies in the direction of expression change could be due to differences in the tissues analysed and breed-related differences. Additionally, a higher expression of humoral immunity molecules in the skin of LR compared to HR cattle has also been reported by Wang et al. [27], however, it follows previous findings that humoral activity might not necessarily protect against tick infestation because it could contribute to immunopathology for the susceptible hosts [15].

Although much of the previous work on tick resistance was initially dedicated to the quantification of gene expression differences between B. indicus and B. taurus breeds, many of the conclusions are breed-specific, suggesting that genetic determinants of the bovine immune function still need to be properly characterised [71]. On the other hand, crossbred models are more likely to minimise breed-related effects, but as reported in this study, there is a wider variation of the individual responses to tick infestation suggesting that more experimental timepoints throughout the experimental trial should be evaluated. It can be appreciated from previous studies that exposure protocols need to be further standardised across studies and the initial exposure state of the reference herds adequately controlled for, i.e. tick-naïve vs. tick-exposed pre-treated with acaricides to determine the contribution of the innate and adaptive immune responses, as well as the many effector mechanisms linking the two. Although the sample herds were relatively small for this experimental approach, this study presents an advantageous opportunity in the recruitment of tick-naïve cattle for elucidating some of the biology behind host responses to tick infestation. While regulatory changes at the transcriptional level may be immediate and short-lived, particularly following tick attachment, the comparisons performed with the chosen timepoints have captured important information that could contribute to improving the experimental design of future artificial tick infestation trials. Lastly, future research should aim toward meta-analyses across tissues and breeds to identify and validate candidate genes for tick burden in large animal populations to promptly build reference populations for genomic-assisted selection to mitigate the impact of cattle ticks.

Conclusions

The present transcriptomic study is the first to evaluate the effect of tick infestation on Brangus cattle with variable phenotypes of host resistance to R. australis. It was found that steers exposed to R. australis ticks for periods of 3 and 12 weeks responded to infestation by expressing leukocyte genes related to chemotaxis, cytokine secretion, and inflammatory response which are represented in immune and inflammatory pathways. IL-17 pathway and cytokine-cytokine interaction pathway appeared to be relevant in protection and immunopathology to tick challenge. Several immune, tissue remodelling, and angiogenesis genes were detected as significantly changed between high and low host resistance steers before the initial infestation, suggesting these mechanisms are relevant in the development of resistance or susceptibility to tick challenge. However, as the timing after repeated infestation progressed, a lower number of differentially expressed genes were found between the two phenotypes. The altered expression of genes from the bovine MHC complex in highly resistant animals at pre- and post- infestation stages also supports the relevance of this genomic region for disease resilience. Overall, this study offers a resource of leukocyte gene expression data that could be further evaluated in studies of ectoparasite resistance and to develop improvements to tick resistance trait selection in composite breed cattle.

Methods

Animals

The study was approved by the Animal Ethics Unit at The University of Queensland (certificate number QAAFI/469/18). All methods were carried out in accordance with the relevant guidelines and regulations. All methods were performed in accordance with the ARRIVE guidelines (Animal Research: Reporting of In Vivo Experiments) [57]. Thirty Brangus steers (6–8 months old and average weight of 200 kg) with no previous exposure to the R. australis were sourced from a tick-free region in Australia (Morven, QLD) and transported to the University of Queensland’s Pinjarra Hills Beef Research Unit (Brisbane, QLD). Animals were vaccinated using the live trivalent tick fever vaccine containing Babesia bovis, Babesia bigemina and Anaplasma centrale as recommended by the manufacturer (Tick Fever Centre, Queensland Department of Agriculture & Fisheries, Wacol, Queensland) at 5 weeks prior to the commencement of the trial.

Cattle tick infestations

The artificial infestation method was used to assess the level of host resistance to ticks on each animal [14]. Approximately 10,000 R. australis larvae of the non-resistant field strain [72] were applied to the base of the tail and dispersed along the back of the animal with a brush (Fig. 7A-B). Tick infestations were undertaken once weekly in the morning (7:30–8:30 am) for 12 consecutive weeks between May and August 2019 (Australian fall/winter), with a total of 13 infestations performed. The average monthly temperature during this period was 23.8 ± 0.9 °C in the Brisbane region [73].

Fig. 7
figure 7

Experimental design of the study. A Representation of the experimental design timeline used in the artificial tick infestation trial undertaken with Brangus steers showing the frequency of larval applications (orange dots), tick scoring (green triangle), sampling (blue square), and chosen sampling timepoints for RNA sequencing (purple star). B Schematic of larval tick application during infestation. C Schematic of tick scoring method. Created with BioRender.com

Tick scoring

Tick burden developing from the initial and subsequent infestations was evaluated by the same observer. Animals were scored in a cattle crush on a weekly basis period starting at week 3 until week 15 post-initial infestation (T3-T15) as represented in Fig. 7A. Scoring on week 7 could not be performed to due to staff unavailability, whereas a blank score was given in instances where animal behaviour was risky during the observation. The tick scoring scale used for this Brangus herd was as follows: 1 = 0–50 ticks, 2 = 50–100 ticks, 3 = 100–200 ticks, 4 = 200–300 ticks, 5 = more than 300 ticks, blank = not scored. Similar studies have used similar scoring methods for tick infestation phenotyping [9, 74]. The scores represent an estimation of the number of female adult ticks (un/semi/fully engorged) found on one side of the animal’s body (Fig. 7C). The mean tick score (MTS) from T8 to T15 timepoint was calculated to rank animals and define host resistance phenotype groups for transcriptomic analyses. Five animals with low MTS (1 < MTS < 2) were classified into the high host resistance (HR) group, whereas five animals with high MTS (3 < MTS < 5) were classified as low host resistance (LR) group. MTS sample values are provided in Additional File 10.

Sample collection and leukocyte RNA extraction

Sample collection was performed before tick challenge (day 0/T0), at 6 hours post-initial infestation, and then at day 22 (T3), 44 (T6), 65 (T9) and 85 (T12) of the trial. For each animal, approximately 2–6 mL of blood were drawn from the jugular vein into a single BD Vacutainer® K2EDTA tube (BD, USA) and kept on ice until transported to the laboratory. Red blood cell lysis for leukocyte isolation was performed within 12–48 hours post-collection with a protocol adaptation from the miRNeasy mini kit guide “Appendix D” (QIAGEN, USA) using 1X RBC lysis buffer (168 mM NH4Cl; 10 mM KHCO3, 0.1 mM EDTA in 500 mL Milli-Q water, pH 7.4). The leukocyte pellets were homogenized with 1200 μL of QIAzol reagent (QIAGEN, USA) and stored at − 80 °C until further use.

Leukocyte RNA was isolated from animals classified as HR (n = 5) and LR (n = 5) from days 0, 21, and 85 post-initial infestation. Samples from these timepoints are referred to in this manuscript as HR/LR-T0, HR/LR-T3, and HR/LR-T12 datasets. In total, 30 samples were extracted with the miRNeasy mini kit (QIAGEN, USA) as per manufacturer’s instructions. RNA samples were treated with DNA-free™ DNAse (Life Technologies, USA) and RNA was quantified with the Nanodrop 2000 spectrophotometer (ThermoFisher, USA). RNA quality analysis (RIN) was performed with the 2100 Bioanalyzer Instrument (Agilent Technologies, USA) by the Institute for Molecular Biosciences Sequencing Facility in St. Lucia, Australia. RIN sample values are provided in Additional File 10.

RNA sequencing

Out of the 30 RNA samples prepared for RNA-sequencing, two samples were removed due to low RNA concentration. Hence, a total of 28 cDNA libraries were prepared with the TruSeq Stranded mRNA kit (Illumina, USA) and sequenced as 100 bp single-end reads in one flow cell lane on a NovaSeq 6000 sequencer (Illumina, USA). The Illumina bcl2fastq 2.20.0.422 pipeline was used to generate the sequence data. Library preparation and sequencing were performed by the Australian Genome Research Facility (AGRF) in Melbourne, Australia.

Genomic estimate of Bos indicus content (BIC)

Tail hair DNA was used to genotype all steers with the GeneSeek® Genomic Profiler™ Bovine 50 K (GGP 50 K) by Neogen Australasia in Gatton, Australia. Bos indicus content was calculated using the 35 K array by comparing each animal to a large purebred B. indicus dataset [75]. Using a GBLUP model, phenotypes were assigned as 1 for B. indicus and 0 if not and the effect of each SNP was back-solved [76]. Prediction equations for B. indicus purebred were then used to estimate BIC in the steers [75]. The Pearson correlation method was used to test the significance between MTS and BIC values obtained for 29 out of 30 Brangus samples. BIC values for DEG covariate models are supplied in Additional File 10.

Bioinformatics pipeline for differential gene expression analysis

The bioinformatics analyses were carried out in The University of Queensland High Performance Computing Cluster [77], Galaxy Australia Server [78], and RStudio [79]. Briefly, a read quality control was performed with FastQC software (version 0.11.4) [80]. Adapters and low quality reads were removed with Trimmomatic software (version 0.35) [81] using parameters for single-end reads including –phred 33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 and Illumina adapter sequence “AGATCGGAAGAGC”. Reads were aligned to the genome assembly Bos taurus ARS-UCD1.2 using HISAT2 (Galaxy version 2.1.0 + galaxy6) [82] with default parameters and read strandness set to “reverse”. Gene count data were generated with featureCounts (Galaxy Version 1.6.4) [83] using default parameters and strand information set to “reverse”. The reference genome and annotation file for bosTau9 were obtained from the UCSC genome browser database [84].

Gene count matrices consisting of 28,787 B. taurus genes as rows and samples as columns were input in RStudio [79]. The edgeR Bioconductor package [85] was used to perform all differential gene expression analyses. Gene filtering was performed with an expression threshold of CPM > 0.5 according to n samples, where n was determined by the condition with the lowest number of samples tested. Libraries were normalised with the function calcNormFactors which uses the trimmed mean of M-values (TMM) method [86]. Under the negative binomial model, the common and tagwise dispersion were calculated with the function estimateDisp. Gene expression models for the analysis of host response to tick infestation (Model 1) and host resistance phenotypes (Model 2:A-B) are listed in Table 5 and further described next.

Table 5 Gene expression models and explanatory variables for leukocyte data

Model 1: host response to tick infestation

Differentially expressed genes were detected between tick-exposed (T3 or T12, n = 10 each) and non-exposed (tick-naive) steers (T0, n = 8) with the quasi-likelihood F-test (edgeRQLF) and the linear model ~IW + MTS + BIC + RIN. Here, infestation week timepoint (IW) was fitted as a factor with two levels for comparisons of timepoint datasets as T3-vs-T0 and T12-vs-T0. The covariates of mean tick score (MTS), B. indicus content (BIC), sample RIN value (RIN) were implemented to account for additional sources of variation in the datasets. Significant differentially expressed genes (DEGs) were identified with a false discovery rate (FDR) threshold at FDR < 5%.

Model 2: host resistance phenotypes

Differentially expressed genes were detected between high (HR) and low (LR) host resistance phenotypes with the likelihood ratio test (edgeRLRT) using two types of linear models: Model 2A (~MTS + BIC + RIN) implemented for pre-infestation (T0) data, and Model 2B (~MTS + TPS + BIC + RIN) implemented for post-infestation data (T3 and T12).

Model 2A evaluated the comparison of the groups LR-T0 (n = 5) and HR-T0 (n = 3) according to mean tick score (MTS) fitted as a continuous variable while accounting for the effects of B. indicus content (BIC) and sample RIN value (RIN).

Model 2B was the same as model 2A but also included a covariate for timepoint tick score (TPS) to account for the effect of timepoint tick burden measured on the T3 and T12 datasets, respectively. The comparisons were performed between LR-T3 (n = 5) and HR-T3 (n = 5), and between LR-T12 (n = 5) and HR-T12 (n = 5) according to the MTS value. Significant DEGs were identified with a threshold at FDR < 5%.

Functional enrichment analysis

Over-representation analysis (ORA) of Gene Ontology (GO) terms and KEGG (Kyoto Encyclopedia Genes and Genomes Pathways) were performed with the clusterProfiler R package [87] using filtered gene lists (FDR < 0.05 and |logFC| > 1) from T3-vs-T0 and T12-vs-T0 datasets. Additional filtering of highly redundant parent-child terms in the ORA-GO output was applied with the simplify function with the default threshold value. Graphics were created with dot plot and category network plot functions from this package. Pathway graphs were rendered with the Pathview R package [88] for Bos taurus organism sourced from the KEGG database [40,41,42].

Availability of data and materials

All relevant data are included in the manuscript and its Supplementary files. The datasets have been deposited with links to BioProject accession number PRJNA802321 in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA802321).

Abbreviations

T0:

Pre-infestation timepoint

T3:

Week 3 post-initial infestation timepoint

T12:

Week 12 post-initial infestation timepoint

HR:

High host Resistance phenotype

LR:

Low host Resistance phenotype

DEG:

Differentially Expressed Gene

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia Genes and Genomes Pathways

FDR:

False Discovery Rate

FC:

Fold change

MTS:

Mean Tick Score

IW:

Infestation Week

TPS:

Timepoint Tick Score

RIN:

RNA Integrity Number

BIC:

Bos indicus content

CPM:

Counts per million

MHC:

Major Histocompatibility Complex

References

  1. Estrada-Peña A, Venzal JM, Nava S, Mangold A, Guglielmone AA, Labruna MB, et al. Reinstatement of Rhipicephalus (Boophilus) australis (Acari: Ixodidae) with redescription of the adult and larval stages. J Med Entomol. 2012;49(4):794–802.

    Article  PubMed  Google Scholar 

  2. White N, Sutherst RW, Hall N, Whish-Wilson P. The vulnerability of the Australian beef industry to impacts of the cattle tick (Boophilus microplus) under climate change. Clim Chang. 2003;61(1):157–90.

    Article  Google Scholar 

  3. Lew-Tabor AE, Rodriguez VM. A review of reverse vaccinology approaches for the development of vaccines against ticks and tick borne diseases. Ticks Tick Borne Dis. 2016;7(4):573–85.

    Article  CAS  PubMed  Google Scholar 

  4. Tabor AE, Ali A, Rehman G, Garcia G, Zangirolamo A, Malardo T, et al. Cattle tick Rhipicephalus microplus-host interface: a review of resistant and susceptible host responses. Front Cell Infect Microbiol. 2017;7(506):1–18.

    CAS  Google Scholar 

  5. Burrow HM, Mans BJ, Cardoso FF, Birkett MA, Kotze AC, Hayes BJ, et al. Towards a new phenotype for tick resistance in beef and dairy cattle: a review. Anim Prod Sci. 2019;59(8):1401–27.

    Article  CAS  Google Scholar 

  6. Rodriguez-Vivas RI, Jonsson NN, Bhushan C. Strategies for the control of Rhipicephalus microplus ticks in a world of conventional acaricide and macrocyclic lactone resistance. Parasitol Res. 2018;117(1):3–29.

    Article  PubMed  Google Scholar 

  7. Abbas RZ, Zaman MA, Colwell DD, Gilleard J, Iqbal Z. Acaricide resistance in cattle ticks and approaches to its management: the state of play. Vet Parasitol. 2014;203(1):6–20.

    Article  CAS  PubMed  Google Scholar 

  8. Jonsson NN, Mayer DG, Green PE. Possible risk factors on Queensland dairy farms for acaricide resistance in cattle tick (Boophilus microplus). Vet Parasitol. 2000;88(1):79–92.

    Article  CAS  PubMed  Google Scholar 

  9. Cardoso FF, Matika O, Djikeng A, Mapholi N, Burrow HM, Yokoo MJI, et al. Multiple country and breed genomic prediction of tick resistance in beef cattle. Front Immunol. 2021;12:620847.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Porto Neto LR, Jonsson NN, D’occhio MJ, Barendse W. Molecular genetic approaches for identifying the basis of variation in resistance to tick infestation in cattle. Vet Parasitol. 2011;180(3–4):165–72.

    Article  CAS  PubMed  Google Scholar 

  11. Kitsou C, Fikrig E, Pal U. Tick host immunity: vector immunomodulation and acquired tick resistance. Trends Immunol. 2021;42(7):554–74.

    Article  CAS  PubMed  Google Scholar 

  12. Park JM, Oliva Chávez AS, Shaw DK. Ticks: more than just a pathogen delivery service. Front Cell Infect Microbiol. 2021;11:739419.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Seifert GW. Variations between and within breeds of cattle in resistance to field infestations of the cattle tick ("Boophilus microplus"). Aust J Agric Res. 1971;22(1):159–68.

    Article  Google Scholar 

  14. Utech KB, Wharton RH, Kerr JD. Resistance to Boophilus microplus (Canestrini) in different breeds of cattle. Aust J Agric Res. 1978;29(4):885–95.

    Article  Google Scholar 

  15. Piper EK, Jonsson NN, Gondro C, Vance ME, Lew-Tabor A, Jackson LA. Peripheral cellular and humoral responses to infestation with the cattle tick Rhipicephalus microplus in Santa Gertrudis cattle. Parasite Immunol. 2017;39(e12402):1–12.

    Google Scholar 

  16. Moré DD, Cardoso FF, Mudadu MA, Malagó-Jr W, Gulias-Gomes CC, Sollero BP, et al. Network analysis uncovers putative genes affecting resistance to tick infestation in Braford cattle skin. BMC Genomics. 2019;20(1):998.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Andreotti R, Barros JC, Garcia MV, Rodrigues VdS, Higa LdOS, Duarte PO, et al. Cattle tick infestation in Brangus cattle raised with Nellore in central Brazil. Semina: Ciênc Agrár. 2018;39(3):1099–114.

  18. Franzin AM, Maruyama SR, Garcia GR, Oliveira RP, Ribeiro JMC, Bishop R, et al. Immune and biochemical responses in skin differ between bovine hosts genetically susceptible and resistant to the cattle tick Rhipicephalus microplus. Parasit Vectors. 2017;10(1):51.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  19. Piper EK, Jackson LA, Bagnall NH, Kongsuwan KK, Lew AE, Jonsson NN. Gene expression in the skin of Bos taurus and Bos indicus cattle infested with the cattle tick, Rhipicephalus (Boophilus) microplus. Vet Immunol Immunopathol. 2008;126(1):110–9.

    Article  CAS  PubMed  Google Scholar 

  20. Piper EK, Jackson LA, Bielefeldt-Ohmann H, Gondro C, Lew-Tabor AE, Jonsson NN. Tick-susceptible Bos taurus cattle display an increased cellular response at the site of larval Rhipicephalus (Boophilus) microplus attachment, compared with tick-resistant Bos indicus cattle. Int J Parasitol. 2010;40(4):431–41.

    Article  PubMed  Google Scholar 

  21. Marima JK, Nel CL, Marufu MC, Jonsson NN, Dube B, Dzama K. A genetic and immunological comparison of tick-resistance in beef cattle following artificial infestation with Rhipicephalus ticks. Exp Appl Acarol. 2020;80(4):569–90.

    Article  CAS  PubMed  Google Scholar 

  22. Porto Neto LR, Bunch RJ, Harrison BE, Barendse W. DNA variation in the gene ELTD1 is associated with tick burden in cattle. Anim Genet. 2010;42(1):50–5.

    Article  PubMed  CAS  Google Scholar 

  23. Porto Neto LR, Bunch RJ, Harrison BE, Prayaga KC, Barendse W. Haplotypes that include the integrin alpha 11 gene are associated with tick burden in cattle. BMC Genet. 2010;11(1):55.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Porto Neto LR, Jonsson NN, Ingham A, Bunch RJ, Harrison BE, Barendse W. The RIPK2 gene: a positional candidate for tick burden supported by genetic associations in cattle and immunological response of knockout mouse. Immunogenetics. 2012;64(5):379–88.

    Article  CAS  PubMed  Google Scholar 

  25. Porto Neto LR, Piper EK, Jonsson NN, Barendse W, Gondro C. Meta-analysis of genome wide association and gene expression studies to identify candidate genes for tick burden in cattle. In: Proceedings of the 9th World Congress on Genetics Applied to Livestock Production. Leipzig: German Society for Animal Science; 2010. p. 0664.

  26. Robbertse L, Richards SA, Stutzer C, Olivier NA, Leisewitz AL, Crafford JE, et al. Temporal analysis of the bovine lymph node transcriptome during cattle tick (Rhipicephalus microplus) infestation. Vaccine. 2020;38(44):6889–98.

    Article  CAS  PubMed  Google Scholar 

  27. Wang YH, Reverter A, Kemp D, McWilliam SM, Ingham A, Davis CA, et al. Gene expression profiling of Hereford shorthorn cattle following challenge with Boophilus microplus tick larvae. Aust J Exp Agric. 2007;47(12):1397–407.

    Article  CAS  Google Scholar 

  28. Constantinoiu CC, Jackson LA, Jorgensen WK, Lew-Tabor AE, Piper EK, Mayer DG, et al. Local immune response against larvae of Rhipicephalus ( Boophilus) microplus in Bos taurus indicus and Bos taurus taurus cattle. Int J Parasitol. 2010;40(7):865–75.

    Article  CAS  PubMed  Google Scholar 

  29. Constantinoiu CC, Jonsson NN, Jorgensen WK, Jackson LA, Piper EK, Lew-Tabor AE. Immuno-fluorescence staining patterns of leukocyte subsets in the skin of taurine and indicine cattle. Res Vet Sci. 2013;95(3):854–60.

    Article  CAS  PubMed  Google Scholar 

  30. Constantinoiu CC, Lew-Tabor A, Jackson LA, Jorgensen WK, Piper EK, Mayer DG, et al. Local immune response to larvae of Rhipicephalus microplus in Santa Gertrudis cattle. Parasite Immunol. 2018;40(4):e12515.

    Article  CAS  PubMed  Google Scholar 

  31. Piper EK, Jonsson NN, Gondro C, Lew-Tabor AE, Moolhuijzen P, Vance ME, et al. Immunological profiles of Bos taurus and Bos indicus cattle infested with the cattle tick, Rhipicephalus (Boophilus) microplus. Clin Vaccine Immunol. 2009;16(7):1074–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kongsuwan K, Josh P, Colgrave ML, Bagnall NH, Gough J, Burns B, et al. Activation of several key components of the epidermal differentiation pathway in cattle following infestation with the cattle tick, Rhipicephalus ( Boophilus) microplus. Int J Parasitol. 2010;40(4):499–507.

    Article  PubMed  Google Scholar 

  33. Marufu MC, Dzama K, Chimonyo M. Cellular responses to Rhipicephalus microplus infestations in pre-sensitised cattle with differing phenotypes of infestation. Exp Appl Acarol. 2014;62(2):241–52.

    Article  PubMed  Google Scholar 

  34. Robbertse L, Richards SA, Clift SJ, Barnard A-C, Leisewitz A, Crafford JE, et al. Comparison of the differential regulation of T and B-lymphocyte subsets in the skin and lymph nodes amongst three cattle breeds as potential mediators of immune-resistance to Rhipicephalus microplus. Ticks Tick Borne Dis. 2018;9(4):976–87.

    Article  PubMed  Google Scholar 

  35. Turni C, Lee RP, Jackson LA. Effect of salivary gland extracts from the tick, Boophilus microplus , on leucocytes from Brahman and Hereford cattle. Parasite Immunol. 2002;24(7):355–61.

    Article  PubMed  Google Scholar 

  36. Kongsuwan K, Piper EK, Bagnall NH, Ryan K, Moolhuijzen P, Bellgard M, et al. Identification of genes involved with tick infestation in Bos taurus and Bos indicus. Dev Biol. 2008;132:77–88.

    CAS  Google Scholar 

  37. Araújo W, Domingues R, de Azevedo MC, da Silva MVGB, Correa G, Facioni SE, et al. Microarray analysis of tick-infested skin in resistant and susceptible cattle confirms the role of inflammatory pathways in immune activation and larval rejection. Vet Parasitol. 2014;205(1–2):307–17.

    Google Scholar 

  38. Maruyama SR, Carvalho B, González-Porta M, Rung J, Brazma A, Gustavo Gardinassi L, et al. Blood transcriptome profile induced by an efficacious vaccine formulated with salivary antigens from cattle ticks. NPJ Vaccines. 2019;4(1):53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Chen KW, Demarco B, Ramos S, Heilig R, Goris M, Grayczyk JP, et al. RIPK1 activates distinct gasdermins in macrophages and neutrophils upon pathogen blockade of innate immune signaling. PNAS. 2021;118(28):e2101189118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.

    Article  CAS  PubMed  Google Scholar 

  43. Jonsson NN, Piper EK, Constantinoiu CC. Host resistance in cattle to infestation with the cattle tick Rhipicephalus microplus. Parasite Immunol. 2014;36(11):553–9.

    Article  CAS  PubMed  Google Scholar 

  44. Junior PB, Rodrigues VS, Garcia MV, Higa LOS, Zimmermann NP, Barros JC, et al. Economic performance evaluation of Brangus and Nellore cattle breed naturally infested with Rhipicephalus microplus in an extensive production system in central-West Brazil. Exp Appl Acarol. 2019;78(4):565–77.

    Article  PubMed  Google Scholar 

  45. Paim TP, Hay EHA, Wilson C, Thomas MG, Kuehn LA, Paiva SR, et al. Genomic breed composition of selection signatures in Brangus beef cattle. Front Genet. 2020;11:710.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Goszczynski DE, Corbi-Botto CM, Durand HM, Rogberg-Muñoz A, Munilla S, Peral-Garcia P, et al. Evidence of positive selection towards Zebuine haplotypes in the BoLA region of Brangus cattle. Animal. 2018;12(2):215–23.

    Article  CAS  PubMed  Google Scholar 

  47. Viney ME, Riley EM, Buchanan KL. Optimal immune responses: immunocompetence revisited. Trends Ecol Evol. 2005;20(12):665–9.

    Article  PubMed  Google Scholar 

  48. Glass E. The molecular pathways underlying host resistance and tolerance to pathogens. Front Genet. 2012;3:263.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Rauw WM. Immune response from a resource allocation perspective. Front Genet. 2012;3:267.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Kelly MN, Kolls JK, Happel K, Schwartzman JD, Schwarzenberger P, Combe C, et al. Interleukin-17/Interleukin-17 receptor-mediated signaling is important for generation of an optimal polymorphonuclear response against toxoplasma gondii infection. Infect Immun. 2005;73(1):617–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Heinze DM, Wikel SK, Thangamani S, Alarcon-Chaidez FJ. Transcriptional profiling of the murine cutaneous response during initial and subsequent infestations with Ixodes scapularis nymphs. Parasit Vectors. 2012;5(1):26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Sun T, Wang F, Pan W, Wu Q, Wang J, Dai J. An immunosuppressive tick salivary gland protein DsCystatin interferes with toll-like receptor signaling by downregulating TRAF6. Front Immunol. 2018;9:1245.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Qian Y, Kang Z, Liu C, Li X. IL-17 signaling in host defense and inflammatory diseases. Cell Mol Immunol. 2010;7(5):328–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Sallusto F, Lanzavecchia A, Mackay CR. Chemokines and chemokine receptors in T-cell priming and Th1/Th2-mediated responses. Immunol Today. 1998;19(12):568–74.

    Article  CAS  PubMed  Google Scholar 

  55. Groom JR, Luster AD. CXCR3 in T cell function. Exp Cell Res. 2011;317(5):620–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Domingues R, Wohlres-Viana S, Reis D, Teixeira H, Ferreira A, Guimarães S, et al. Expression of immune response genes in peripheral blood of cattle infested with Rhipicephalus microplus. Genet Mol Res. 2014;13(2):4013–21.

    Article  CAS  PubMed  Google Scholar 

  57. Raza A, Schulz BL, Nouwens A, Jackson LA, Piper EK, James P, et al. Serum proteomes of Santa Gertrudis cattle before and after infestation with Rhipicephalus australis ticks. Parasite Immunol. 2021;43(7):e12836–n/a.

    Article  CAS  PubMed  Google Scholar 

  58. Behl JD, Verma NK, Tyagi N, Mishra P, Behl R, Joshi BK. The major histocompatibility complex in bovines: a review. ISRN Vet Sci. 2012;2012:872710.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Martinez ML, Machado MA, Nascimento CS, Silva MVGB, Teodoro RL, Furlong J, et al. Association of BoLA-DRB3.2 alleles with tick (Boophilus microplus) resistance in cattle. Genet Mol Res. 2006;5(3):513–24.

    CAS  PubMed  Google Scholar 

  60. Machado MA, Azevedo ALS, Teodoro RL, Pires MA, Peixoto MGC, de Freitas C, et al. Genome wide scan for quantitative trait loci affecting tick resistance in cattle (Bos taurus × Bos indicus). BMC Genomics. 2010;11(280):1–11.

    Google Scholar 

  61. Otto PI, Guimarães SEF, Verardo LL, Azevedo ALS, Vandenplas J, Soares ACC, et al. Genome-wide association studies for tick resistance in Bos taurus × Bos indicus crossbred cattle: a deeper look into this intricate mechanism. J Dairy Sci. 2018;101(12):11020–32.

    Article  CAS  PubMed  Google Scholar 

  62. Acosta-Rodríguez R, Alonso-Morales R, Balladares S, Flores-Aguilar H, García-Vazquez Z, Gorodezky C. Analysis of BoLA class II microsatellites in cattle infested with Boophilus microplus ticks: class II is probably associated with susceptibility. Vet Parasitol. 2005;127(3):313–21.

    Article  PubMed  CAS  Google Scholar 

  63. Charo IF, Ransohoff RM. The many roles of chemokines and chemokine receptors in inflammation. N Engl J Med. 2006;354(6):610–21.

    Article  CAS  PubMed  Google Scholar 

  64. Bonetti MI. Microfibrils: a cornerstone of extracellular matrix and a key to understand Marfan syndrome. Ital J Anat Embryol. 2009;114(4):201–24.

    PubMed  Google Scholar 

  65. Sato K, Yomogida K, Wada T, Yorihuzi T, Nishimune Y, Hosokawa N, et al. Type XXVI collagen, a new member of the collagen family, is specifically expressed in the testis and ovary. J Biol Chem. 2002;277(40):37678–84.

    Article  CAS  PubMed  Google Scholar 

  66. Ross FP, Chappel J, Alvarez J, Sander D, Butler W, Farach-Carson M, et al. Interactions between the bone matrix proteins osteopontin and bone sialoprotein and the osteoclast integrin alpha v beta 3 potentiate bone resorption. J Biol Chem. 1993;268(13):9901–7.

    Article  CAS  PubMed  Google Scholar 

  67. Bagnall N, Gough J, Cadogan L, Burns B, Kongsuwan K. Expression of intracellular calcium signalling genes in cattle skin during tick infestation. Parasite Immunol. 2009;31(4):177–87.

    Article  CAS  PubMed  Google Scholar 

  68. DeHahn KC, Gonzales M, Gonzalez AM, Hopkinson SB, Chandel NS, Brunelle JK, et al. The alpha4 laminin subunit regulates endothelial cell survival. Exp Cell Res. 2004;294(1):281–9.

    Article  CAS  PubMed  Google Scholar 

  69. Gonzalez AM, Gonzales M, Herron GS, Nagavarapu U, Hopkinson SB, Tsuruta D, et al. Complex interactions between the laminin α4 subunit and integrins regulate endothelial cell behavior in vitro and angiogenesis in vivo. PNAS. 2002;99(25):16075.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Chao J, Shen B, Gao L, Xia C-F, Bledsoe G, Chao L. Tissue kallikrein in cardiovascular, cerebrovascular and renal diseases and skin wound healing. Biol Chem. 2010;391(4):345–55.

    Article  CAS  PubMed  Google Scholar 

  71. Robbertse L, Richards SA, Maritz-Olivier C. Bovine immune factors underlying tick resistance: integration and future directions. Front Cell Infect Microbiol. 2017;7:522.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. Stewart NP, Callow LL, Duncalfe F. Biological comparisons between a laboratory-maintained and a recently isolated field strain of Boophilus microplus. J Parasitol. 1982;68(4):691–4.

    Article  CAS  PubMed  Google Scholar 

  73. Australian Bureau of Meteorology. https://www.bom.gov.au/climate/data/index.shtml. Accessed 11 Apr 2022.

  74. Abeysinghe P, Turner N, Peiris H, Vaswani K, Cameron N, McGhee N, et al. Differentially expressed extracellular vesicle, exosome and non-exosome miRNA profile in high and low tick-resistant beef cattle. Front Cell Infect Microbiol. 2021;11:780424.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Hayes B, Fordyce G, S L. Genomic predictions for fertility traits in tropical beef cattle from a multi-breed, crossbred and composite reference population. In: Proceedings of the 23rd Conference of the Association for the Advancement of Animal Breeding and Genetics (AAABG). Armidale: Animal Production Science; 2019. p. 282–85.

  76. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. QRIScloud University of Queensland. https://www.qriscloud.org.au. Accessed 30 Oct 2019.

  78. Galaxy Australia. https://usegalaxy.org.au. Accessed 20 Nov 2019.

  79. RStudio team. RStudio: integrated development for R. Boston, USA. 2020. http://www.rstudio.com/. Accessed 25 Aug 2020.

  80. Andrews S. FastQC A quality control tool for high throughput sequence data. Babraham Bioinformatics. http://www.bioinformatics.babraham.ac.uk/projects/fastqc. Accessed 30 Oct 2019.

  81. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2013;30(7):923–30.

    Article  PubMed  CAS  Google Scholar 

  84. UCSC Genome Browser. ftp://hgdownload.soe.ucsc.edu/goldenPath/bosTau9//bigZips. Accessed 2 Mar 2020.

  85. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009;26(1):139–40.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  86. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  87. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.

    PubMed  PubMed Central  Google Scholar 

  88. Luo W, Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. 2013;29(14):1830–1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgments

The authors thank Tom Conolly and Alison Moore for assistance with herd management and tick infestations at the UQ Pinjarra Hills Farm. The authors also acknowledge Dr. Tim Bruxner from the UQ Institute for Molecular Bioscience for providing services in the Bioanalyser 2100 Instrument, as well as Dr. Russell Lyons and Dr. Karen Schutt from Neogen Australasia for providing animal genotyping services.

Funding

This research has been funded by Meat and Livestock Australia Donor Company (MDC) project P.PSH.0798. EFMV acknowledges the Australian Government Research Training Program (RTP) Scholarship for financial support.

Author information

Authors and Affiliations

Authors

Contributions

EFMV collected and analysed all transcriptomic data and wrote the manuscript. ER, BJ and NJ provided support in model development and data analyses; NJ and AET assisted with the interpretation of results. AR performed tick scoring. BJ performed calculations for genomic estimates of Bos indicus content. AR, MK, MNN and EFMV undertook field work and sample collection. PJ and AET lead the associated grant. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Emily F. Mantilla Valdivieso, Nicholas N. Jonsson or Ala E. Tabor.

Ethics declarations

Ethics approval and consent to participate

The study was conducted under approval of The University of Queensland Animal Ethics for project “Biomarkers and accurate phenotyping for selecting resistance to ectoparasites” with certificate number QAAFI/469/18. All methods were carried out in accordance with the relevant guidelines and regulations. All methods were reported in accordance with the ARRIVE guidelines (Animal Research: Reporting of In Vivo Experiments).

Consent for publication

Not Applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mantilla Valdivieso, E.F., Ross, E.M., Raza, A. et al. Transcriptional changes in the peripheral blood leukocytes from Brangus cattle before and after tick challenge with Rhipicephalus australis. BMC Genomics 23, 454 (2022). https://doi.org/10.1186/s12864-022-08686-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08686-3

Keywords