Blood transcriptome profile induced by an efficacious vaccine formulated with salivary antigens from cattle ticks

Ticks cause massive damage to livestock and vaccines are one sustainable alternative for the acaricide poisons currently heavily used to control infestations. An experimental vaccine adjuvanted with alum and composed by four recombinant salivary antigens mined with reverse vaccinology from a transcriptome of salivary glands from Rhipicephalus microplus ticks was previously shown to present an overall efficacy of 73.2% and cause a significant decrease of tick loads in artificially tick-infested, immunized heifers; this decrease was accompanied by increased levels of antigen-specific IgG1 and IgG2 antibodies, which were boosted during a challenge infestation. In order to gain insights into the systemic effects induced by the vaccine and by the tick challenge we now report the gene expression profile of these hosts’ whole-blood leukocytes with RNA-seq followed by functional analyses. These analyses show that vaccination induced unique responses to infestations; genes upregulated in the comparisons were enriched for processes associated with chemotaxis, cell adhesion, T-cell responses and wound repair. Blood transcriptional modules were enriched for activation of dendritic cells, cell cycle, phosphatidylinositol signaling, and platelets. Together, the results indicate that by neutralizing the tick’s salivary mediators of parasitism with vaccine-induced antibodies, the bovine host is able to mount normal homeostatic responses that hinder tick attachment and haematophagy and that the tick otherwise suppresses with its saliva.


INTRODUCTION
Rhipicephalus microplus, the cattle tick, causes heavy infestations in taurine breeds of cattle, severely affecting their health and causing huge production losses where they are employed. 1 Vaccines represent an alternative to acaricides to control tick infestations. TickGard 2 and GAVAC 3 are the vaccines to control cattle ticks that went to market. Both are based on a single antigen, the tick gut glycoprotein Bm86. The efficacy of both vaccines is variable [4][5][6][7] and since the antigen is hidden from the hosts immune system during parasitism, the memory they induce is short-lived, thus requiring at least one boost a year, 8,9 which makes them unattractive to producers. An approach to obtain a level of vaccine efficacy compatible with ease of management and good production in the field consists of formulating a multicomponent vaccine combining many antigens, in an attempt to decrease the proportion of poor responders to the vaccine and to also target the tick's multiple mediators of parasitism. In addition, since these mediators are mostly in the saliva that ticks inoculate in their hosts, memory may be boosted when cattle are naturally exposed to ticks.
In a previous study, we tested an experimental vaccine formulation composed by four recombinant salivary antigens from R. microplus, named Rm39, Rm76, Rm180, Rm239, which were predicted in silico as a glycine-rich cement protein, an immunoglobulin binding-protein, a thrombin inhibitor, and a metalloprotease, respectively. These antigens were chosen due to their roles in parasitism of the host: attachment of the tick (Rm39), abrogation of antibody responses (Rm76), anticoagulation (Rm180) and disruption of the extracellular matrix to create the blood-feeding pool (Rm239). The vaccine caused a significant decrease of tick loads in artificially tick-infested, immunized heifers; this decrease was accompanied by increased levels of antigen-specific IgG1 and IgG2 antibodies, including during a challenge infestation. 10 In order to gain insight into the mechanisms that must be induced by an efficacious vaccine and to identify potential biomarkers of efficacy, we examined the transcriptional profile elicited in blood leukocytes by this multicomponent tick vaccine. For this, we performed an RNA-seq analysis of blood cells from vaccinated and control animals of a genetically tick-susceptible breed, the target for such a vaccine. The present study presents a bovine transcriptome sequencing data obtained in the context of Systems Vaccinology 11,12 applied to understanding effector mechanisms of immunological-based control of cattle ticks.

