Nymphs of O. brasiliensis were collected in the field from a site previously implicated in tick parasitism of travelers in Brazil (Dall’Agnol et al., 2019) and maintained unfed in the laboratory for ~2 months before dissection. Salivary glands were dissected from 10 nymphs and placed in RNA later before storing at −70°C. Total RNA was extracted using the RNeasy Protect Mini Kit (QIAGEN Group). Briefly, glands were suspended in 500 μl RLT buffer and disrupted by 10X passage through an 18G needle followed by 10X passage using a 24G needle. Residual genomic DNA was removed with DNase I digestion. Total RNA quantification was performed using the Qubit fluorimeter 2.0 (Life Technologies, Carlsbad, CA).

For library preparation, 1.0 μg purified total RNA was used with the TruSeq stranded mRNA sample preparation kit (Illumina, San Diego, CA). Poly-A mRNA was isolated, fragmented (for 3 min), and converted to double-stranded cDNA, adapters ligated, and PCR amplified for 12 cycles. Amplified bands were size selected from 450 to 1,200 bp. Bands were excised, purified, and sequenced using the Illumina MiSeq system (300 bp×300 bp). Raw sequence reads were submitted to GenBank under BioProject PRJNA719007 with small read archive accession number: SRR14139641.

Raw Illumina reads were quality trimmed (0.001 quality limit) and TruSeq adapters removed using CLC Genomics Workbench. Reads were imported as single or paired end reads. The paired end reads were merged to produce a merged dataset (Merged), while the single reads were used as unpaired (Single) and to produce a single-merged dataset (SM). Duplicates were also removed from these datasets to produce three duplicate removed datasets (Mddup, Sddup, and SMddup; Pienaar et al., 2021). These datasets were used to assemble the transcriptome using Trinity v2.4.0 and CLC Genomics Workbench v 20.0. Trinity was used with default parameter settings of kmer size 25. For CLC Genomics Workbench, kmer sizes were used in step sizes of 5 starting at 15 up to 60 and an additional assembly using kmer 64 (11 assemblies) with assembly parameters: mismatch cost-2, insertion cost-3, deletion cost-3, length fraction-0.9, similarity-0.9, minimum contig length-240, kmer size-variable, and bubble size-automatic. Given the different dataset structures used, a total of 72 assemblies were produced.

The mitochondrial genome was identified in the assemblies by BLASTN analysis (Altschul et al., 1993), using the previously published mitochondrial genome for O. brasiliensis (Burger et al., 2014). The mitochondrial genome was annotated using the MITOS server to identify tRNA genes (Bernt et al., 2013). Protein coding and rRNA genes were identified using BLAST analysis (Altschul et al., 1993). The translated COI, CYTB, ND1, ND2, and ND4 proteins were used for phylogenetic analysis as previously described (Mans et al., 2015, 2019, 2021).

Open reading frames (ORFs) were extracted using a Perl-script and chimeric and duplicate sequences removed by clustering at 90% identity using CD-HIT (Li and Godzik, 2006). The single dataset was mapped against the clustered ORFs using CLC Genomics Workbench and ORFs with TPM>1 (transcripts per million) were selected. BLASTX analysis against the ACARI database reduced the dataset further by selecting ORFs with E-values below 0.004 for further analysis. The transcriptome quality was measured for accuracy, completeness, contiguity, and chimerism using the Benchmarking Universal Single-Copy Orthologs (BUSCO; Simao et al., 2015). The final set of ORFs was submitted to GenBank under BioProject PRJNA719007.

