Candidate genes for tick resistance in cattle: a systematic review combining post-GWAS analyses with sequencing data

ABSTRACT

Rhipicephalus microplus causes huge losses in cattle. Host genetic background greatly affects the immune efficiency in resistance or susceptibility to tick infestation, which is one of the many factors that play a role on that trait. We performed a systematic review of genome-wide association studies (GWAS) for tick resistance in cattle resulting in 1353 candidate genes for post-GWAS analyses. From those, genes showing possible structural variants from the bovine genome were classified by the Variant Effect Predictor from Ensembl. Ninety-two candidate genes showed potential structural variants in 5′ UTR and coding region and were used for functional annotation. Enriched biological processes (e.g. regulation of eosinophil chemotaxis, RIG-I signalling pathway and monocyte differentiation) and candidate genes (e.g. DAPK2, PUM1, ACIN1, INPP5D) linked with immune system function were identified and thus associated with tick resistance. Besides, gene-transcription factors (TFs) networks were obtained from TFs associated with immune system (FOXO3, PPARG, STAT3, NFKB1, GATA3 and ARNT) and the candidate genes associated with tick resistance in cattle highlighted (e.g. OR4L1, PNP, LRRIQ1, GIMAP8, MYO6, MEP1A and LRFN2). Thus, promising candidate genes with a possible functional role for tick resistance in cattle are presented for further in vitro and/or in vivo analyses.

KEYWORDS:

• Ectoparasite
• genome
• immune system
• gene networks
• dairy
• beef

1. Introduction

Cattle production is greatly impacted by the ectoparasite Rhipicephalus microplus infestation, which deeply impacts the animal health resulting in losses of meat, milk and leather production (Gueretz et al. 2020). Infestation control is largely based on the use of chemicals, which may include compounds with high toxicity and environmental persistence (Bai and Ogbourne 2016), and consequently harmful to human, animal and environment health. In addition, the recurrent and massive use of these compounds cause ticks to acquire resistance against the active principles used (Garcia et al. 2019), hindering the viability of their control. In this way, alternative strategies aiming to mitigate the losses caused by ectoparasites infestation must be offered for both, dairy and beef production.

It is known that Bos indicus breeds are genetically more resistant to ticks than Bos taurus breeds, probably due to naturally selected genes throughout their evolutionary process (Ayres et al. 2015). This indicates that resistance to tick infestation may be a potential factor to be studied in animal breeding programmes (Cardoso et al. 2014). Hence, genome-wide association studies (GWAS) identifying genomic regions associated with tick resistance (Neto et al. 2010; Porto Neto et al. 2010; Mapholi et al. 2016; Sollero et al. 2017; Otto et al. 2018) have been carried out. However, analyses of post-GWAS covering biological knowledge of the genes underlying this trait, as well as the possible potentially associated variants, are still scarce.

Post-GWAS analyses using gene networks have been carried out in the main livestock species including bovine (Verardo et al. 2016; Otto et al. 2018; Littiere et al. 2020). However, the use of these tools combined with sequencing data may be a promising way to identify the most candidate genes. After the advent of the next-generation sequencing technology, genome resequencing has increased the chances of finding variants in genes associated with important phenotypes (Schneeberger and Weigel 2011), and the structural variants are a promising class. Structural variants are usually determined as a DNA region of approximately 50 base pairs (bp) (Conrad et al. 2010) to 1 kb or more (Freeman et al. 2006). They can be classified into inversions and translocations, segmental duplications (SDs) and copy number variants (CNVs) which are unbalanced genomic events (Feuk et al. 2006).

A previous study described that CNVs might strongly influence gene expression through deletions/duplications changes in the position and function of cis-regulatory elements (Spielmann et al. 2018). For example, deletions variants located in genes involved with the human innate immune system affected the clinical infection caused by the Andean orthohantavirus (Ribeiro et al. 2019). Since the immune system is one of the important mechanisms that confers bovine resistance or susceptibility to tick infestation (Turner et al. 2022), it is possible that candidate genes underlying the immune response may show structural variants. Thus, in this study, we identified candidate genes for tick resistance in cattle that may present structural variants to perform functional annotations via gene networks (biological processes and transcription factors networks) and highlight the most candidate genes for further studies.