RESULTS
Characteristics of the blood transcriptome and most significant differentially expressed genes Details of the immunization scheme and challenge with ticks are reported in our previous study 10 and depicted in Fig. 1. Blood samples from that experiment were employed in this study to perform RNA-seq analysis; 24 blood samples were collected from calves in the control (adjuvant only) and vaccinated groups (n = 4 in each group) at three different time points: before vaccination (BA, adjuvant only group; BV, vaccinated group), on day 7 after the third and last vaccine dose (AA, adjuvant only group; AV, vaccinated group) and on day 17 after a challenge with tick larvae, when they have molted to adults and the heifers have thus been exposed to most of the parasite's salivary antigens, including those in the vaccine (CHA, adjuvant only group; CHV, vaccinated group).
The 24 Illumina RNA-seq libraries were sequenced in four lanes with an average of 31 million single-end reads and 3.1 Gb per library. After quality assessment of the libraries, the reads were mapped to the bovine genome and quantified at the gene level. Subsequently, the differentially expressed genes (DEGs) across the experimental conditions were determined. In total 13,952 genes were expressed across all 24 samples, presenting a low biological coefficient of variation (BCV = 16.2%) among the biological replicates.
The DEGs were calculated in a comparative analyses to answer the following questions: (a) which genes respond to vaccination (i.  Table 1. Some of these DEGs, such as TIEG2 (Krueppel-like factor 11) and BT.64205 (antigen WC1.1 precursor, also named BoWC1.1, WC1 isolate CH149, CD163 molecule-like 1), were differentially expressed in two or more comparisons. Many uncharacterised proteins were highly differentially expressed. Other possible comparisons were also performed: BA vs AA, CHA vs AA and CHV vs AV, resulting in 985, 79 and 70 differentially expressed genes, respectively. All DEGs are described in 'Supplementary Data 1, sheets a-dDEG'.
Profiles of DEGs according to the experimental conditions To explore how the patterns of DEGs are associated across the groups during different time points of the experiment, hierarchical clustering (HCL) analyses were performed and the results were plotted as a heatmap (Fig. 2). To this end, RPKM values of the DEG lists were obtained through the normalisation of read counts according to gene length. Thus, we compared the levels of gene expression between genes. The results obtained in the heatmap showed that the expression profile of DEGs clustered according to the experimental conditions. Interestingly, the 74 DEGs identified in CHV vs. CHA pairwise comparisons presented an expression profile in which the infested vaccinated group (CHV) was clustered in an isolated position in the heatmap (Fig. 2 right group). These HCL analyses suggest that the multicomponent vaccine induced a unique transcriptional profile regulation of blood leukocytes from bovines responding to tick infestations. The identified DEGs will be useful for the elucidation of the molecular mechanisms involved in the response induced in vaccinated animals that reduces tick infestations.
Functional analyses of transcriptome To identify the biological pathways, networks and GO (Gene Ontology) processes associated with the genes affected after vaccination and tick infestation, the MetaCore™ platform for functional genomics data from Clarivate Analytics was used. The bovine gene IDs were converted in orthologous human gene IDs using BioMart. 13 Approximately 90% of bovine genes (from DEG lists) demonstrated homology with human genes. From the DEG list of AV, CHV and CHA vs. CHV comparisons, 380, 152 and 71 human orthologues, respectively, were retrieved from BioMart. The functional analyses of the genes affected after vaccination and tick infestation were enriched for several biological functions associated with immune responses (AV comparisons, CHV comparisons and CHV vs. CHA comparisons, Tables 2-4, respectively). The genes upregulated in the three comparisons were enriched for processes associated with chemotaxis, cell adhesion, the T-cell immune response, and organ regeneration, this latter process most likely associated with wound healing of the skin. The downregulated genes were enriched in other categories of immune processes, such as inflammation through the complement system and interferon signalling, and cell cycle processes.
Regarding pathway analysis in Metacore outputs, the pathway maps are graphic images drawn to represent biochemical pathways or signalling cascades; the genes in the dataset that are present in the given pathway map are indicated in a thermometer pictogram that is either in blue (negative logFC,  downregulated gene) or red (positive logFC, upregulated gene). The biological pathway "Role of IL-8 angiogenesis" (Fig. 3) was highly enriched in the CHV comparison. Another functional insight comes from the outputs generated by Metacore for networks, which comprise graphical representations displaying objects (molecular entities) connected with other molecules in reactions, processes or relationships, each having a number of biological attributes with symbolic representations in the outputs. The "Network Statistics" feature in Metacore calculates the most prevalent networks involved in a given dataset, combining objects from different networks in a unique chart. Network objects (genes) presented in the dataset (AV and CHV comparisons) are indicated as either a blue (negative logFC, downregulated gene), red (positive logFC, upregulated gene) or mixed-coloured (when positive and negative values are present for the same object) circle. The two most significant network compilation images were generated through the "Network Statistics" analysis using the DEGs from AV and CHV comparison-enriched processes; these networks pertained to regulation of the immune response, response to wounding and leukocyte chemotaxis (Fig. 4a, b).

Evaluation of simultaneous increases and decreases of expression of sets of genes
To gain further insights about the transcriptional responses to vaccination and tick challenge, we performed gene set enrichment analysis (GSEA) using the framework of Blood Transcription Modules (BTM) as gene sets. 14 GSEA accounts for simultaneously increased and decreased expression of a set of genes in a given module, the members of which may not have been identified as differentially expressed at the level of single-gene analysis. Over 90 BTMs were significantly associated with the entire dataset for at least one statistical comparison ('Supplementary Data 1, sheet eBTMs'). BTMs were grouped into high-level annotation based on common pathways or cell lineages to facilitate the interpretation and visualization (Fig. 5). Compared to baseline gene expression, adjuvant injection induced higher expression of genes related to extracellular matrix, intracellular organelles, B cells and T cells, while downregulated genes were related to cell cycle and division, type I interferon response, inflammatory signaling, dendritic cells (DCs), monocytes, neutrophils and platelets. Vaccination also induced higher expression of genes related to B and T cells, but not extracellular matrix. Still compared to baseline levels, challenge of control group induced higher expression of genes associated with cell cycle and division, intracellular organelles and platelets, whereas genes related to T cells were downregulated. Genes from vaccinated and challenged animals were only negatively associated with high-level annotation BTMs when compared to baseline levels. These include BTMs related to cell cycle and division, type I interferon, inflammatory signaling and DCs. Of note, comparisons between experimental groups, CHA vs AA or CHV vs AV, using log CPM values resulted in only seven significant associations with BTMs. The lack of statistical significance in GSEA suggested high individual variability in gene expression between animals composing the same experimental groups. Therefore, to remove individual confounding effects, we used baseline-subtracted log CPM values to perform GSEA, as previously described 15 Tick challenge after adjuvant injection induced positive associations with most high-level annotation BTMs, while genes related to inflammatory signaling were downregulated. Furthermore, tick challenge after vaccination induced higher expression of genes related to cell cycle and division, inflammatory signaling, intracellular organelles, monocytes and platelets, while negative associations included BTMs for

DISCUSSION
The RNA-seq analysis of whole blood from vaccinated calves indicated sets of genes that were differentially expressed and likely involved in the molecular mechanisms elicited during vaccination and tick infestation. Several genes associated with immune responses were differentially expressed in vaccinated calves. The transcriptional profile of 74 DEGs from the CHV vs. CHA comparison showed the efficient segregation of the vaccinated group during tick infestation, suggesting that this set of genes represents a transcriptional signature for protective vaccination against ticks. 10 Several genes presented immune-related functions, such as WC1, CD7, BLK, CXCL10, CXCL12 and CCR10. They are functionally highlighted and discussed below (Fig. 2).
The bovine genome contains 13 members of the WC1 gene family, organised into two loci on chromosome 5. 16 Two of these members (ENSBTAG00000039470 and ENSBTAG00000026531) were significantly upregulated (FDR < 0.01) in the CHV group. The WC1 gene family encodes co-receptors exclusively detected in gammadelta T cells, which play important roles in immune responses to develop recall responses to antigens. 17 CD7 is another T-cell receptor transcriptionally upregulated in the CHV group. This protein is expressed in NK cells and mature CD4+ and CD8+ peripheral T cells, playing an essential role in T-cell and T-cell/B-cell interactions during early lymphoid development. 18 Expression of the gene encoding the chemokine CXCL10 was downregulated in the CHV group. In two other studies that evaluated correlates of protective immune responses to inactivated influenza vaccines with a systems vaccinology approach, levels of serum CXCL10 correlated with good humoral responders, 19,20 i.e. production of neutralizing antibodies. Two other genes encoding proteins involved in chemotaxis, CXCL12 and CCR10, were upregulated. The CXCL12 chemokine, also known as stromal cell-derived factor 1 (SDF-1), is a strong chemoattractant for lymphocyte migration from the blood stream into sites of inflammation and the efficient activation of these cells under proinflammatory environments. 21 CCR10 is a chemokine receptor whose ligands, CCL27 and CCL28, are expressed on melanocytes, plasma cells and, of interest for protection against a skin ectoparasite, skin-homing T cells. 22 Interestingly, previous studies have found a greater accumulation of CD3 + T lymphocytes and WC1 + T lymphocytes in the infested skin of genetically tickresistant bovines, [23][24][25] which also recognize many more tick salivary proteins with IgG antibodies than the tick-susceptible breed   employed in the present study. 26 These results could signify that by neutralizing components of tick saliva with antibodies induced by the multicomponent vaccine, a tissue repair response is able to proceed and hinder tick attachment and blood feeding. The GSEA, in turn, was enriched for modules reflecting antigen presentation and phosphatidylinositol signaling of immune cells, among others. These modules were also reported in other studies to be enriched by vaccination with antigens from Plasmodium falciparum 27 or by intradermal challenges with sporozoites 28 of this parasite. In fact, there was great concordance between upregulated and downregulated transcriptional modules between this latter study and the present study. Phosphatidylinositol signaling is also significantly associated with neutralizing antibody titers after measles vaccination. 29 Finally, another blood transcriptional module enriched in the GSEA of the present study, Mismatch Repair (I) (M22.0 in Fig. 2), was also found to be correlated with serum neutralizing antibody titers in sheep vaccinated with inactivated foot and mouth disease virus. 30 Of Fig. 3 Metacore Pathway obtained with DEGs in blood from vaccinated calves. The pathway map "Role of IL-8 angiogenesis" obtained by functional analysis of differentially expressed genes in vaccinated calves: Thermometers 1 and 2 next to gene/name pictograms depict the genes present in AV and CHV comparisons, respectively, which are highlighted with opened dashed circles. The red thermometer indicates upregulated genes and the blue thermometer indicates downregulated genes; the thermometer level represents the intensity of the log 2 foldchange. Image generated using Metacore.
interest for the present study, this pathway is one of the mutagenic DNA damage responses that participate in somatic hypermutation to increase antibody affinity. 31 In conclusion, all these genes and pathways should play important roles in controlling tick infestation in vaccinated animals with the participation of molecules involved in skin immune responses. Finally, the bovine transcriptome data based on a Systems Vaccinology approach presented herein, will provide findings and insights to help the rational vaccine design for development of protective vaccines for ticks that rely on generating neutralizing antibodies.

RNA isolation
Whole peripheral blood samples were stabilised in PAXgene Blood RNA tubes, followed by isolation and purification using the PAXgene Blood RNA Kit (PreAnalytix) according to the manufacturer's instructions. Briefly, the blood samples were defrosted and centrifuged, and subsequently, the supernatants were removed. The pellets containing whole-blood cell lysate were washed and resuspended in the buffer provided in the kit, followed by the binding, DNase I treatment, washing and eluting of total RNA using the column systems and buffers provided in the kit. Total RNA samples were quantified using a NanoDrop 1000 spectrophotometer, and the RNA quality was analysed using a lab-on-chip Agilent 2100 Bioanalyzer. The RNA samples were precipitated and submitted to a sequencing service.

RNA sequencing
Sequencing data were generated on an Illumina platform at the High-Throughput Sequencing and Genotyping Unit of the Roy J. Carver Biotechnology Center, University of Illinois, Urbana-Champaign, Illinois, USA. All procedures were performed according to the manufacturer's instructions. The libraries were prepared using a TruSeq RNA Sample Prep kit and sequenced using a TruSeq SBS kit v3-HS on an Illumina HiSeq2000 system. Briefly, after sample quality control analysis, the poly (A) + mRNA samples were purified from total RNA using oligo(dT) magnetic beads, followed by fragmentation, cDNA conversion using random primers, end repair, 3′-end adenylation and adaptor (barcode) ligation. The 24 libraries were combined into six pools, quantified through qPCR, and sequenced in four lanes for 100 cycles. Single-end reads of 100 nucleotides were analysed using Casava 1.8 software (pipeline 1.9) and the de-multiplexed compressed fastq files downloaded from the server of the sequencing service unit.

RNA-seq data analysis
The quality control of the fastq files containing the raw sequence reads was performed using FastQC (Babraham Institute; available at: http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). No filtering step was performed with the fastq files, as the quality of the results from FastQC reports showed a high level of average quality reads. Next, the reads were aligned to the bovine reference genome (Bos taurus genome assembly UMD3.1 downloaded from Ensembl) using TopHat2 mapper, version 2.0.9. 32 The quantification of mapped reads was performed using HTSeq version 0.5.4, 33 whose read count outputs were used as inputs for differential expression analysis calculated with edgeR package version 3.2.4 34 using the generalised linear model (GLM) likelihood ratio test. A threshold of FDR < 0.05 were applied to obtain the differentially expressed genes in all comparisons. Because RNA-seq measures absolute numbers of transcripts and because qRT-PCR correlates poorly with those genes presenting either low or high levels of expression in transcriptomes, 35 we proceeded to validate the data biologically by seeking functional correlations with the bioinformatic strategies described in the next section, with the relevant responses measured in our previous studies, in particular the study that generated the samples employed herein, as well as with relevant responses measured in studies by other investigators; these correlations will be presented in the Results and Discussion sections. The RPKM (reads per kilobase per million) was calculated to normalise the read counts according to gene length (sum of exons for a transcript) based on edgeR user's guide instructions. The gene length was obtained using an inhouse perl script, using information from the genomic features of Ensembl IDs. For genes with more than one transcript, the larger transcript was considered. Hierarchical clustering (HCL) to examine the gene expression pattern was performed with the RPKM values using a standard statistical algorithm previously described. 36 HCL employed Euclidean distance for metric calculations, and the average linkage method was used to indicate the cluster-to-cluster distance in the hierarchical trees, which were displayed as heatmaps. HCL analyses were performed using MeV (Multi Experiment Viewer) software. 15 For functional analyses, the differential expression profiles obtained from all comparisons were functionally enriched using MetaCore™ from Clarivate Analytics, a web-based software for functional analyses based on high-quality, manually curated databases. Pathway and network analyses and Gene Ontology (GO) processes were retrieved using an integrated workflow available in MetaCore™. Because the databases of this software only contain information from humans and rodents, the bovine gene IDs were converted into orthologous human gene IDs using the BioMart database 13 through the Ensembl website. Gene set enrichment analysis (GSEA) was performed using Blood Transcription Modules (BTMs) as described elsewhere. 14 Log CPM values were used for comparisons to baseline gene expression. To determine significant associations of BTM with tick challenge phenotypes after adjuvant injection or vaccination, we used baseline-normalized gene expression, as previously described. 37 Parameters included signal to noise ranking metric and weighted enrichment statistic with 1000 permutations.