To identify potential secretory peptides translated ORFs were submitted to SignalP (Petersen et al., 2011), while TMHMM and Phobius (Krogh et al., 2001; Kall et al., 2004) were used to identify membrane proteins. Potential housekeeping and secretory proteins were identified by BLASTP analysis against an ACARI database annotated using the KEGG database and GhostKOALA (Kanehisa et al., 2016a,b), the TSFAM database (Ribeiro and Mans, 2020), and an in-house annotation of secretory protein families (de Castro et al., 2016). To identify functional orthologs of proteins with experimentally verified functions, BLASTP analysis of secretory proteins was performed against the NCBI non-redundant database before phylogenetic analysis was performed to confirm clustering in functional clades with high bootstrap support. Protein families were aligned using MAFFT (Katoh and Standley, 2013), and maximum-likelihood analysis performed using IQ-Tree2 v 1.6.12 (Minh et al., 2020) with a standard 1,000,000 bootstraps.

A total of 58,490,632 paired reads were generated that yielded 2,285,451 merged reads (Merged) after quality trimming and merging and 10,505,566 single reads (Single) after quality trimming, resulting in a combined 12,791,017 reads (SM). Removal of duplicate reads resulted in 1,301,316 merged reads (Mddup), 6,841,083 single reads (Sddup), and combined 8,142,399 reads (SMddup). Assembly resulted in 100,397 contigs after clustering with CDHIT that was further reduced to 44,052 contigs with TPM>1. This resulted in 16,908 contigs with E-values >0.004 and 14,957 contigs after manual curation. BUSCO analysis indicated 95.0% completeness with 92.6% as single genes, 2.4% as duplicated, 2.7% fragmented, and 2.3% missing from a set of 1,066 conserved genes. This compares well with other tick transcriptomes sequenced thus far (Figure 1A). Comparison to two other soft tick salivary gland transcriptomes (O. rostratus and Ornithodoros turicata which belong to the Neotropic and Nearctic Pavlovskyella), for which protein sequence data are available in the public repositories, indicates that O. brasiliensis generally has longer contigs although it also has higher numbers of short contigs, which may indicate that some contigs may be truncated (Figure 1B). Reciprocal best hit analysis indicated that O. brasiliensis shares 7,441 orthologs in total with these transcriptomes (Figure 1C). This can be considered a minimum, since this is limited by the numbers of contigs submitted for the other transcriptomes, i.e., O. rostratus (n = 6,602) and O. turicata (n = 7,544). As such, the percentage of orthologs shared for each transcriptome is 89% of O. rostratus and 82% of O. turicata. These measures were taken to indicate a well-represented high quality transcriptome. BLASTP analysis of ACARI database using the transcriptome retrieved as highest number of hits, proteins from related soft ticks, such as O. rostratus, O. turicata, O. erraticus, and O. moubata (Figure 2). This may be expected but is also a good measure of transcriptome quality.

Previously, O. brasiliensis ticks shown to cause toxicoses and for which a mitochondrial genome was sequenced were collected at a site distant from the current study (Reck et al., 2013a; Burger et al., 2014). Since the collection sites were approximately 60 km apart, the question was raised regarding the genetic relationship of the ticks collected in this study and that of the mitochondrial genome previously sequenced for this species. The full-length mitochondrial genome was assembled in the transcriptome and was 99% identical to the previously sequenced genome and was deposited under accession number MW864544. Phylogenetic analysis indicated that both genomes cluster together with high support in a clade shared with O. rostratus in what has been described as the Neotropic Pavlovskyella clade (Figure 3). The sister-clade with high bootstrap support is that of the Afrotropic Ornithodoros sensu stricto group suggesting as previously indicated that these clades probably diverged during continental breakup of Gondwanaland ~127 MYA (Mans et al., 2019). This places the transcriptome within a phylogenetic context for comparative analysis and suggests that the transcriptomes of Afrotropic Ornithodoros and Neotropic Pavlovskyella should share extensive similarities with regard to orthologs and function. In addition, full-length 18S and 28S rRNA sequences were retrieved from the study and were deposited under GenBank accession numbers MW857182 and MW877711, respectively. These commonly used phylogenetic markers can therefore also be retrieved from transcriptome sequencing projects.

