Introduction {#Sec1}
============

Melanoma is a type of cancer arising from the malignant transformation of melanocytes, pigment producing-cells found predominantly in the basal layer of the epidermis and eyes.  Cutaneous melanoma is the most aggressive and treatment-resistant form of skin cancer responsible for the vast majority of skin cancer-related deaths in the Caucasian population^[@CR1]^. The global incidence of melanoma continues to increase at an alarming rate, despite decades of public prevention programs in many countries.  Around 232,000 new cases of skin cancer were recorded worldwide in 2012, accounting for 1.6% of all new cases of cancer back then, while over 300,000 new cases of melanoma were diagnosed worldwide in 2018, according to the World Cancer Research Foundation^[@CR2],[@CR3]^. In Brazil, 1,547 deaths were recorded in 2013 due to melanoma, with around 5,690 new cases reported back then, while around 6,260 new cases were expected due in 2018, according to the National Cancer Institute (INCA)^[@CR4]^.

Cutaneous melanoma usually affects a higher proportion of patients, in the age range 40--60 years.  They can be treated by surgical excision when detected in the early stage (0, I, II and resectable III), however, in the later stages (unresectable III, IV and recurrent melanoma) the treatment options are chemotherapy, target therapy (BRAF/MEK pathway inhibitors), immunotherapy (checkpoint blockade CTLA-4 receptor inhibition, PD-1 ↔ PD-L1 axis inhibition, and interferon-gamma immunotherapy), or a combination of them.  Death in most patients is caused by metastatic disease which may have evolved from the primary tumor.  Therefore, there is a need for new strategies in melanoma therapy because many patients do not respond to these treatments due to mutations in their genome, microbiota profile, response to neo-antigens, and resistance to drugs^[@CR5]^.

In horses, skin tumors are the most common among neoplasms.  Horse melanomas are encapsulated complex systems surrounded by a hard peel and represent 5 to 14% of equine skin neoplasms^[@CR6],[@CR7]^. These melanomas are usually benign tumors; however, they may have an unpredictable malignancy, progressing to malignant forms, and metastasize^[@CR6]^. Dark-haired horses feature higher melanoma malignancy as they get older compared to white-haired horses^[@CR7]^.

Malignant transformations are caused by driver mutations.  The acknowledged mutated genes in human melanoma are *BRAF* and *NRAS*, as well as the newly discovered *PPP6C*, *RAC1, SNX31, TACC1, STK19*, and *ARID2*^[@CR8]^. In horses, the only recognized driver mutation hitherto is caused by a 4.6 Kb duplication in intron 6 of Syntaxin 17 (*STX17*)^[@CR9],[@CR10].^ As is known, spontaneous horse melanoma tumors have a different natural history and many anatomical differences, compared to human melanoma tumors, but for some^[@CR11],[@CR12]^ it is considered a good translational model. 

In the present study, we treated horse melanoma tumors with Amblyomin-X, a Kunitz-type recombinant protein identified in a cDNA library from the tick *Amblyomma sculptum* salivary glands^[@CR13]^. This molecule has the ability to inhibit Factor Xa in the blood coagulation cascade and triggers apoptosis by activating its intrinsic pathway in tumor cells^[@CR14]--[@CR16]^. In previous studies, Amblyomin-X showed more avidity in the recognition of tumor cells, and rapid excretion in healthy murines^[@CR17]^. We have already shown that Amblyomin-X causes cell death via induction of both endoplasmic reticulum (ER) stress and proteasome inhibition (PI) in renal adenocarcinoma (murine RENCA), in melanoma (murine B16F10 and human SK-MEL-28)^[@CR18]^, as well as in pancreatic (human Mia-Paca-2) tumor cell lines^[@CR17],[@CR19],[@CR20]^.

Over the past ten years, preclinical studies have shown the potential of proteasome inhibition as an effective therapeutic strategy for treating different types of cancers^[@CR21]--[@CR24]^ including melanoma^[@CR25]--[@CR29]^. Bortezomib and Amblyomin-X are proteasome inhibitors (PIs) activating different innate immune pathways, the former via MAVS-NOXA^[@CR25]^, and the latter via RIG-I-MAVS, as suggested by our transcriptome analysis.  Trials with Bortezomib^[@CR30]--[@CR35]^ alone have shown insufficient clinical efficacy, and in combination with paclitaxel and carboplatin, limited clinical benefit and significant toxicity^[@CR30]^. Therefore, the effort to study proteasome inhibitors in different cancer therapies progresses^[@CR31],[@CR36]^.

Herein, we investigated the effect of Amblyomin-X intratumoral injection in horse melanomas, followed by interactome and transcriptome analyses conducted to understand genes modulation and respective enriched pathways.  The data suggest that Amblyomin-X modulates several pathways responsible for tumor cell death and its immune microenvironment.  Moreover, transcriptome analysis suggests potentially yet unknown DEGs and novel target pathways, mainly in relation to the responses of the immune system. 

Results {#Sec2}
=======

Equine melanoma regression upon treatment with Amblyomin-X {#Sec3}
----------------------------------------------------------

We injected Amblyomin-X intratumorally, at 1 mg/kg of the tumor mass, on each third day, during 28 days, and followed the tumor volume evolution and clinical animal conditions by over a period of 5 months.  We treated some selected tumors in the ventral tail of four animals by following three groups: treated (Amblyomin-X), and two control groups, vehicle (PBS) and untreated (see Materials and Methods).  Tumor dimensions were recorded before each injection and calculated the relative volume means related to the initial volume (Fig.  [1](#Fig1){ref-type="fig"}).  In all cases Amblyomin-X treatment promoted a reduction of at least 75% in tumor volume at the end of the treatment period (one month).  Tumor \#5 (from animal number 417) completely disappeared after 12 days under treatment, while tumors from another animal (number 795) continued to show volume reduction, even one month after the end of treatment.  Interestingly some animal control tumors (untreated and PBS-injected vehicle), that were in proximity to tumors treated with Amblyomin-X, also showed regressions in tumor size, which might suggest an indirect effect that could be related to the lymphatic system or recruitment of the immune system cells.  In two other animals, control tumors maintained or increased their initial volumes.  The evolution curves of each tumor per animal are shown in Fig.  [2](#Fig2){ref-type="fig"} (see also, Table [S1](#MOESM1){ref-type="media"}).  It is worth to mention, that each horse tumor is presented as an independent capsule.  We believe that long-term development and shrinkage are similar to all tail horse melanoma tumors, but experiments showed that each tumor responds differently, and we don't know the reason. Figure 1*In vivo* animal 797 tumor treatment and excision.  We can see on the left (**a**) animal 795 tail and chosen tumors: 1 and 8 as control, 2, 3 and 5 as vehicle, 4, 6, and 7 as treatment; (**b**) horse tail and tumor evolution on day 0, 15, 26 and 90 (excision); (**c**) veterinarian injecting Amblyomin-X solution or PBS in a tumor; and (**d**) T6, one of the excised tumors.  Tumors were chosen related to distance and size.  Three tumors, with a diameter of around 1.5 cm, were chosen for treatment.  About three farthest tumors were chosen as vehicle and other three tumors as control, whenever possible. Figure 2This plot shows the percentage volume curves related to the initial volume (0 h).  Every three days, Amblyomin-X was injected for about one month's time.  Animals 795 and 797 and their tumors were followed for at most 5 months.  Treated tumor volume curves were colored red, control tumor volumes green, and vehicle tumor volumes blue, for animals 417 in (**a**), 438 in (**b**), 795 in (**c**) and 797 in (**d**).  Each treated curve shows a decrease in tumor volumes.  Vehicle and control animals also showed a decrease in tumor volumes for animals 417 and 438.  Animals 795 and 797 showed either a constant or increased volume related to the initial volume.  Since each tumor has its own history and fate, the rate of decrease is unpredictable. 

Clinical evaluation and biochemical parameters {#Sec4}
----------------------------------------------

Global biochemical parameters *i.e*.  urea, creatinine, TGO (AST, aspartate aminotransferase), GGT (gamma-glutamyl transferase), BT (total bilirubin), BD (direct bilirubin), BI (indirect bilirubin), FAL (alkaline phosphatase) and ALB (Albumin) were all within normal limits during and after treatment termination.  We have observed no adverse effects in the animals used in our study (Table [S2](#MOESM1){ref-type="media"}). 

Histological analysis {#Sec5}
---------------------

Histological analysis of untreated and treated tumors showed proliferation of atypical and hyperpigmented melanocytes, presence of numerous macrophages phagocytizing melanin, but only a few perivascular lymphocytes (Fig.  [3](#Fig3){ref-type="fig"}).  Tumor treated for 28 days (animal 797, tumor 6), showed regression areas represented by the absence of atypical melanocytes, only remaining melanophages (Fig.  [3d](#Fig3){ref-type="fig"}).Figure 3*In vivo* histology.  Distinct histopathological slides are presented. (** a**) animal 438, tumor 2, vehicle excised on day 28 -- an untreated tumor showing melanophages and many atypical melanocytes. (** b**) and (**c**) animal 797, tumor 3, vehicle excised at day 28 -- an untreated tumor showing less atypical melanocytic cells. (** d**) animal 797, tumor 6, treated with Amblyomin-X and excised at day 28 -- regression areas showing only melanophages can be seen. 

Bioinformatics {#Sec6}
--------------

### Mapping and differential expression {#Sec7}

Alignment and read count resulted in 269,991 transcripts, displayed in a matrix with horse Ensembl IDs as rows and 18 samples as columns, of which we removed two with low library size, both related to 12 h samples.  We also excluded lowly expressed genes.  We used edgeR to calculate the normalized expression table in CPM and BioMaRt to find horse-to-human ortholog coding genes, which yielded 14.414 transcripts, of which 13,138 transcripts had valid gene symbols for horse and 13,943 for human.  There were more symbols for human genes due to paralog-ortholog conversion.  For horse transcripts, there were 546 DEGs for 6hx0h and 259 DEGs for 12hx0h (Table [S3](#MOESM1){ref-type="media"}).  Those DEGs having human orthologs were 546 for 6hx0h and 259 for 12hx0h.  We could not find DEGs for the 12hx6h comparison, probably because gene expressions between these time points are somehow very close. 

### Pathway analyses {#Sec8}

Metacore pathway enrichment analysis found 196 pathways for 6hx0h and 67 pathways 12hx0h.  The same analysis can be seen using Reactome, String-db with KEGG, and String-db with GO Biological Processes (Table [S4](#MOESM1){ref-type="media"}). 

We classified the enriched pathways according to the words in each pathway description name that may represent a function or have a gene symbol (Fig.  [S1a](#MOESM1){ref-type="media"}).  We also ranked these pathways regarding how strongly modulated they were (delta of the sum of LFC = sum of LFC for 6hx0h minus sum of LFC for 12hx0h, for all genes in the pathway), resulting in positive values for pathways more strongly modulated for 6hx0h compared to 12hx0h, negative if more strongly modulated for 12hx0h compared to 6hx0h, and similar if the modulations were close to each other (Table [S5](#MOESM1){ref-type="media"}).  For instance, the Innate Immune System (IIS) class comprised 26 pathways, of which 19 had genes more strongly modulated for 6hx0h, 3 for 12hx0h, and 4 were equally modulated (Fig.  [S1b--d](#MOESM1){ref-type="media"}).  According to the pathway classification, we have found that "immune system" was the class with more enriched pathways (80 pathways), according to Metacore, followed by "signaling" (55 pathways), "lung" (37 pathways), "canonical pathway" (31 pathways), "innate immune system" (26 pathways), "inflammation" (23 pathways), "cancer" (22 pathways), "muscle" (18 pathways), and "adhesion/ECM/cytoskeleton" (15 pathways). 

Figure [S1b](#MOESM1){ref-type="media"} shows 6hx0h higher modulated pathway classes, where immune system, pathway signaling, (lung) inflammation, innate immune system, canonical pathway, cancer, disease, migration/invasion/motility, adhesion/ECM/cytoskeleton, apoptosis, muscle, immunogenic cell death, virus responses, and stress, stand out as the most important pathway classes for this condition. 

Figure [S1c](#MOESM1){ref-type="media"} shows 12hx0h higher modulated pathway classes, however only a few pathways can be seen such as adhesion/ECM/cytoskeleton, lung, muscle, immune system, receptor pathways, signaling pathways, innate immune system, migration/invasion/motility, and cancer, resulting in a weaker response for 12hx0h when compared to 6hx0h. 

Finally, fig.  [S1d](#MOESM1){ref-type="media"} shows classes with similar responses for 6hx0h and 12×0h, including immune system, pathway signaling, canonical pathways, development, complement, muscle, cancer, dermis, and innate immune system. 

In summary, most of the pathways were related to cell stress, innate immune response/inflammation, cytokines, cell damage/endothelial, neutrophil, apoptosis, ECM, cell adhesion, muscle, T-help, complement, and coagulation.  Therefore, pathway classes were divided into three groups: a) first responses (ER-stress and IIS); b) confounding factors; and c) secondary responses (Table [1](#Tab1){ref-type="table"}).  Herein, confounding factors are denoted as enriched pathways, and part of the expression of their modulated genes, observed not probably due to the drug action, but mainly due to the stress and wounds caused by the injections.  Next, we summarize the results highlighting pathways.  In the text below, all pathway classes are written in bold. Table 1Selected enriched pathways that represent genes and important functions related to the *ex vivo* experiment.a -- first responses**FunctionalPathwayDEGs(6hx0h)DEGs(12hx0h)**ER-stressHCV-mediated liver damage and predisposition to HCC via cell stressCALR, CASP3, CYCS, HSP90B1, HSPA5, HSPD1, IL1R1, IL6, IL6ST, ITPR1, SOD2, STAT3, XBP1Endoplasmic reticulum stress response pathwayCYCS, HSP90B1, HSPA5, ITPR1, SOD2, XBP1Role of PKR in stress-induced antiviral cell responseCASP3, EIF2AK2, IL1B, IL1R1, IL6, IL8, STAT1, TLR2EIF2AK2, IL1B, IRF7, STAT1, TLR2TLR (IIS)HSP60 and HSP70 - TLR signaling pathwayHSPA5, HSPA6, HSPA8, HSPD1, IL1B, IL6, IL8, TLR2RLR (IIS)Innate immune response to RNA viral infectionDDX58, DHX58, IFIH1, IRF7OSM (IIS)Oncostatin M signaling via JAK-StatIL6ST, OSMR, SERPINA3, SOCS3, STAT1, STAT3, TIMP1SOCS3, STAT1, TIMP1Oncostatin M signaling via MAPKIL6ST, LDLR, MMP3, OSMR, STAT1, TIMP1OAS/RNase L (IIS)Antiviral actions of InterferonsEIF2AK2, IRF9, OAS2, STAT1, WARSEIF2AK2, IRF9, MX1, OAS1, OAS2, STAT1**b -- confounding factorsFunctionalPathwayDEGs(6hx0h)DEGs(12hx0h)**NeutrophilInhibition of neutrophil migration by proresolving lipid mediators in COPDACTN2, C5AR1, CD34, CXCR2, FPR2, IL1B, IL1R1, IL8, ITPR1, PRKCQ, PTAFR, TLR2ACTN2, CXCR2, IL1B, PRKCQ, TLR2EosinophilEosinophil chemotaxis in asthmaC3, CCL7, CXCL10, CXCR2, NGFCell adhesionECM remodelingIGFBP4, IL8, MMP3, PLAUR, SERPINE1, TIMP1Integrin inside-out signaling in neutrophilsCXCL1, CXCR2, IL8, ITPR1, LYN, PTAFR, RASGRP2, SELE, SELPMuscleMuscle contraction: Delta-type opioid receptor in smooth muscle contractionITPR1, MYH8, PENKMYH8, PENK**c -- secondary responsesFunctionalPathwayDEGs(6hx0h)DEGs(12hx0h)**InflammationRelease of pro-inflammatory factors and proteases by alveolar macrophages in asthmaCXCL1, IL1B, IL6, IL8, MMP3, STAT1, TIMP1, TLR2CXCL10, IL1B, STAT1, TIMP1, TLR2Release of pro-inflammatory mediators and elastolytic enzymes by alveolar macrophages in COPDCTSL, IL1B, IL6, IL8, STAT1, TLR2CTSL, CXCL10, IL1B, STAT1, TLR2Inflammatory mechanisms of pancreatic cancerogenesisCD46, CXCL1, CXCR2, IL1B, IL1R1, IL6, IL8, STAT1, STAT3IL-1IL-1 signaling pathwayBIRC3, CCL7, CXCL1, IL1B, IL1R1, IL6, IL8, PTGES, STAT1, ZC3H12AThe innate immune response to contact allergensHSPA5, HSPA6, HSPA8, IL1B, IL1R1, IL6, IL8, TLR2Cell migrationChemotaxis: CCL16-, CCL20-, CXCL16- and CCL25-mediated cell migrationITPR1, MMP11, MMP3, MMP8Cytoskeleton remodelingCytoskeleton remodeling: Regulation of actin cytoskeleton organization by the kinase effectors of Rho GTPasesACTN2, MYH8ACTN2, MYH8Apoptosis & SurvivalEndoplasmic reticulum stress response pathwayCYCS, HSP90B1, HSPA5, ITPR1, SOD2, XBP1Apoptosis and survival: CXCR3-B signalingCASP3, CYCS, HMOX1, ITPR1, RYR1, STAT3Apo-2L(TNFSF10)-induced apoptosis in melanomaCASP3, CASP4, CYCS, IL8, XBP1Interferon-alpha/betaIFN-alpha/beta signaling via JAK/STATEIF2AK2, IFIT3, IL1RN, IRF9, ISG15, SOCS3, STAT1, STAT3CXCL10, EIF2AK2, IFIT3, IFITM1, IL1RN, IRF7, IRF9, ISG15, ISG20, OAS1, SOCS3, STAT1, XAF1IFN-alpha/beta signaling via MAPKsEIF2AK2, IFIT3, IRF9, ISG15, PRKCQ, STAT1, ZBTB16CXCL10, EIF2AK2, IFIT3, IRF7, IRF9, ISG15, PRKCQ, STAT1The four columns are the functions, pathways, enriched DEGs for 6hx0h comparison and for 12hx0h comparison, respectively.  The table was divided into A, B, and C, for pathways classified as "ER-Stress and IIS", "confounding factors", and "secondary responses", respectively. 

**ER-stress** was inferred due to the enriched pathways such as "HCV-mediated liver damage and predisposition to HCC via cell stress", "Endoplasmic reticulum stress response pathway", and "Role of PKR in stress-induced antiviral cell response".  For **Immune System**, six innate immune system responses were found: TLR pathway, RIG-I-like Receptors (RLR), Oncostatin-M pathway, OAS/RNase L pathway, Neutrophil, Eosinophil pathways. 

**TLR** was inferred due to the enriched "HSP60 and HSP70 - TLR signaling pathway"; **RLR** due to the enriched "Innate immune response to RNA viral infection" pathway; and **Oncostatin-M** due to the enriched "Antiviral actions of Interferons" pathway. 

For the present study, **Neutrophil, Eosinophil** and **Complement** are termed as confounding factors, since these biological functions are related to the damage caused by the intratumoral injections, activating "endothelial damage" and "wound healing" pathways.  Another confounding factor is **Muscle contraction**, to be discussed later.  All data related to these pathways can be seen in Table [1B](#Tab1){ref-type="table"}. 

In addition, **ER-stress** and **Innate Immune Responses** should activate many secondary responses like Inflammation, IL-1, Cell migration, Cytoskeleton remodeling, Apoptosis & Survival, and Interferon-alpha/beta. ** Inflammation** was inferred due to the enriched pathways related to macrophages: "Release of pro-inflammatory factors and proteases by alveolar macrophages in asthma", "Release of pro-inflammatory mediators and elastolytic enzymes by alveolar macrophages in COPD", "Inflammatory mechanisms of pancreatic cancerogenesis", "Immune response_IL-6-induced acute-phase response in hepatocytes", and "Immune response_TREM1 signaling pathway"; and **IL-1** was inferred due to enriched pathways: "IL-1 signaling pathway" and "The innate immune response to contact allergens". 

Moreover, **Cell Migration** was inferred due to the enriched "Chemotaxis: *CCL16*-, *CCL20*-, *CXCL16*- and *CCL25*-mediated cell migration" pathway; and **Cytoskeleton Remodeling** was inferred due to the enriched "Cytoskeleton remodeling: Regulation of actin cytoskeleton organization by the kinase effectors of Rho GTPases" pathway.  Finally, **Apoptosis & Survival** was inferred due to enriched "Endoplasmic reticulum stress response pathway", "Apoptosis and survival: CXCR3-B signaling", and "Apo-2L(TNFSF10)-induced apoptosis in melanoma" pathways; and **Interferon-alpha/beta** was inferred due to the enriched "Immune response: IFN-alpha/beta signaling via JAK/STAT" and "Immune response: IFN-alpha/beta signaling via MAPKs" pathways. 

Network analyses {#Sec9}
----------------

With the application of String-db, 93 enriched pathways for 6hx0h and 59 for 12hx0h were found.  String-db allows for the calculation of protein-protein interaction (PPI) networks, and with igraph^[@CR37]^ it was possible to figure out the connectivity index (k) and betweenness centrality (g) (Table [S6](#MOESM1){ref-type="media"}).  Next, Gephi was used to cluster genes according to the modularity index (community detection algorithm) with a resolution of 1.2 and the use of weights.  The result can be seen in Table [S7](#MOESM1){ref-type="media"} and respective networks can be seen in Fig.  [4](#Fig4){ref-type="fig"}, where the main enriched pathways could be achieved using all clustered DEGs and Reactome to find the enriched pathways.  Node colors were obtained using modularity classes for nodes clusterization, their sizes are proportional to the degree of connectivity and the label size is proportional to the betweenness centrality. Figure 4Network analysis.  Gephi clusterized DEGs according to the modularity index with resolution of 1.2 and using weights, and respective enrichment analyses were done using Reactome.  On the left (**a**) are all DEGs related to 6hx0h comparison, and 9 clusters could be defined related to Immune System, Transport & Cytochrome C, Antiviral pathways (RLR), GPCR pathway, UPR, rRNA processing & metabolism of nucleotides, Cell junction & cell-cell communication, muscle & ion homeostasis, and Circadian pathway & Ubiquitination.  On the right (**b**) are DEGs related to 12hx0h comparison, resulting in a smaller network, with seven clusters: Innate Immune System responses, Cell cycle pathways, Neutrophil Degranulation, ESR-mediated signaling, Muscle, Circadian Clock, and rRNA processing. 

*In-silico* crosstalk simulation {#Sec10}
--------------------------------

### Crosstalk between pathways {#Sec11}

Five different possible crosstalks were simulated: a) adhesion-ECM-cytoskeleton versus remodeling versus stress, b) ICD versus apoptosis versus ER-chaperone-Golgi, c) cancer versus complement versus inflammation, d) IIS versus apoptosis" versus autophagy, and e) angiogenesis-vascular versus hypoxia versus stress.  For this purpose we used Circlize R package. 

All five crosstalks showed a strong interaction between genes related to different biological functions; additionally, there were genes that participate simultaneously in different functions (pleiotropy).  Only autophagy was not seen, possibly because it is not based on a phenomenon of transcriptional modulation.  Crosstalk analysis between Innate Immune System x Apoptosis x Autophagy can be seen in Fig.  [5](#Fig5){ref-type="fig"}, and respective enriched gene modulations in Fig.  [6](#Fig6){ref-type="fig"} (other crosstalks can be seen in Fig.  [S2](#MOESM1){ref-type="media"}, and respective heatmap, dotplot, and spaghetti plots in Fig.  [S3](#MOESM1){ref-type="media"}).  There is a stronger interaction between genes for 6hx0h compared to 12hx0h.  All circlize plots have more genes and are more connected for 6hx0h compared to 12hx0h, demonstrating a strong response at 6 h, diminishing at 12 h.Figure 5Crosstalk between three pathway classes.  To better visualize the crosstalk between DEGs, we downloaded the String-db PPI table and filtered only interactions with scores greater than or equal to 0.4.  On the left (**a**) is the crosstalk between DEGs related to innate Immune system, apoptosis and autophagy, with all 59 DEGs well connected, which leads to the inference that IIS and Apoptosis have common genes to work properly, or to be activated.  Only autophagy was a less significant pathway, not observed until 12 h by transcriptome analysis.  We can see on the right (**b**) all 28 DEGs for 12hx0h, with fewer genes and interactions comparing to 6hx0h.  Near each gene are its possible pathways like 'iis' (innate immune system), 'ap' (apoptosis), and 'aut' (autophagy).  In supplementary material other four distinct crosstalks can be seen (Fig.  S2).Figure 6Heatmap and dotplot for the innate immune system (**a,b**), apoptosis (**c,d**), and autophagy (**e,f**).  There are 59 genes related to the Innate Immune System and most of them were upregulated.  Here, hierarchical cluster analysis shows 6 h samples mixed with 12 h samples, but completely apart form 0 h. Genes like *CXCL10, IRF7, OAS21, IFH1, IFITM1, MX1, OAS2, XAF1, DDX58, IRF9*, and *ISG20* are more upregulated at 12 h than 6 h, meaning that RLR and OAS pathways responses are increasing in time.  Genes like *MYL1, MYH8*, and *ACTN2* are downregulated, the first two with higher modulation at 6 h, and the latter at 12 h. Apoptosis split very well the 3 time points, besides sample A4T2h6 sample, related to 6 h, is closer to 12 h. Many apoptotic genes were upregulated, except *RYR1*, with a high modulation for *SOD2*, *HSPA6* and *IL1β*.  Genes related to autophagy could not enrich the "autophagy pathway", but *IL6, IL1β, CASP3, IL6ST, ALPL*, and *IL8 were* upregulated DEGs, and *MYL1* and *MYH8* were downregulated DEGs. 

### CEMiTool/WGCNA coexpression analysis {#Sec12}

Performing coexpression of genes modulated by Amblyomin-X treatment using CEMiTool, yielded 33 modules named M1 to M32, and an additional one called not-correlated (Table [S8](#MOESM1){ref-type="media"}).  Enrichment analysis was calculated on these modules using KEGG, Reactome, WikiPathways and Gene Ontology (GO), which resulted in a predominance of signaling pathways related to muscle, interleukins, and response to virus infection.  Interestingly, some of the genes enriched for these pathways were positioned just after the DEG definition threshold. 

For instance, 7 genes (*ISG15*, *CASP10*, *CASP8*, *TNF*, *DHX58*, *TRAF3*, and *IRF7*) were found, related to KEGG RIG-I-like receptor signaling pathway/RIG-I-like Receptor Signaling WP3865 at M13, although only 3 (*ISG15*, *IRF7* and *DHX58*) out of those 7 were DEGs for 12hx0h and only one (*ISG15*) for 06hx0h.  Since these genes are clustered in the same module, they showed the same expression pattern. 

In this same module M13, seven genes (*IRF1, CASP10, CASP8, TNF, CDKN2A, TRAF3*, and *IRF7*) were related to apoptosis and suggesting that Amblyomin-X has triggered this pathway. * HSH2D*, coding for a protein that may be a modulator of the apoptotic response through its ability to affect mitochondrial stability^[@CR38]^, was one of the main intramodular hubs in this module.  Another intramodular hub in module M13 is *ORM1*, coding for a member of the lipocalin family of proteins, with a function in modulating the innate immune system in events of acute phase inflammation, and known to bind synthetic drugs^[@CR39]^. *ORM1* is also a DEG for 06hx0h and 12hx0h comparisons.  The interaction hub in module M13 was *SMURF1*, which promotes ubiquitination and subsequent proteasomal degradation of *MAVS*^[@CR40]^ and plays a role in dendrite formation by melanocytes^[@CR41]^.

In module M9, 13 genes (*H2AFZ, HLA-B, HLA-C, HLA-F, PSMA7, PSMB1, PSMB2, PSMB3, PSMB6*, *PSMC3, PSMC4, PSMD8*, and *UBE2D2*) were found, enriched by WikiPathways for "Proteasome Degradation" (WP183 from WikiPathways), and 16 (*ASB1*, *ASB8*, *KEAP1 PSMA7*, *PSMB1*, *PSMB2*, *PSMB3*, *PSMB6*, *PSMC3*, *PSMC4*, *PSMD8*, *RNF138*, *RNF217*, *UBE2D2*, *UBE2M*, and *TRAF7*) enriched by Reactome for "Antigen processing: Ubiquitination & Proteasome degradation".  None of which are DEGs, but still following the same pattern of increase in expression over time.  We expected such a phenomenon with Amblyomin-X treatment^[@CR19]^.

In module 9, the major intramodular hub is *PRELID1*, involved in the modulation of the mitochondrial apoptotic pathway by ensuring the accumulation of cardiolipin (CL) in mitochondrial membranes.  One of the interaction hubs in this module is *UBE2M*, one of the ubiquitin-conjugating enzymes.  The encoded protein is linked with the ubiquitin-like protein *NEDD8*, which can be conjugated to cellular proteins, such as Cdc53/Culin. 

### Cellular interaction partners of Amblyomin-X {#Sec13}

Identification of cellular interacting partners is critical to understand the initial process originated by Amblyomin-X, resulting in modulation of pathways at the later stage of the route.  Therefore, we obtained the equine melanoma cellular interactomics profile for Amblyomin-X by co-precipitation.  Multiple potentially interacting proteins were identified in eluent fraction that bound to immobilized Amblyomin-X. The list of proteins is shown in Table [S9](#MOESM1){ref-type="media"}, including unique identifiers and the number of peptides.  Interestingly, we found proteins related to biological effect previously reported in Amblyomin-X studies, such as apoptosis, mitochondrial dysfunction, regulation of cell migration and protein clearance, which includes cytochrome-c, plasminogen, actin cytoplasmic 1, fibronectin, heat shock protein HSP 90-alpha, E3 ubiquitin-protein ligase Mdm2, mitochondrial superoxide dismutase, and vimentin^[@CR15],[@CR18],[@CR42],[@CR43]^. Furthermore, we identified immunogenic proteins such as Toll-like receptor 2 (*TLR2*) and T-cell surface antigen CD2 (*CD2*) as interacting partners of Amblyomin-X.

### Validation of RNA-Seq results through qRT-PCR {#Sec14}

This study selected 56 DEGs related to innate immune response, apoptosis and inflammation for validation with qRT-PCR.  Next, these genes were classified into 3 groups, according to the median expression: highly expressed (one of the median expressions greater or equal to 50 CPM), moderately expressed (not lowly neither highly expressed) and lowly expressed (all expressions less than 5 CPM) (Table [S10](#MOESM1){ref-type="media"}).  Forty-two (42) out of these 56 genes were chosen for qRT-PCR validation (Fig.  [S4](#MOESM1){ref-type="media"}).  The correlation between qRT-PCR LFC for 6hx0h and RNA-Seq LFC for 6hx0h (Table [S11](#MOESM1){ref-type="media"}) was 80.2%, and the correlation between qRT-PCR LFC for 12hx0h and RNA-Seq LFC for 12hx0h was 81.5% (Fig.  [S5](#MOESM1){ref-type="media"}). 

Discussion {#Sec15}
==========

Amblyomin-X treatment reduced equine melanoma tumor size, after which we investigated modulated genes and enriched pathways.  The study considered as confounding factors: the pathways related to needle tissue injury as wound healing, the increase in tumor volume as hypoxia, and the vein injury as vascular endothelial cell damage.  Although real, these processes were probably not the reason for the decrease in tumor volumes since many tumors treated with vehicle did not show a decrease in size. 

Our results validated previous uses of Amblyomin-X against several tumor cell lines including skin melanoma^[@CR15]--[@CR18],[@CR42]^ like ER-stress, Mitochondria-stress, and Apoptosis.  It is important to note that the transcriptome results presented herein are the mean cell expression values in the tumor's core, a complex system composed of many distinct normal cells surrounding tumor cells.  It is also necessary to observe that pleiotropic effects can deceive the expression of some DEGs related to these functions.  For instance, Actins and Calcium ion play an important role in pathways such as Endoplasmic Reticulum Stress (ER-stress), Cytoskeleton remodeling, and Muscle contraction, and IL-8 plays a role in endothelial damage and inflammation response.  Since we are interested in the orchestrated responses between Immune System, Inflammation, Apoptosis, ER-stress and Cytoskeleton remodeling, any actin or calcium ion release, belonging only to muscle processes, would confound our inferences.  Taken the above into consideration, we presented the most statistically significant genes and pathways modulated by Amblyomin-X (see also supplementary material). 

Figure [7](#Fig7){ref-type="fig"}, a simplified model of the tumor microenvironment, shows the Amblyomin-X molecule in the extracellular medium.  That allows us to hypothesize on how some pathways are interrelated and evolve in time.  Amblyomin-X is able to enter cancer cells that expose phosphatidylserine on the external membrane, entering by endocytosis, and promoting ER-stress through the proteasome inhibition (PI) or other unknown processes.  Next, a series of related phenomena are observed, like Unfolded Protein Response (UPR), Mitochondria-stress and Apoptosis.  We also hypothesize that in parallel or due to the cell stress, a second phenomenon of exosomes transport^[@CR44]^ occurs between these cells or, less likely, the activation of retrovirus transcription in cancer cells^[@CR45]^. Both processes release viral RNA that interact with viral RNA sensors (*RIG-I, LGP2, MDA5*, and *DDX60*) activating the RLR pathway, which has three outcomes: a) apoptosis through *PMAIP1* (*NOXA*); b) interferon transcription, through *TBK1, IRF3*, and *IRF7;* and c) cytokines transcription (*IL1β, IL-6, IL-8, IL12, IP10*, and *TNF-α*) through *NF-κB*.  Related to the NF-κB path, *IKB* must be degraded by the proteasome in order to *NF-κB* to migrate to the nucleus, however, we just hypothesized that the proteasome was inhibited.  Taking this into account, the transcription of these cytokines should decrease in melanoma cells.  However, other cells, like macrophages, will be able to transcribe these cytokines after RLR signaling activation.  Therefore, some melanoma cells will be unable to sustain the cytokine transcription, but paracrine signaling will be present and RLR sensors transcription will be boosted, as well as the inflammatory response, and probably M1 and M2 polarization starts.  As mentioned before, these signaling crosstalk should increase the inflammatory response, and also *TREM1* and IL-6-induced acute-phase response (downstream *ORM1, SERPINE1, CXCL1, and C3)* pathways will be activated. Figure 7Schematic overview of a simplified tumor environment model.  Melanoma, Macrophage, and Fibroblast (stroma) are possibly the main cells interacting through cytokines, exosomes, and other signaling proteins.  Herein, we postulate that some Melanoma cells enter in ER-stress state because Amblyomin-X was internalized, next mitochondria stress and cytoskeleton remodeling occur, and finally, apoptosis and survival duel are the final outcomes.  Concomitant, IIS pathways are activated in melanoma cells, first by the RLR responses, possibly sensing RNA molecules originating from fibroblast exosomes^[@CR44]^ or being transcribed from endogenous retroviral elements^[@CR45]^. IFN type I and NOXA transcriptions were not seen, but many cytokines (*IL1β, IL-6, IL-8, IL12, IP10*, and *TNF-α*) were DEGs, and their proteins were likely to have been produced and released to extracellular medium.  Many macrophage signaling pathways were enriched, and their transcribed cytokines are *IL1β, IL-6, IL-8, CCL2, CXCL1, CXCL10 (IP10)*, and *TNF-α*, all DEGs.  Possibly the production of some of these proteins results in M1 and M2 polarization, since *STAT1, STAT3*, and *SOCS3* were also DEGs.  These macrophage cytokine releases should lead to many processes, including a crosstalk to Fibroblast, but mainly, the feedback to RLR pathway and the inflammatory response, see text.  All these inferences were made according to Metacore algorithms, references, and database. 

Interesting evidence was the simultaneous activation of different innate immune system pathways, increasing the expression of many DEGs related to RLR, OAS, TLR and Oncostatin-M pathways.  Observing Fig.  [6a,b](#Fig6){ref-type="fig"}, most of the DEGs related to IIS have a maximum fold-change at 6 h, decreasing their expressions at 12 h, but still remaining DEGs; besides *CXCL10, IRF7, OAS1, IFIH1, DHX58, MX1, OAS2, XAF1, DDX58, IRF9, ISG20, CCL7* (or *CCL2*, ortholog), *TIMP1, ISG15, ZBTB16, MYH8, RASGRP2*, and *CD34* increasing their expression after 6 h, perhaps boosting the IIS response.  In Fig.  [6c,d](#Fig6){ref-type="fig"}, most of the genes related to apoptosis have also a higher expression at 6 h; therefore, we cannot infer causality between IIS and Apoptotic genes and pathways. 

We hypothesize that *PKR (EIF2AK2)* is a central hub related to the innate immune system pathways, supporting the transcription of inflammatory cytokines.  As mentioned before, macrophages continuously induce the transcription of *IL1β, IL-6, IL-8, CXCL1*, *CXCL10*, and *TNF-α*, the last of which was moderately expressed. * CXCL10* probably stimulates monocytes, NK cells, and T-Cells or plays a role in modulating adhesion genes.  Interferons type 1, *IFNα* and *IFNβ*, should be transcribed through one of the RLR pathway-branches - via *TBK1, IRF3*, and *IRF7 -* but we could not observe any transcript expression^[@CR46],[@CR47]^. It is also important to highlight that *TMEM173* (*STING*) was moderately expressed and not modulated, and *PMAIP1* (*NOXA*), an important gene related to the Bortezomib pathway, was not transcribed.  These results demonstrated that Amblyomin-X follows different routes. 

Another important pathway was Cytoskeleton remodeling, enriched for 6hx0h.  Interestingly Öhman *et al*.^[@CR48]^ described the crosstalk between many RLR proteins and Actin and Tubulin proteins close to the mitochondria.  In our results *ACTN1* and *ACTR3* are upregulated DEGs; *ACTA1* and *ACTN2* are highly downregulated DEGs; and regarding Tubulin family, *TUBB* and *TUBB6* are upregulated DEGs. 

"Oncostatin M signaling via JAK-Stat" or "Oncostatin M signaling via MAPK" pathways are related to the ligand *OSM* (Oncostatin M) and their receptors. * OSM* is a cytokine member of the leukemia inhibitory factor/oncostatin-M (*LIF*/*OSM*) family and a growth regulator with important roles in inflammation and tumor growth inhibition. * OSM* effector was not a DEG showing low expression, but its receptors were DEGs (*OSMR* and *OSM receptor*).  Downstream, *TIMP1*, a Metalloproteinase inhibitor, seems to inhibit *MMP1* trying to stop ECM remodeling; *CCL2* (or *CCL7*, its horse ortholog gene) increases inflammation signal and *SERPINE3*, a plasma protease inhibitor, probably regulates inflammatory response. 

Amblyomin-X modulates ER-stress, mitochondria dysfunction, and apoptosis {#Sec16}
------------------------------------------------------------------------

Many apoptotic signals were high at 6 h (*CASP3*, *CASP4*, and *CYCS*) and their expression lowered at 12 h. Survival signals showed the same behavior (*BIRC2*, *BIRC3*, *ATF6*, and *SOD2*); besides; *SOD2* (up for 6hx0h, and up for 12hx0h), was the only pro-survival signal which is upregulated at both times.  At 6 h cytochrome-c (*CYCS*) was upregulated, which shows a possible mitochondrial damage.  Calpain 2 (*CAPN2*) and *HSP60*, the latter was a DEG, and *BAX* (verified in previous *in vitro* experiments) support this hypothesis.  We assume that ER-stress shows a peak close to 6 h since *GRP78*, *IP3R1*, *ATF6*, *XBP1*, *PDIA6*, and *HSP90B1* were upregulated at this time.  Since the treatment can kill some cells and allows others to survive, we hypothesized that it also activates immunogenic cell death (ICD) mechanisms via ER-stress with the release of *HMGB1* and *CALR* (calreticulin), an important DEG at 6 h, which can be found close to ER, in the inner side of the cell membrane, and also on the external side.  According to Jheng *et al*.^[@CR49]^, we could map autophagy upstream signals related to ER-stress with the following genes: *CASP4*, *CASP12*, *JNK*, *ATF6*, *CHOP*, *ATF4*, *EIF2AK3* (*PERK*), *EIF2A*, and *GADD34*.  Many of them were well expressed, but not modulated.  Therefore, autophagy may start later, since it was not enriched, or may even not start (see Fig.  [6](#Fig6){ref-type="fig"}).  Since Amblyomin-X can inhibit the proteasome machinery, some proteins accumulate nearby ER and new chaperones were transcribed (*HSPD1* (*HSP60)*, *HSPA4* (*HSP70)*, *HSP90AA1* (*HSP90*)) resembling what the UPR, an important mechanism belonging to ER-stress pathway.  HMGB1-RAGE signaling pathway was also enriched, indicating that it is possibly another effector inducing the transcription of *IL6*, *IL1β*, *CXCL8 (IL8*), *IL1RN*, *SERPINE1* (*PAI1*) and *CHGA* (Chromogranin A).  The data was supported by the interactomics profile of Amblyomin-X, in which the binding partners identified were proteins such as cytochrome-c, plasminogen, actin cytoplasmic 1, fibronectin, heat shock protein HSP 90-alpha, E3 ubiquitin-protein ligase, mitochondrial superoxide dismutase, and vimentin. 

In conclusion, Amblyomin-X is an antitumor molecule because it acts on multiple cellular targets to induce tumor cell death.  It selectively enters only cancer cells via endocytosis, binding to the external membrane attracted by phosphatidylserine affinity, and thereafter, inside the cell, modulates its microenvironment.  Amblyomin-X activates the ER- and mitochondria-stress pathways followed by apoptosis and survival responses.  There is also a crosstalk between proteasome inhibition and ER-stress that must be further studied.  The activation of the innate immune pathways inducing pro-inflammatory cytokines and chemokines was observed, but the initial cause is still unknown.  Interferons type I transcription were expected to be seen, but an unknown process inhibits it.  We believe that a more detailed data acquisition between 0 h and 6 h will elucidate some processes and their causes.  In summary, the horse melanoma tumor model shed light on new genes, pathways, and crosstalks among them, showing to be a useful tool in which our new outcomes must be validated, in the future, in human melanoma experiments. 

Methods {#Sec17}
=======

The whole experiment was divided in: a) *in vivo* melanoma tumor treatment for a month and follow-up with clinical lab tests for up to five months; b) excision of *in vivo* melanoma tumors for histology; c) e*x vivo* excision melanoma tumors after 0 h (control), after 6 h treatment, and after 12 h treatment followed by interactomics and proteomics (Fig.  [S6](#MOESM1){ref-type="media"}). 

Experimental conduct and ethical compliance {#Sec18}
-------------------------------------------

Horses were housed and treated on the Sao Joaquim farm (Araçariguama, SP, Brazil) of the Butantan Institute (Sao Paulo, SP, Brazil).  All animal procedures were previously approved by the internal ethics committee of the Butantan Institute (CEAUIB/BUTANTAN, protocol numbers 1999061115 and 8644050716) and were conducted in compliance with Good Clinical Practice guidelines (European Union Committee for Veterinary Medicinal Products), with the Brazilian National Council of Animal Experimentation Control (CONCEA, DBCA v.12/2016), and with the Brazilian Ministry of Agriculture, Livestock and Supply (MAPA) regulations.  Animal husbandry was in compliance with the requirements of national legislation and local animal welfare guidelines (Brazilian Federal Law N° 11794/08).  These studies were conducted under veterinary supervision, with veterinary attention available at all times. 

Expression and purification of recombinant Amblyomin-X {#Sec19}
------------------------------------------------------

Recombinant Amblyomin-X was produced using an *E. coli* expression system as described previously^[@CR50]^.

*In vivo a*nimal experimental conditions and tumor treatment {#Sec20}
------------------------------------------------------------

Four male gray horses were the subjects used in the *in vivo* experiments herein studied.  The horses, aged between 10 and 20 years, and weighing between 320 and 430 kg, were housed in stables with access to food (ration and green) and water *ad libitum*.  Each animal was submitted to a general clinical exam before its inclusion in the experiment.  Tumor sizes were recorded and melanoma lesions, which appear spontaneously at the perianal and ventral tail regions, were mapped.  Tumor dimensions were measured with a caliper and tumor volumes calculated by approximating a sphere volume, *i.e*.,  using the formula for sphere volume = 4/3 π x radius 1 x radius 2 x radius^3^ (mm^3^)^[@CR51]^. Deeply infiltrating tumors were rejected for treatment due to possible measurement errors.  Volume measurements and excised tumor weight were used to estimate the mean tumor density value, which in turn were used to calculate the Amblyomin-X doses as 1 mg/Kg of tumor mass. 

On 4 animals, we selected about 9 tumors, of which 3 were injected with Amblyomin-X every 3 days over a period of 28 days; 3 were injected with the same volume of PBS as a vehicle as often; and 3 were not treated (Table [S12](#MOESM1){ref-type="media"}).  Some animals had only 8 tumors available, of which we selected 3 for treatment, 2--3 for control, and the remaining 2--3 for the vehicle. 

Clinical evaluation {#Sec21}
-------------------

Body weights were recorded immediately before treatment and then after every 15 days until the end of the experiment, using a load-bar scale (TOLEDO^®^ MGR-4000).  Blood samples were collected from the jugular vein for hematological and biochemical analysis at time zero and after every 3 weeks until the end of the protocol.  Hematological and biochemical parameters were obtained using automated analyzers (BC-2800Vet, Midnray and Cobas-Mira Plus, Roche, respectively).  Biochemical parameters included urea, creatinine, TGO (AST, aspartate aminotransferase), GGT (gamma-glutamyl transferase), BT (total bilirubin), BD (direct bilirubin), BI (indirect bilirubin), FAL (alkaline phosphatase) and ALB (Albumin).  Reagents used were obtained from Kovalent, Brazil. 

Tumor resection and preservation {#Sec22}
--------------------------------

For tumor resection, animals were sedated with 10% xylazine hydrochloride (0.53 mg/kg) and detomidine (10 µg/kg), followed by local anesthesia with 10 mL subcutaneous injection of lidocaine.  The incision site was shaved and disinfected it by wiping the skin surface with 0.2% chlorhexidine gluconate alcohol, followed by spraying the local with 79^o^ GL alcohol.  After the tumor had been completely removing, the cutaneous wound was sutured and cleaned daily with gentamicin sulfate solution (100 mg/mL) and application of a gentamicin based ointment (Vetagloss, Vetenil Inc., Brazil). 

Histological analysis {#Sec23}
---------------------

For histological analysis, excised tumors were fixed in 10% buffered formalin solution, dehydrated and enclosed in synthetic paraffin, followed by hematoxylin and eosin staining. 

*Ex vivo* experiment: transcriptome analysis {#Sec24}
--------------------------------------------

Nine horses were selected for the *ex vivo* experiment.  Melanoma tumors, in the ventral tail, were selected for transcriptome analysis.  Treated tumors were injected with a single dose of Amblyomin-X (1 mg/kg of tumor), being excised 6 hours and 12 hours after the treatment, and controls were untreated tumors being injected PBS at 0 h. For each time-point, two tumors were removed from three different animals, yielding six tumors per time-point for analysis.  The transcriptomic design can be seen on Table [S13](#MOESM1){ref-type="media"}. 

See supplementary material documentation for RNA extraction and library preparation, Quality and filtering fastq Reads, cDNA synthesis, and qRT-PCR

Protein interactions with Amblyomin-X {#Sec25}
-------------------------------------

Purified and lyophilized Amblyomin-X (6 mg) was dissolved in 1 ml of 0.2 M NaHCO~3~ containing 0.5 M NaCl, pH 8.3.  Proteins were immobilized on a HiTrap™ NHS-activated HP 1 ml column (GE Healthcare Life Sciences) according to the manufacturer's instructions.  After this, tissue extract from untreated equine tumor samples were applied to the HiTrap™ Amblyomin-X affinity column.  The column was washed with 60 ml of 20 mM Tris-HCl buffer, pH 8.3, to eliminate any non-specific protein interaction.  Bound proteins were eluted with 200 mM glycine containing 0.5 M NaCl, pH 4.0.  The eluent was extensively dialyzed with 25 mM ammonium bicarbonate and dried in a speed vacuum evaporator (Thermo Scientific).  Dried samples were stored at −80 °C or dissolved in 50 mM ammonium bicarbonate, containing 10 mM CaCl~2~ for MS/MS analyses.  For protein identification, trypsin hydrolysates samples were injected into an EASY-nano LC system (Proxeon Biosystems) coupled online to an ESI-LTQ-OrbitrapVelos mass spectrometer (Thermo Fisher Scientific), which was operated in a positive mode of ionization using the data-dependent automatic (DDA) survey MS scanned tandem mass spectra acquisition.  Peaks Studio 7.5 (Bioinformatics Solutions Inc. Canada) was employed for data acquisition, processing and analyses. 

Bioinformatics and systems biology {#Sec26}
----------------------------------

### Expression analysis {#Sec27}

In order to map and quantify transcripts, we used featureCounts from Subread package^[@CR52]^ resulting in a table of gene IDs versus samples with raw read counts.  Then, we calculated differentially expressed genes (DEGs) using the edgeR^[@CR53]^ package.  Genes with low expression (*i*.*e*.,  all sample groups having sum of counts less than one CPM, count per million in each) were filtered out.  DEGs were calculated with the following comparisons: a) 6 h after treatment (6 h) x control (0 h); b) 12 h after treatment (12 h) x control; and also c) 12 h x 6 h, here named as 6hx0h, 12hx0h, and 6hx12h, respectively.  Let define LFC as the log2 of the fold change between two conditions and FDR as the false discovery rate.  DEGs were defined as genes with absolute value of LFC \> = 1 and FDR \< 0.05.  All transcripts were mapped using horse Ensembl Gene ID, followed by BioMaRt tool^[@CR54]^ mapping orthologs between horse gene IDs and human gene IDs (*Homo sapiens* reference genome GRCh38 and annotation version 89). 

### Enrichment, network, and coexpression analysis {#Sec28}

Three different methods were used for *in-silico* enrichment analysis (Fig.  [S7](#MOESM1){ref-type="media"}).  Pathway Enrichment Analysis (PEA) which uses Over-Representation Analysis (ORA) statistics and databases like Metacore^[@CR55]^, Reactome^[@CR56],[@CR57]^, WikiPathways^[@CR58]^, KEGG^[@CR59]^ and Enrichr^[@CR60],[@CR61]^; Fast Gene Set Enrichment Analysis (fGSEA)^[@CR62]^, and String-db^[@CR63]^ version 11.0 enriching against KEGG, GO^[@CR64]^, and pFam^[@CR65]^.

For network analysis String-db, a PPI network, was used, setting the minimum score to 0.40.  Gephi version 0.9.2^[@CR66]^ was also applied, with DEGs as nodes and edge weights as PPI scores obtained from String-db.  In addition, circlize^[@CR67]^ was employed to simulate possible crosstalks between enriched pathways. 

WGCNA^[@CR68]^/CEMiTool^[@CR69]^ was used for coexpression analysis, as an unbiased technique based on gene-gene expression correlation, setting a minimum of 20 genes per module, and merging modules with correlation coefficient above 0.75.  We resorted to BioGrid database^[@CR70]^, including only genes with enough reads, in order to build a PPI network for each module.  In each module, intramodular and interaction hubs were respectively identified as coexpression and interaction. 

Supplementary information
=========================

 {#Sec29}

Supplementary Information.  S Table 01.  S Table 02.  S Table 03.  S Table 04.  S Table 05.  S Table 06.  S Table 07.  S Table 08.  S Table 09.  S Table 10.  S Table 11.  S Table 12.  S Table 13.  S Table 14. 

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 

Supplementary information
=========================

is available for this paper at 10.1038/s41598-020-63275-2. 

This project was supported by Fapesp grant numbers 2015/50040-4, 2016/13588-4, and 2017/02052-9, São Paulo Research Foundation and GlaxoSmithKline grant 2015/50040-4, BNDES grant 13.2.0711/2013, FAPESP-CeTICS grant 2013/07467-1, the National Council of Technological and Scientific Development (CNPQ) grant 303197/2017-0, and Coordination for the Improvement of Higher Education Personnel (CAPES). 

F.L. and C.E.M.T. -- bioinformatics, system biology, writing and discussion; A.I. and K.L.P.M. -- interactomics, writing and discussion; S.E.A.L.W. and R.V.B. - *in vivo* experiments; C.D.P. - qRT-PCR, RNA sample preparation and discussion; M.B.G. -- sample preparation, writing and discussion; R.C. - project idea, histological analysis, discussion and writing revision; J.G.S. - discussion; R.J.M.N. - histological analysis; F.I.A.A. - bioinformatics and discussion; J.M.A. -- bioinformatics; A.R.M.L. - discussion and writing; I.J.A., M.Y.N.J. and U.C.O. -- sample preparation and transcriptome; A.M.C.T. project idea and discussion. 

Transcriptomic data were submitted to GEO/NCBI with the accession number GSE145292.  Recombinant Amblyomin-X is patent granted (WO/2006/029492). 

The authors declare no competing interests.