2. Materials and methods

2.1. Systematic review and identification of candidate genes

A systematic review was carried out in May 2020 to identify candidate genes related to tick resistance in cattle. We searched for articles using the Google Scholar search engine (https://scholar.google.com.br) in terms to maximize the search by not discriminating specifics databases, and thus assuming to consider all available articles database. The search consisted of combinations of keywords (in English) with the following criteria terms: (A) term related to the evaluated trait (‘tick resistance’, ‘tick’, ‘inflammatory response’); (B) type of association test (‘GWAS’, ‘genome wide association’) and (C) species/breed (‘Holstein’, ‘Gir’, ‘Girolando’, ‘cattle’). It is important to note that the use of breed term was first intended to search for a Brazilian national breed without compromising the final results, in terms that we used the ‘cattle’ keyword, and then focused on studies of all breeds. Two independent reviewers (1 and 2) searched for peer-reviewed papers using the queries described. Discrepancies in the judgment were resolved by consensus among the reviewers. The first step in selecting articles was to check whether the article contained (in the abstract, title or keywords) the keywords used in each combination. The articles that met these criteria were selected, and if duplicated, articles were removed. In addition, we confirmed that: (1) the article was published from a peer-reviewed journal; (2) GWAS was performed for traits related to tick resistance; (3) cattle were used as a model; (4) the consistent position of the gene and the methodology used were provided and (5) full text available. In addition, articles that did not provide enough information about the associated markers or windows (e.g. marker names and / or ‘rs’ and / or complete genomic coordinates) were also excluded from this step.

Each article was defined as a group for the identification of candidate genes. The positions of significant markers in each group were updated according to the ARS-UCD1.2 genome assembly, and the search for candidate genes carried out following the flanking region defined in each study. The search for the genes was carried out through the NCBI database (https://www.ncbi.nlm.nih.gov/gene/?term = bovine). A distance of 22 kb upstream and 22 kb downstream from each significant marker was considered for studies that did not show the flanking region for candidate genes identification, since it is the half of the average distance between markers in the SNPChip used in the studies. By checking this distance, we identified the presence of any gene that could be related with a QTL or a significant marker, as proposed by Verardo et al. (2016).

2.2. Identification of structural variants associated with genes identified in systematic review

The public Ensembl database was used (‘ftp://ftp.ensembl.org/pub/release-101/variation/vcf/bos_taurus/’, accessed on November 18, 2020) to search bovine structural variants. Aiming to select target variants, we first classified the variants according to their potential function using the web interface of Ensembl Variant Effect Predictor (VEP) tool (McLaren et al. 2016), prioritizing structural variants located in the coding and non-coding regions of the bovine genome. Only variants associated with known genes identified in systematic review were selected in this step.

2.3. Functional analysis of candidate genes for tick resistance

From the sets of candidate genes that may present structural variants, we selected only the genes showing variants in the 5′ untranslated region (5′ UTR) and coding regions according to VEP classification. Thus, only genes containing such target variants were used for the construction of gene networks.

We first obtained a biological process network. From each group, we used the ClueGO application of Cytoscape (Bindea et al. 2009) to highlight biological processes associated with tick resistance. The analyses of ClueGO were based on a hypergeometric test with Bonferroni correction.

In addition, we also built a gene-transcription factor (gene-TF) network. Thus, the promoter sequences of the candidate genes of each group were analysed using the TFM-Explorer programme (http://bioinfo.lifl.fr/TFM/TFME) to search for transcription factors binding sites (TFBS). This programme uses weighting matrices from the JASPAR vertebrate database (Sandelin et al. 2004) to detect potential TFBS, calculating a scoring function as described in Touzet and Varré (2007). From each group, the promoter region was defined by selecting the sequences at 3000 bp upstream and 300 bp downstream (FASTA format) from the transcription start site, based on the bovine genome assembly ARS-UCD1.2. This data was used as an input to TFM-Explorer. The list of transcription factors (TFs) obtained was analysed in Cytoscape 3.7.2 (Shannon et al. 2003) via BinGO plug-in (Maere et al. 2005) to determine which terms of gene ontology were significantly over-represented, assuming the standard analysis and multiple correction test.

Based on the biological processes associated with tick resistance, besides the literature review, the main TFs (key TFs) in each group related with the studied trait were selected. By this way, key TFs with known roles in biological processes related with tick resistance trait such as immune system were used to build gene-TFs networks. This approach has already been used in other species, including bovine (Verardo et al. 2016; Otto et al. 2018; Littiere et al. 2020), being that this judicious review of each TF is taken to avoid false associations between TF and the studied trait.

Through the Network Analyzer tool in Cytoscape v.3.7.2, promising candidate genes were highlighted according to the number of TFBS and, consequently, the number of connections/lines (edges) in each node (gene and TF), and the most connected genes in the gene-TF network were determined. Genes with more TFBS for the most representative key TFs were highlighted in the gene-TF network. Thus, the networks of biological processes and gene-TF highlighted candidate genes associated with tick resistance in cattle. A summary of whole analysis is shown in Figure 1.

Figure 1. Flowchart summarizing the steps towards identification of most candidate genes associated with tick resistance in cattle.

3. Results

3.1. Systematic review

The systematic review process is summarized in Figure 2. The searches performed with 24 keyword combinations returned 43 articles, 23 found by reviewer 1 and 20 found by reviewer 2 (Supplementary Material – Table S1). Of those, 36 were duplicated between combinations of keywords and removed. The remaining seven articles proceeded to the full reading stage and three of them (Kim et al. 2015 ; Neto et al. 2011 ; Turner et al., 2010) were removed by not meeting the selection criteria (e.g. not be a tick resistance GWAS study or did not provide proper information about marker names and / or ‘rs’ and / or complete genomic coordinates). Thus, four articles were selected and defined as groups (from 1 to 4) for functional analysis, as shown in Table 1. Information about genotyping and quality control of each selected paper are presented in Supplementary Material – Table S2. Based on the significant SNPs/windows and their respective positions (bp) on the chromosomes described in each selected article, a total of 1353 candidate genes were identified (Supplementary Material – Table S3).

Figure 2. Flowchart showing the pipeline used to select articles based on specific criteria for systematic review.

Table 1. Information of country, breed, biological samples, and traits selected from articles related to tick resistance in cattle.

Article (Published year)	Group	Trait	Biological sample	Breed	Country
Mapholi et al. (2016)	1	Counting of ticks based on body location in the animal and tick species	Blood and hair	Nguni	South Africa
Sollero et al. (2017)	2	Counting of ticks	Blood, hair and semen	Hereford and Braford	Brazil
Mota et al. (2018)	3	Counting of ticks	Blood and semen	Hereford and Braford	Brazil
Otto et al. (2018)	4	Counting of ticks	Blood	Gir × Holstein F2	Brazil

3.2. Structural variants associated with genes

The Ensembl public database showed 462,845 structural variants in cattle. From them, 448,045 were classified according to their effect, using the VEP, which were associated with 14,593 genes. These variants were in gene regions such as 5′UTR, coding regions, introns and intergenic. The genes identified with these variants were used in the next step of the functional analyses.

3.3. Functional analysis of candidate genes for tick resistance

From the candidate genes identified on systematic review (1353) and those with possible structural variants identified by the VEP (14,593), a merged list was obtained with 626 remaining genes (Supplementary Material – Tables S4–S7). From them, considering the four groups, a total of 92 candidate genes showed possible target structural variants in the coding region and 5′UTR (Supplementary Material – Tables S8–S11). The target variants were classified as ‘modifier’ and ‘high’, showing impacts on the protein according to the VEP classification. In addition, we observed that the target variants located in candidate genes for tick resistance were defined as duplication, tandem duplication, gain of copy number, loss of copy number and deletion.

Each group of candidate genes for tick resistance displaying possible target structural variants were used to obtain genes-biological processes networks. These networks highlighted different processes associated with tick resistance in cattle (Figure 3). Some of highlighted processes were the regulation of eosinophil chemotaxis associated with the DAPK2 gene (identified in group 2), positive regulation of RIG-I signalling pathway associated with PUM1 gene (identified in group 2) and regulation of monocyte differentiation that was in common between ACIN1 and INPP5D genes (identified in group 1).

Figure 3. Functional networks between genes and biological processes associated with tick resistance in cattle. Main interaction networks between biological processes and genes (nodes) are shown in zoom. Blue nodes are associated with group 2 and red nodes are associated with group 1. Grey node represents process shared between the genes of the two groups. The size of the node characterizes the enrichment of the process according to ClueGO Cytoscape plug-in (Bindea et al. 2009). The terms that are most enriched by subnet are in bold.

According to TFM-Explorer programme, 65 TF were related to the candidate genes (Supplementary Material – Tables S12–S15). Based on biological processes (Supplementary Material – Tables S16–S19) and literature evidence, six key TF associated with tick resistance in cattle were selected (Table 2). These key TF were used to generate gene-TF networks for each group (Supplementary Material – Figures S1–S4). Based on that, a merged network was built (Figure 4) enabling to highlight promising candidate genes. Thus, from the 92 candidate genes with possible target variants, promising candidate genes for tick resistance in cattle were observed (e.g. CACNA1A, LRRIQ1, NKD2, GIMAP8, WIPF3, LRRC4C, LOC112449525, LOC536229, MYO6, SPHKAP, OR4K2, OR4F1, NOR4F1, NOR4F1, NOR469, SORCS3, PNP, MEP1A, LRFNE). A summary of the main highlighted candidate genes number in each group after gene network analysis is presented as a workflow (Supplementary Material – Figure S5)

Figure 4. Gene-transcription factors network for tick resistance in cattle. Transcription factors (yellow diamond nodes) associated with candidate genes (circle nodes in red: group 1; blue: group 2; light pink: group 3; green: group 4) for tick resistance in cattle. The yellow diamond node with the purple border is associated with three groups (1, 2 and 3) while the others are only associated with the group of the border colour. Dark pink square nodes represent the biological processes associated with transcription factors. The size of the nodes represents the enrichment as a function of transcription factors binding sites number.

Table 2. Key transcription factors (TF) for tick resistance in cattle. Most representatives TF associated with genes identified for each group, based on their biological processes and evidence in the literature.

TF	Group	Biological Process	Literature evidence
FOXO3	1	Regulation of erythrocyte differentiation	Promote differentiation a function of T cells (Qi et al. 2020).
PPARG	1	Myeloid leukocyte differentiation and monocyte differentiation	Immune cells differentiation (Christofides et al. 2020).
STAT3	1,2 and 3	JAK-STAT cascade	Development and function of innate and adaptative immunity (Yan et al. 2018).
NFKB1	2	Regulation of interleukin-12 production and negative regulation of cytokine biosynthetic process	Development of innate immunity (Dev et al. 2010).
GATA3	4	Regulation of interleukin-4 production	Development and function of adaptative immunity cells (Wan 2014).
ARNT	4	Regulation of cytokine production	Control of adaptative immune response (Quintana and Sherr 2013).

Note: *The literature evidence cited is only a sample of the vast available literature.

4. Discussion

4.1. Systematic review

According to the systematic literature review results, we observed that GWAS for tick resistance in cattle is little explored still, in terms that only four scientific articles were identified meeting the established criteria in our survey. Among them, we observed three studies in Brazil (Hereford and Braford breeds, and Gir × Holstein F2 animals) and one in South Africa with Nguni breed. The locally adapted Nguni breed is classified as Bos taurus africanus (Nyamushamba et al. 2017) and known by a good fertility and resistance to diseases (Mapiye et al. 2020). When compared to exotic temperate breeds, it has lower levels of nematode infestation (Ndlovu et al. 2009), tick count (Spickett et al. 1989; Muchenje et al. 2008; Marufu et al. 2011) and seroprevalence for Anaplasma marginale and Babesia bigemina (Marufu et al. 2010). In Brazil, the Gir × Holstein F2 cattle are mostly used in milk production, as they adapt to different production systems and are quite efficient (Ribeiro et al. 2017). According to Otto et al. (2018), candidate genes for tick resistance in these crossbred cattle were identified and, according to the origin of allele, resistant cattle showed two alleles of Gir breed, while the susceptible ones possess alleles from Holstein breed.

In addition, Hereford and Braford cattle are predominant beef breeds in the southern region of Brazil (Piccoli et al. 2020). Studies using these breeds have shown that tick resistance is a trait with a moderate heritability estimate (0.13–0.19) and shows a positive genetic correlation with production traits (Cardoso et al. 2015; Biegelmeyer et al. 2017). Knowing that ticks cause enormous damage to cattle, especially in tropical and subtropical regions, where the climate is favourable for their development, such as Brazil and South Africa, studies of the genetic architecture of tick resistance trait is extremely important.

Thus, from the systematic review of this study, it was possible to shed light on the GWAS for tick resistance in cattle, from different regions and breeds of beef and dairy cattle. Several genomic regions associated with tick counts in cattle have been identified (Mapholi et al. 2016; Sollero et al. 2017; Mota et al. 2018; Otto et al. 2018). Thus to better understand the genetic architecture of this complex trait, a post-GWAS analysis of candidate genes combined with identification of possible target structural variants was performed for the first time in the literature.

4.2. Functional analysis of candidate genes for tick resistance with possible structural variants

From 1353 candidate genes obtained via literature review, 626 showed possible variants in several gene regions. Among them, 92 genes showed target structural variants in the 5′ UTR and coding regions. French and Edwards (2020) demonstrated that genetic variants in the 5′ UTR region have a potential effect on affecting the expression or function of genes, as they can alter regulatory elements that influence their interaction with 5′ UTR, and consequently impact on post-transcriptional and translation processes. In addition, variants located in coding regions hold the potential to alter the protein sequence, which may lead to loss and alteration of the protein’s function (Zhao et al. 2019). Since genes identified from GWAS can contribute to the identification of inter-individual variant differences in phenotypes (Wang et al. 2019), it is possible that target structural variants may be located in the 92 candidate genes and can have the potential to influence the expression of tick resistance in cattle. Thus, we performed functional analyses of these genes aiming to highlight the most candidate genes associate with the studied trait.

From the biological processes networks, it was possible to observe that the ACIN1 and INPP5D genes from group 1 and DAPK2 and PUM1 genes from group 2, were enriched and associated with biological processes with probable involvement in cattle tick resistance. The ACIN1 gene encodes the Apoptotic Chromatin Condensation Inducer 1 protein, which has been shown to be important in the nuclear matrix (Mendonça et al. 2018) and apoptosis processes, contributing to several functions in the immune system, mainly in the adaptive response (Nagata and Tanaka 2017). For example, in humans, cells that are potentially cancerous and contaminated by viruses are eliminated from the body through apoptosis (Bedoui et al. 2020). The INPP5D gene belonging to the inositol Polyphosphate-5-Phosphatase family that is expressed in haematopoietic cells (Ware et al. 1996) and plays a role in several pathways, such as microglial activation, neuroinflammation, and immune response (Malik et al. 2015; López González et al. 2016; Yoshino et al. 2017). Moreover, Kong et al. (2019) integrated proteomics and microRNA from the plasma of Jersey cattle, in response to high altitude hypoxia (HAH), inferred that miR-155-5p regulates the immune response of B cells by targeting INPP5D to HAH stress adaptation. In general, both ACIN1 and INPP5D shared important biological processes of the immune system in our analysis, such as regulation of monocyte differentiation. Monocytes perform different functions in the course of infection, inflammation, injury and healing (Okabe and Medzhitov 2016). When these undergo differentiation, their effector functions are potentiated, being able to induce long-term adaptive immune responses (Witte et al. 2018). Thus, ACIN1 located in BTA10 and INPP5D in BTA 3 are potential candidate genes for tick resistance in cattle by playing important roles in the adaptive immune system.

Considering the selected genes from group 2, we also observed biological processes that influence the immune system, such as regulation of eosinophil chemotaxis associated with the DAPK2 gene and positive regulation of RIG-I signalling pathway associated with the PUM1 gene. DAPK2 belongs to the death-associated protein kinase family of serine–threonine kinase, encoding the protein kinase 2 associated with cell death (Humbert et al. 2017). It shows functions in the regulation of apoptosis, autophagy and inflammation (Geering 2015), and is mainly expressed in neutrophil and eosinophil cells in humans (Humbert et al. 2017). It also plays an important role in granulocytes differentiation (Rizzi et al. 2007; Humbert et al. 2014). Granulocyte chemotaxis is one of the essential elements of innate immune response to injury or infection, making the response more efficient, in terms that this process recruits and attracts the necessary cells to the area of ⁣inflammation (Geering et al. 2014). Studies involving the immune responses in the skin of cattle infested with ticks have indicated a strong hypersensitivity reaction with increased resistance of cattle to ticks, including the infiltration of eosinophils and the concentration of histamine in the tick fixation site (Schleger et al. 1976; Kemp and Bourne 1980; Schleger et al. 1981). In a histological study of the tick bite location, a high infiltration of basophils, mast cells and eosinophils were observed in the skin of resistant Bos indicus cattle, while in susceptible Bos taurus animals, a greater recruitment of neutrophils was observed (Tabor et al. 2017).

Moreover, the PUM1 gene encodes the pumilio RNA binding family member 1 protein, and acts as an important regulator of key genes of the innate immune system (Liu et al. 2017). In addition, the biological process positive regulation of the RIG-I signalling pathway found in our study is essential to the innate immune system to identify cells that have been contaminated by viruses (Kell et al. 2015 ). Interestingly, the flavivirus is one pathogen that can be transmitted by R. microplus tick (Souza et al. 2018), and RIG-I signalling pathway seems to play an important role in the immune system regarding the recognition and elimination of these pathogens by preventing their invasion (Guo et al. 2018).

In this study, we also identified TF related to candidate genes for tick resistance. Among the 65 identified TF, six were highlighted as key TF associated with the immune system (FOXO3, PPARG, STAT3, NFKB1, GATA3 and ARNT) which were used to build a gene-TF network. Among them, STAT3 was identified in common to groups 1, 2 and 3 and, consequently, in three breeds: Nguni (group 1), Hereford and Braford (groups 2 and 3). The STAT3 is a downstream effector of several cytokines, an important regulator of immune responses via JAK/STAT signalling (Yan et al. 2018). From group 1, the FOXO3 was found associated with the regulation of erythrocyte differentiation and PPARG was associated with the differentiation of myeloid leukocytes and monocytes. In a study with mice, modulation of FOXO3 activity showed that it could increase or prevent antibody-mediated immune response, playing a role in regulating the differentiation and function of T helper cells (Qi et al. 2020). In addition, PPARG plays a role in differentiating several immune system cells such as macrophages, monocytes, dendritic cells, T and B lymphocytes and platelets (Croasdell et al. 2015; Christofides et al. 2020). Piper et al. (2009) suggested that resistant cattle (Bos indicus) show a greater tendency to potentiate a consistent T cell-mediated response against R. microplus.

Another well-enriched TF in the gene-TF network was the NFKB1, which was associated with genes of the group 2. This TF belongs to nuclear factor-κB (NF-κB) family that regulates a diversity of genes with functions in the development of immune system, inflammation, innate and adaptive immune responses (Mulero et al. 2019). Moré et al. (2019), analyzing differentially expressed genes in Braford cattle, demonstrated that in susceptible cattle, members of the NF-kB family were the mainly expressed. This gene family has also been observed to be expressed in tick fixation sites in Holstein cattle (Piper et al. 2008). In our study, NFKB1 was associated with biological processes such as regulation of the production of interleukin-12 and negative regulation of cytokine biosynthetic process, what could be a potential TF associated with candidate genes for tick resistance in cattle.

In the group 4, GATA3 was associated with the regulation of interleukin-4 (IL4) production process. This TF is known to act as a regulator of the innate and adaptive immune system through the differentiation of memory cells such as T CD4 + and T CD8 +, Th2 and innate lymphocytes (Wan 2014). Rodriguez-Valle et al. (2013) demonstrated that susceptible Holstein cattle exposed to R. microplus showed lower expression of cytokine Th2 related to IL4 if compared with the resistant Brahman cattle. In addition, IL4 is a key regulator of humoral and adaptive immune responses (Tabor et al. 2017). Also observed in group 4, ARNT is expressed by several cells of the immune system and it is known to integrate actions of the environment and metabolism in the immune response (Gutiérrez-Vázquez and Quintana 2018). In mice, increased expression of ARNT mainly in myeloid cells, in addition to contributing to the immune system, also indicates the possibility of an excellent therapeutic strategy to improve wound healing (Scott et al. 2014). The wound healing process is regulated by several growth factors and cytokines released at the injury site (Kurahashi and Fujii 2015). In resistant x susceptible Braford cattle, after the tick infestation, the proteolysis processes in the degradation of the connective tissue and remodelling of the extracellular matrix, inflammation through cytokine signalling and the platelet-endothelium-leukocyte interactions were enriched only in resistant cattle (Moré et al. 2019). Thus, it is suggested that GATA3 and ARNT through the regulation of cytokines production may be influencing the healing process, in response to injury caused by ticks in cattle, thus contributing to their resistance to ticks.

Additionally, we also observed genes enriched by binding sites for the highlighted TFs (e.g. OR4L1, PNP, LRRIQ1, GIMAP8, MYO6, MEP1A and LRFN2). From them, OR4L1 and PNP were identified on group 1. The OR4L1 gene belongs to the family of olfactory receptors (Malnic et al. 2004). Lee et al. (2018) identified this gene among the 30 main genes encoding proteins with differential methylation due to Rh2 Ginsenoside, which is known to increase the activity of cells in the immune system. In addition, the purine nucleoside phosphorylase (PNP) gene, in BTA10, is well recognized in several studies about its function in the immune system (Ghodke-Puranik et al. 2017; Albar et al. 2020; Tecle et al. 2021). In humans, variants found in this gene that cause the reduction or absence of its activity, were found in autoimmune diseases, gradual immune T-cell insufficiency, B-cell dysfunction, and linked to high susceptibility to infections (Albar et al. 2020). Tecle et al. (2021) indicated that purine homeostasis in the host regulates resistance to contamination from intracellular and extracellular pathogens. Thus, OR4L1 and PNP are highlighted here as promising candidate gene for tick resistance in cattle.

In the group 2, the GIMAP8 and LRRIQ1 genes were highlighted. GIMAP8 is one of the members of GTPase family of immunity-associated proteins (GIMAPs), strongly expressed in the terminal stages of T and B cell development (Berg et al. 2021). Robbertse et al. (2018) demonstrated that in bovine tick fixation areas, B lymphocytes are important mediators of immune responses, due to their considerable influx and significant increase in CD3 + T lymphocytes, observed in tick-resistant breeds. The LRRIQ1 gene encodes the leucine-rich repeat and IQ motif containing 1 protein, which indicates evidence that this gene acts on the immune system, recognizing pathogen effector proteins and triggering localized cell death, to limit infection in animal and plant cells (Horsefield et al. 2019). Moreover, proteins containing leucine-rich repeats are also known to function as intracellular immune receptors, protecting hosts against pathogen invasion (Mermigka et al. 2020). In this way, we suggest that GIMAP8, in BTA4, and LRRIQ1, in BTA5, may be promising candidate genes since they play an important role in cattle tick resistance.

The myosin VI gene (MYO6) was identified in the group 3 and is involved in a wide range of cellular processes such as endocytosis, exocytosis, autophagy and regulation of actin dynamics (Jonge et al. 2019). In inflammation situations, it plays a role in the epithelial barrier through endosome/lysosome fusion in epithelial cells (Liao et al. 2013). Mota et al. (2018) also highlighted MYO6 as a potential candidate for tick resistance in cattle, due to its possible action on the epithelial barrier, which corroborates our findings in gene-TF network results. Thus, we suggest that the MYO6 gene in BTA9 is a potential candidate gene for resistance to ticks in cattle by playing a role in the epithelial barrier of these animals.

In the group 4, the meprin A subunit alpha (MEP1A) gene is expressed in several mammalian cells and organs, such as kidney, intestine, lung, bladder, skin, cancer and immune cells, proteolytically processing various inflammatory mediators such as cytokines and chemokines (Herzog et al. 2019). In Brangus cattle, MEP1A was identified in the bovine leukocyte antigen (BoLA) region on BTA23, which is responsible for encoding genes associated with the adaptive immune system and may be related to the adaptation of these animals to subtropical environments (Goszczynski et al. 2018). In addition, in the same group 4, the LRFN2 gene was initially described to be expressed in the brain (Brouwer et al. 2019), but studies have identified its expression in non-neural tissues and associated with the risk of lung cancer (Jin et al. 2012; Bhat et al. 2020). In addition, a member of his family, the LRFN4 gene, has been found to function in immune cells in monocyte/macrophage migration (Konakahara et al. 2011). Thus, we suggest that the MEP1A and LRFN2 genes in BTA23 are promising candidate genes for tick resistance in cattle.

It is important to note that even though we have focused on immune system-related processes in this study, there are other factors such as skin thickness and coat colours that may also play a role in tick resistance. Foster et al. (2007) suggested that the skin thickness appears to play an important role on host resistance to ticks. Moreover, Gasparin et al. (2007) analyzing contemporary group, revealed that lighter coloured animals are more resistant than dark coloured ones. However, in terms to address the physiological role in host tick resistance, we focused on immune system processes.

In cattle, the immune system is one of the important mechanisms that confer resistance or susceptibility to tick infestation. Ribeiro et al. (2019) observed that deletions variants found in genes involved with the human innate immune system affected the clinical infection caused by the Andean orthohantavirus. The variants used in this present study were deletion, duplication, tandem duplication, gain in copy number and loss in copy number, besides the possibility to be located in genes associated with the immune system, we may infer that these genes are strong candidates to be associated with tick resistance. Based on that, we present a robust post-GWAS analyses approach via the systematic review followed by functional analyses of candidate genes with possible variants for important economical traits in livestock.

5. Conclusions

Candidate genes associated with tick resistance in cattle presenting possible target variants were identified. The gene networks allowed the identification of enriched biological processes related to the immune system and promising candidate genes associated with tick resistance in cattle. Furthermore, the gene-transcription factors networks allowed the identification of putative main TF and candidate genes for tick resistance in cattle. Thus, promising candidate genes with a possible functional role for tick resistance in cattle are presented for further in vitro and/or in vivo analyses.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

We would like to thank the support of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES) – Financing Code 001, and Fundação de Amparo à Pesquisa do Estado de Minas Gerais – FAPEMIG (process APQ-01834-18).