The transcriptome can be divided into housekeeping (n = 12,471), secretory (n = 368), and unknown categories (n = 2,118; Figure 4A). Housekeeping proteins are defined as all proteins not part of the secretory class that is involved in general cellular or organismal housekeeping functions, while secretory proteins are defined as those with secretory peptides, considered to be secreted into the feeding site during feeding, while unknown proteins refer to proteins with significant BLASTP hits in the database, but which has no annotation. These classifications have been used in all tick transcriptome studies dating back to Ribeiro et al. (2006). Housekeeping proteins account for ~83.3% of the transcriptome, secretory proteins for ~2.5%, and unknown proteins for ~14.1%. Reads mapped back to the transcriptome indicate that housekeeping proteins account for ~42% of the coverage, secretory proteins for ~53%, and unknowns for ~5%. This would suggest that secretory proteins are present at higher concentrations in the salivary glands relative to the housekeeping proteins. To corroborate this, the first 39 proteins with highest coverage are secretory proteins and comprise 43% of the total TPM coverage (80% of the secretory protein contribution; Figure 4B). The high abundance of secretory proteins was previously observed in soft tick salivary glands and may be explained by the fact that soft ticks synthesize secretory proteins and store them in large secretory granules that occupy most of the space in the salivary gland cells until secretion (Mans and Neitz, 2004; Mans et al., 2004). Abundance of secretory transcripts was also observed in conventional Sanger sequenced argasid salivary gland transcriptomes, as well as in those sequenced with next-generation sequencing technologies (Francischetti et al., 2008a,b; Mans et al., 2008a; Araujo et al., 2019).

The housekeeping proteins can be divided into various functional classes as classified by the KEGG database (Figure 5). These include proteins involved in metabolism, protein synthesis (transcription and translation), protein folding, sorting, degradation and excretion (FSDE), environmental sensing processes, such as signal transduction systems, cellular processes, such as cell growth, death, and motility, and organismal processes, such as the circulatory, developmental, digestive, endocrine, excretory, immune, nervous, and sensory systems. Signal transduction processes are notably enriched with regard to number of proteins. The most abundant housekeeping classes include those involved in protein synthesis (transcription and translation), FSDE (Figure 6A), that cumulatively account for 54% of TPM coverage for housekeeping proteins. This may be expected for an organ, such as salivary glands whose primary function is the synthesis, folding, and sorting of secretory proteins.

Some housekeeping functions with specific reference to tick biology may be highlighted. Of ~270 proteins identified in vertebrates to function as part of the protein secretory pathway (Gutierrez et al., 2020), 216 orthologs could be identified with confidence (Supplementary Table 1). The reason for many of those not identified may be due to multiple isoforms in the vertebrate secretory system which may be only represented by single proteins in ticks. This suggests that a large portion of the secretory system is present in the current transcriptome and is also conserved in ticks. This may be expected since the secretory pathway may be considered one of the Lineages of Life processes conserved in all Metazoa (Mans et al., 2016). Similarly, several orthologs of proteins previously implicated in neuronal control of salivary gland secretion in ticks (Šimo et al., 2014) were identified (Supplementary Table 2). This included hormones, such as bursicon, crustacean hyperglycemic hormone, eclosion hormone, elevenin, insulin, and orcokinin, as well as hormone receptors, such as allatostatin receptor, calcitonin gene-related peptide type 1 receptor, corazonin receptor 1, dopamine receptor 1, dopamine receptor 2-like, dopamine D2-like receptor, elevenin receptor 2, insulin-like growth factor 1 receptor, tachykinin-like peptides receptor, and periviscerokinin/Cap2b receptor.

It was previously indicated that the heme biosynthesis pathway is incomplete in ticks with ixodids lacking hemA-hemE, but possess hemF-hemH (Braz et al., 1999; Mans et al., 2016; Perner et al., 2016). Argasids possessed hemB and hemF-hemH (Mans et al., 2016). BLASTP analysis using hemA-hemH for Metaseiulus occidentalis indicated orthologs for hemB, hemF, and hemG, consistent with previous observations for argasids.

