Expression network analysis of bovine skin infested with Rhipicephalus australis identifies pro-inflammatory genes contributing to tick susceptibility

Abstract

The skin is the primary feeding site of ticks that infest livestock animals such as cattle. The highly specialised functions of skin at the molecular level may be a factor contributing to variation in susceptibility to tick infestation; but these remain to be well defined. The aim of this study was to investigate the bovine skin transcriptomic profiles of tick-naïve and tick-infested cattle and to uncover the gene expression networks that influence contrasting phenotypes of host resistance to ticks. RNA-Seq data was obtained from skin of Brangus cattle with high (n = 5) and low (n = 6) host resistance at 0 and 12 weeks following artificial tick challenge with Rhipicephalus australis larvae. No differentially expressed genes were detected pre-infestation between high and low resistance groups, but at 12-weeks there were 229 differentially expressed genes (DEGs; FDR < 0.05), of which 212 were the target of at least 1866 transcription factors (TFs) expressed in skin. Regulatory impact factor (RIF) analysis identified 158 significant TFs (P < 0.05) of which GRHL3, and DTX1 were also DEGs in the experiment. Gene term enrichment showed the significant TFs and DEGs were enriched in processes related to immune response and biological pathways related to host response to infectious diseases. Interferon Type 1-stimulated genes, including MX2, ISG15, MX1, OAS2 were upregulated in low host resistance steers after repeated tick challenge, suggesting dysregulated wound healing and chronic inflammatory skin processes contributing to host susceptibility to ticks. The present study provides an assessment of the bovine skin transcriptome before and after repeated tick challenge and shows that the up-regulation of pro-inflammatory genes is a prominent feature in the skin of tick-susceptible animals. In addition, the identification of transcription factors with high regulatory impact provides insights into the potentially meaningful gene-gene interactions involved in the variation of phenotypes of bovine host resistance to ticks.

Keywords: Bovine; Cattle; Gene expression; Host resistance; RNA-Seq.

Figure 1

Schematic representation of the bioinformatic pipeline implemented in the analysis of total RNA-Seq data from skin of Brangus steers before and after tick challenge. Shaded boxes represent software packages. The timepoint tick scores, mean tick scores, RNA integrity numbers (RIN), and Bos indicus content (BIC) values for the selected high and low resistant animals are summarised in Supplementary File 1.

Figure 2

Differentially expressed genes (DEGs) in skin from tick-infested compared to tick-naïve steers. (A). Multidimensional scaling (MDS) plot shows clustering of samples according to sampling timepoint (T0: tick-naive; T12: week12 post-initial infestation). (B) Volcano plot of DEGs where x-axis represents fold change and y-axis represents statistical significance. Gene symbols are shown for DEGs with fold change |logFC|> 2 and FDR < 0.05. (C) Heatmap of the top 50 DEGs (25 most upregulated and 25 most downregulated) showing gene expression levels (logCPM) across samples from each timepoint.

Figure 3

Differentially expressed genes in skin samples of low compared to high host resistance steers. (A) MDS plot shows clustering of samples according to phenotype within sampling timepoint (T0: tick-naive; T12: week12 post-initial infestation). (B) Volcano plot of DEGs in T12 comparison where x-axis represents fold change and y-axis represents statistical significance. Symbols are shown for DEGs with fold change |logFC|> 2 and FDR < 0.05. (C) Heatmap of the top 50 DEGs (25 most upregulated and 25 most downregulated) showing gene expression levels (logCPM) across samples from each timepoint.

Figure 4

Distribution of significant DEGs and TFs in skin of Brangus steers and their functional enrichment. Venn diagram showing the distribution of genes identified as DEG, expressed TFs, and significant TFs (from RIF analysis) in the comparison of (A) tick-infested vs. tick-naïve steers (T12-vs-T0) and (B) high vs. low host resistance steers at T12 (LR_T12-vs-HR_T12). Enriched GO biological process terms in skin DEGs and TFs (FDR < 0.05, RIF metrics > 1.96, logFC > 1) from (C) T12-vs-T0 timepoint comparison and (D) LR_T12-vs-HR_T12 phenotype comparison. Enriched KEGG pathways in skin DEGs and TFs (FDR < 0.05, RIF metrics > 1.96) from (E) T12-vs-T0 timepoint comparison and (F) LR_T12-vs-HR_T12 phenotype comparison. The dot colour represents significance of the term (p-adjusted < 0.05), and dot size (GeneRatio) represents the number of significant genes in that pathway.