Secretory protein families identified in the transcriptome included the majority of protein families generally found in tick salivary glands (Figures 6B, 7). The most abundant protein families in terms of sequence coverage (TPM) were the lipocalin, basic tail secretory family (BTSP), 7 disulfide bond domain (7DB), Kunitz-BPTI domains, 7 cysteine domain, and reprolysin, comprising ~98% of the total (Figure 6B). The lipocalin, basic tail, reprolysin, and Kunitz-BPTI families were also the most abundant in terms of the number of family members (Figure 7). This has previously been observed in soft tick salivary gland transcriptomes (Francischetti et al., 2008a,b; Mans et al., 2008a; Araujo et al., 2019). In addition, a number of peptides with secretory signals but no BLASTP hits were also identified as unknown proteins.

Approximately 120 protein functions involved at the tick-host interface has been experimentally validated (Mans et al., 2019). BLASTP analysis using these proteins found 73 potential orthologs for 13 functions (Table 1). This includes the 5' nucleotidase family member apyrase that inhibits platelet aggregation by hydrolyzing ADP (Ribeiro et al., 1991; Mans et al., 1998; Stutzer et al., 2009).

For the lipocalin family, a number of orthologs were found for proteins with known function, including leukotriene B₄ (LTB₄) scavengers (Mans and Ribeiro, 2008a). The LTB₄ scavengers did not have the motifs conserved for complement C5 inhibition or thromboxane A₂ scavenging and probably do not possess these functions (Nunn et al., 2005; Mans and Ribeiro, 2008a). Orthologs were also found for the biogenic amine scavengers and these orthologs possessed the biogenic amine-binding motif of the lower-binding site, while some also possessed the conserved amino acid residues involved in the upper-binding site of TSGP1 (Mans et al., 2008b). These orthologs therefore likely bind both histamine and serotonin. Orthologs for leukotriene C₄ scavengers were also found (Mans and Ribeiro, 2008b).

In addition, an ortholog for the lipocalin (savicalin) was found. Savicalin was implicated in antimicrobial activity (Cheng et al., 2010). Other potential secretory antimicrobial orthologs present are defensins, microplusins, and trypsin inhibitor-like domains (Nakajima et al., 2001; Fogaça et al., 2004, 2006; Sasaki et al., 2008).

Four orthologs for adrenomedullin (vasodilation) were found in the transcriptome. Adrenomedullin was possibly acquired by soft ticks from the genus Ornithodoros by horizontal gene transfer from mammals and has thus far been found in O. coriaceus, O. moubata, O. parkeri, and O. rostratus (Iwanaga et al., 2014; Araujo et al., 2019). It would therefore seem as if the horizontal gene transfer occurred in the ancestral lineage to the Pavlovskyella and Ornithodoros (Figure 3). The presence of adrenomedullin in O. brasiliensis is therefore not surprising.

For members of the Kunitz-BPTI family, three orthologs of savignygrin were found. Single domain orthologs to the savignygrins were found that possessed the integrin RGD motif on the substrate binding presenting loop of the BPTI fold (Mans et al., 2002a). These inhibitors target the fibrinogen receptor integrin α_IIbβ₃ and inhibit platelet aggregation. Orthologs have also been found in A. monolakensis suggesting that these inhibitors were evolved in the ancestral argasid lineage (Mans et al., 2008c). A double-domain Kunitz-BPTI protein (Obras12157) with the same integrin RGD recognition motif located in the second Kunitz domain was also found. Such double Kunitz domain proteins were also observed in the transcriptomes of O. coriaceus and O. parkeri (Francischetti et al., 2008a,b). To date, no experimental evidence exists that these target integrin α_IIbβ₃, but it is likely that they target an integrin. The savignygrins are the Kunitz-BPTI proteins with the highest coverage for this family, as may be expected for an inhibitor that targets highly abundant platelet receptors in the host (Mans, 2019). Of interest is that no orthologs were found for the thrombin inhibitors like monobin, ornithodorin, or savignin (van de Locht et al., 1996; Nienaber et al., 1999; Mans et al., 2002b, 2008c), or the fXa inhibitors (Waxman et al., 1990; Gaspar et al., 1996; Joubert et al., 1998). This was also observed for Nearctic tick species, such as O. coriaceus and O. parkeri (Francischetti et al., 2008a,b). It was suggested that fXa inhibitors evolved exclusively within the genus Ornithodoros while thrombin inhibitors evolved in the ancestral argasid lineage (Mans et al., 2008c). Absence of these inhibitors in other Ornithodorinae lineages may indicate gene losses probably associated with host preferences. Possible orthologs to longistatin may be involved in blood coagulation modulation by activating plasminogen and in the immune system by targeting the receptor for advanced glycation end products (Anisuzzaman et al., 2011, 2014). Another ortholog associated with fibrinolysis is carboxypeptidase inhibitor that targets plasma carboxypeptidase B (thrombin-activatable fibrinolysis inhibitor) leading to fibrinolysis (Arolas et al., 2005). As such, a number of inhibitors that targets the blood coagulation cascade are present.

Orthologs of cystatin-2 from O. moubata were also detected. These cystatins inhibit cathepsin L and S allowing for modulation of the inflammatory responses in the host (Grunclová et al., 2006; Salát et al., 2010). Orthologs of serpins were also found that may inhibit elastase and cathepsin G acting as immunomodulators and platelet aggregation inhibitors (Prevot et al., 2006; Chmelar et al., 2011). It should be noted that the orthologs to which functions could be ascribed comprise ~19% of all the secretory proteins found in the transcriptome, suggesting that numerous undescribed functions still exist for this tick species.

The toxicoses and bite of O. brasiliensis are accompanied by pruritis, local edema and erythema, pain, and blisters, while histopathology of the feeding site indicates extensive subcutaneous edema and hemorrhage (Reck et al., 2013a, 2014). Possibly, salivary proteins that may contribute to this clinical outcome would include the clotting and platelet aggregation inhibitors that would induce hemorrhage. On the other hand, the LTB₄ and LTC₄ scavengers may inhibit edema and erythema (Mans and Ribeiro, 2008b). Additional contributors to edema and hemorrhage may be the secretory proteolytic enzymes found in the transcriptome that includes a cathepsin, two serine proteases, and 38 metalloproteases (astacin, gluzincin, and reprolysin). These may also impact in wound-healing processes leading to prolonged recovery times.

The salivary gland transcriptome of O. brasiliensis shows a number of evolutionary features previously found for argasids. This includes the presence of the major salivary gland secretory protein families conserved in all ticks (Mans et al., 2008a), as depicted in Figure 7. A number of orthologs for functions thus far conserved in argasids were also identified and include biogenic amine, leukotriene B₄, and leukotriene C₄ scavengers, as well as apyrase, savignygrin, and defensins (Mans and Ribeiro, 2008b). These functions are thus far conserved in all argasid species that feed on blood (Mans et al., 2016) and underscore their important role in blood-feeding. The only argasid salivary gland transcriptome sequenced to date that did not show these conserved features of protein families and conserved functions, were the salivary gland transcriptome for adult A. delacruzi (Ribeiro et al., 2012), which lacked all conserved features. This tick does not feed on a host in the adult phase and it is likely that this difference is due to differential expression in different life stages, and that these conserved features will be found in its blood-feeding stages. Since Antricola groups well within the Ornithodorinae, it may be expected that their lack of conserved protein families and functions in non-feeding phases is a derived trait. Conversely, the conserved protein families and functions as observed in O. brasiliensis and other argasids, indicate evolution of these traits in the last common ancestor to the Argasidae, or even the Ixodida (Mans et al., 2016).