Skip to main content
  • Research article
  • Open access
  • Published:

Unraveling the role of male reproductive tract and haemolymph in cantharidin-exuding Lydus trimaculatus and Mylabris variabilis (Coleoptera: Meloidae): a comparative transcriptomics approach



Meloidae (blister beetles) are known to synthetize cantharidin (CA), a toxic and defensive terpene mainly stored in male accessory glands (MAG) and emitted outward through reflex-bleeding. Recent progresses in understanding CA biosynthesis and production organ(s) in Meloidae have been made, but the way in which self-protection is achieved from the hazardous accumulation and release of CA in blister beetles has been experimentally neglected. To provide hints on this pending question, a comparative de novo assembly transcriptomic approach was performed by targeting two tissues where CA is largely accumulated and regularly circulates in Meloidae: the male reproductive tract (MRT) and the haemolymph. Differential gene expression profiles in these tissues were examined in two blister beetle species, Lydus trimaculatus (Fabricius, 1775) (tribe Lyttini) and Mylabris variabilis (Pallas, 1781) (tribe Mylabrini). Upregulated transcripts were compared between the two species to identify conserved genes possibly involved in CA detoxification and transport.


Based on our results, we hypothesize that, to avoid auto-intoxication, ABC, MFS or other solute transporters might sequester purported glycosylated CA precursors into MAG, and lipocalins could bind CA and mitigate its reactivity when released into the haemolymph during the autohaemorrhaging response. We also found an over-representation in haemolymph of protein-domains related to coagulation and integument repairing mechanisms that likely reflects the need to limit fluid loss during reflex-bleeding.


The de novo assembled transcriptomes of L. trimaculatus and M. variabilis here provided represent valuable genetic resources to further explore the mechanisms employed to cope with toxicity of CA in blister beetle tissues. These, if revealed, might help conceiving safe and effective drug-delivery approaches to enhance the use of CA in medicine.


Insects use a wide array of defensive toxic compounds which may be produced autogenously by de novo synthesis or acquired by dietary sequestration of secondary compounds derived from plants, ingested prey, parental transfer, or endosymbionts. The strategy for using these toxins greatly differs among insects. ‘Reflex bleeding’ or ‘autohaemorrhaging’ [1] is the release of toxin-containing haemolymph from integumental ruptures when individuals are threatened or exposed to direct physical attack. Among all defensive strategies, this is probably the most physiologically costly and extreme to serve a harmful compound to a predator [2,3,4]. Among insects, Plecoptera, Orthoptera and Hemiptera include some taxa exhibiting reflex-bleeding defensive mechanisms [5]. In contrast, autohaemorrhaging is common in several species from different families of Coleoptera (e.g. Erotylidae, Lampyridae, Coccinellidae), but represent a typical trait of blister beetles (Meloidae) [6].

Meloidae is a widespread family including almost 3000 species [6,7,8] well-known for secreting cantharidin (C10H12O4; hereafter named CA), a toxic terpene. Blister beetles, if disturbed, naturally exude CA in yellowish oily haemolymph droplets from leg and antennal joints to defend themselves from predators. The concentration of the secreted CA vary from 0.03 to 0.79 mg for each gram of exuded haemolymph [9, 10], but largely depends on species, environmental and physiological conditions. Males produce more CA than females [9,10,11,12,13,14,15] and accumulate large amounts of this terpene in their reproductive organs and, particularly, in male accessory glands (MAG) [13, 16,17,18,19]. During mating, males transfer large quantities of CA to females [6, 20,21,22,23] which, in turn, use the received compound for protecting eggs from potential predators or parasites [13, 16, 21, 24].

CA is well-known in popular pharmacology for its traditional use as a sexual stimulant or antiphlogistic agent [6, 25,26,27]. However, CA is also a potent blistering compound causing many adverse effects after ingestion, such as severe damages to gastrointestinal, kidney and urinary tracts [6, 27]. The large use of CA in modern medicine is hindered by its extreme toxicity that renders its employment limited to the topical treatment of warts under strict legislative regulations [27, 28].

Recently, CA and its derivatives have regained popularity as alternative compounds for anti-cancer treatments [29,30,31,32,33]. Such a renewed interest has led to a growing literature devoted to reveal the molecular basis of CA de novo biosynthesis in blister beetles, focusing on species belonging to the (species-rich) Meloinae subfamily. Both gene-expression analyses of key-enzymes in different organs of Epicauta sibirica Pallas, 1773 (tribe Epicautini; all reported as chinensis Laporte, 1849) [18, 19, 34, 35] and a de novo transcriptomic approach comparing the relative abundance of transcripts in males vs. females in Hycleus cichorii (Linnaeus, 1767) (tribe Mylabrini; reported as Mylabris) [36] were coherent with previous studies indicating that: (i) CA may be synthetized via the mevalonate (MVA) pathway [37], (ii) farnesol may act as an intermediate [38,39,40], and (iii) a juvenile hormone (JH) metabolite could be involved in one of the latest biosynthetic steps [35, 37, 41].

However, there is still a gap of knowledge on some crucial steps of the CA biosynthetic pathway, as well as the specific site of its in vivo production. At first, the third pair of MAG was designated as the production site [16, 42], but recently, both high CA content and transcriptional level of the 3-hydroxy-3-methylglutary-CoA reductase (HMGCR, an essential enzyme of the MVA pathway) observed in the fat body of E. sibirica suggested that CA biosynthesis may more likely occur in this tissue [19].

McCormick et al. were the first to suggest that MAG represented the storage organ for CA and not the site of its biosynthesis [21, 41]: these authors hypothesized that the newly synthetized CA could first proceed to the haemolymph, from which it would have been removed by a ‘cantharidin kidney’ and stored in MAG [21]. Whatever the mechanism, a not-yet thoroughly addressed question regards the safe modes of CA circulation and storage. It is in fact still unclear how self-protection from the toxicity of a high CA content could be prompted in male reproductive organs, in which mechanisms able to inactivate or mild the action of the (incoming and accumulating) CA might be operating. Furthermore, still unidentified mechanisms should allow CA to be transported through the haemolymph from the fat body (or another production sites) to the storing organs, and then discharged outside the body without causing self-injuries.

A third, more general (and evolutionary) question concerns the identification of proteins in the haemolymph repertoire that should support such a metabolically expensive defensive strategy. In fact, yet highly effective, a haemolymph-mediated chemical defence requires a rapid renewal of water [43] and components of the immune system (e.g. alkaloids, haemocytes or antimicrobial peptides) after bleeding. Moreover, both efficient haemolymph-coagulating and integument repairing mechanisms are expected, as these should be essential requirements of autohaemorrhaging insects to limit haemolymph loss and guarantee their survival.

With this study, we attempted for the first time to shed light on modes of mobilization, storage, and deactivation of CA in blister beetles which shall prevent its adverse effects on organs/tissues where CA is stored and/or regularly circulates. To do so, we used a de novo assembly transcriptome approach and a differential gene expression analysis targeting the reproductive organs and exuded haemolymph (gathered during autohaemorrhaging) of males of two Meloinae species, i.e. Lydus trimaculatus (LT) [Fabricius, 1775 (tribe Lyttini)] and Mylabris variabilis (MV) [Pallas, 1781 (tribe Mylabrini)]. The choice of these two species was dictated by their: 1) proved ability to produce high titres of CA in natural conditions [10]; 2) relatively high abundance on field [10], crucial to acquire the amount of haemolymph necessary to gain reasonable RNA yields, and 3) ascription to two different tribes of Meloinae owning the proper level of phylogenetic divergence [6, 8] to observe a co-occurrence of abundantly expressed candidate detoxification genes. Besides, our approach allowed examining the expression levels of the thus-far-identified genes related to CA biosynthesis in male reproductive tracts (MRT) of LT and MV. We also provided a catalogue of over-represented protein families in haemolymph involved in the costly physiology of the autohaemorrhaging strategy.


De novo assembly and annotation of LT and MV transcriptomes

Our final dataset consisted of 252,190,732 and 301,963,686 (50–150 bp) quality-filtered reads for LT and MV, respectively (Table 1). Both de novo assemblies received a high score of completeness (BUSCO ‘C’ values = ~ 97%) and produced for LT a total of 190,214 assembled transcripts, ranging from 401 to 31,991 bp, and corresponding to 111,614 Trinity unigenes, and for MV a total of 165,191 assembled transcripts, ranging from 401 to 26,915 bp, and corresponding to 94,199 Trinity unigenes (Table 1). We found a higher percentage of duplicated BUSCOs (D) in LT than in MV (Table 1). This can be considered a possible artifact originated by the high number of allelic variants (= high heterozygosity) in LT populations, being interpreted as multiple unigenes/isoforms. Since this did not affect our main conclusions, we decided to not de-replicate the obtained transcriptomes and maintain the whole information. By selecting the longest transcript isoform for each Trinity unigene cluster and using the Annocript pipeline (see Methods), we observed that 47.7% (LT) and 33.4% (MV) could be annotated with BLAST searches against UniRef [mostly from Tribolium castaneum Herbst, 1787; Coleoptera, Tenebrionidae], whereas 30.0% (LT) and 20.8% (MV) were annotated against Swiss-Prot [mostly from Drosophila melanogaster Fallén, 1823; Diptera, Drosophilidae]. Finally, 30.2% (LT) and 21.3% (MV) of the longest transcripts matched with functional protein domains in InterProScan with at least one match in one of the following repositories: SignalP, TMHMM, Pfam, SMART, Tigr and ProfileScan.

Table 1 De novo transcriptome assembly statistics, quality control and annotation for LT and MV

Detailed annotation information for individual LT and MV longest transcripts is available in Additional file 1. Classification and relative abundance of top 15 conserved Protein families (Pfam) and Gene Ontology (GO) terms are reported in Fig. 1. Residue motifs mediating protein-protein interactions and signalling, such as leucin rich (Pf13855) and ankyrin (Pf12796) repeats, as well as immunoglobulin-related domains (Pf07679, Pf13927) were among the most abundant Pfam in transcriptomes of both species (Fig. 1). ‘Protein kinase’ (Pf00069), ‘Cytochrome P450’ (Cyp450) (pf00067), ‘Reverse Transcriptase’ (pf00078) and ‘Trypsin’ (Pf00089) were also well-represented in both species (Fig. 1A), as well as protein families involved in movement of small solutes, such as ‘Major Facilitator Superfamily’ (Pf07690) and ‘Sugar (and other) transporters’ (Pf00083) (Fig. 1A). Overall, the relative abundances of GO were also coherent between the two analysed species. Among the most significant ‘GO-biological processes’, three DNA-related terms (‘transcription’, ‘regulation of transcription’, ‘DNA integration’) were highly frequent, but interestingly, also processes associated to transport (e.g. ‘transmembrane transport’ and ‘intracellular protein transport’) were highly represented (Fig. 1B). The most represented locations of gene products in cellular structures (‘GO-cellular components’) were, in order, ‘integral to membrane’, ‘nucleus’ and ‘cytoplasm’ (Fig. 1C). Finally, among ‘GO-molecular functions’, ‘ATP-binding’, ‘nucleic acid binding’, ‘DNA binding’ were the most recurrent, but, noteworthy, a number of other GO binding-related terms - such as metal-ion binding (zinc, calcium, iron) - were also rather abundant (Fig. 1D).

Fig. 1
figure 1

Classification and relative abundance of top 15 protein families (Pfam) (A) and Gene Onthology (GO) terms in each category (B, C, D) in LT and MV transcriptomes. Pfam and GO differentially expressed between species are marked with * (p < 0.05)

Analysis of differentially expressed genes (DEGs)

Overall, the expression profiles observed for both targeted tissues in LT and MV were comparable in terms of relative abundance of differentially expressed transcripts (Fig. 2). Following a conservative approach (logFC > 2, FDR < 0.001), 1800 (LT) and 568 (MV) transcripts appeared upregulated in MRT when compared to the whole body, whereas 1275 (LT) and 1471 (MV) transcripts showed an enhanced expression in haemolymph (Fig. 2). Among the upregulated transcripts in these two tissues, ~ 15% did not display a significant similarity with known proteins in the non-redundant (nr) database after BLAST search. Comparison of DEGs between the two targeted tissues showed 5005 (LT) and 6139 (MV) transcripts with enhanced expression in MRT with respect to haemolymph, whereas 2186 (LT) and 2527 (MV) were upregulated in the haemolymph respect to MRT (Fig. 2). The complete list (and annotation) of DEGs (i.e. upregulated) in MRT and haemolymph of each species is reported in Additional file 2 (a list of downregulated genes in both tissues is available upon request).

Fig. 2
figure 2

Tissue-specific gene expression of MRT and haemolymph in LT and MV. Volcano plots display the relative abundance levels of upregulated transcripts in MRT and haemolymph compared to the whole body. The x-axis represents the log2 of the expression ratio (FC = fold change) for each transcript (tissue specific logCPM: whole body logCPM, where CPM stands for Counts Per Million reads); the y-axis represents the log10 of the p-value corrected for the false discovery rate (FDR). Red dots represent differentially expressed transcripts with logFDR < 0.001 and at least 2-fold difference in logCPM; black dots if logFDR > 0.001. Negative logFC values indicate transcripts enhanced in the target tissue subsets, while positive values indicate transcripts upregulated in the whole body

Significantly over-represented GO terms in DEGs (when compared to the whole assemblies) are reported in Additional file 3. Among the most abundant ‘GO-biological processes’ already observed in the two whole transcriptomes, transmembrane and intra-cellular protein transportations were also specifically enriched in MRT and haemolymph of both species (Additional file 3). Enriched-GO shared between species in MRT and haemolymph are reported in Tables 2 and 3, respectively.

Table 2 Most significant (p < 0.05) GO terms enriched in MRT of both LT and MV. ALL = count number of all GO retrieved in transcriptome; UP = count number of enriched GO
Table 3 Most significant (p < 0.05) GO terms enriched in haemolymph of both LT and MV. ALL = count number of all GO retrieved in transcriptome; UP = count number of enriched GO

In MRT, dynein, cilium- and microtubule-related GO terms (likely associated with cytoskeletal dynamics in spermatogenesis) were the most over-represented, whereas others were related to ubiquitination, transcription regulation and phosphorylation (i.e. ‘serine/threonine kinase’) (Table 2). In haemolymph, most of shared enriched GO terms pertained to DNA replication, cytoskeleton, cell adhesion, actin binding (probably related to clotting), signal transduction and vesicle mediated-transport, but responses to oxidative stress and activities of G-protein coupled receptors were also recorded (Table 3). MRT- and haemolymph-specific Pfam upregulated in each species (adj. p < 0.05) are reported in Additional file 4. The subsets of Pfam (adj. p < 0.05) shared between LT and MV are reported in Table 4 (MRT) and Table 5 (haemolymph), respectively.

Table 4 Pfam shared and enriched in MRT of both LT and MV (adj p < 0.05). ALL = count number of all Pfam retrieved in transcriptome; UP = count number of enriched Pfam
Table 5 Pfam shared and enriched in haemolymph of both LT and MV (adj p < 0.05). ALL = count number of all Pfam retrieved in transcriptome; UP = count number of enriched Pfam

In MRT most of shared Pfam were related to motor functions of cilia and flagella, microtubule dimerization and transport of organelles and vesicles along microtubules (i.e. dynein-related subfamilies, e.g. pf12774 and pf03028; Table 4). Other shared and highly enriched domains were possibly involved in structural and/or regulation processes (i.e., pf13868; pf00412; pf15921; Table 4). Noteworthily, some significantly enriched domains in MRT of LT (but also detected in MV, yet with p > 0.05) were related to protein phosphorylation (pf00069) or, possibly, detoxification, such as glutathione synthases (pf03917) and zinc carboxypeptidases (pf00246) (Additional file 4: Table S11-S12). Additional species-specific protein families were linked to other processes, such as aminopeptidases (pf00883, pf02789) in LT or carboxylases (pf05090) and hydrolases (pf00723) in MV (Additional file 4: Table S11-S12). In haemolymph, most of shared enriched-Pfam were related to actin cross-linking proteins (e.g. pf00630, pf06268), vesicle transport (e.g. pf00191) or cell/membrane adhesion (e.g. pf00052) (Table 5). Other shared domains were related to oxygen transport (hemocyanin; pf03723), DNA replication processes (Minichromosome maintenance; pf00493, pf14551, pf17207), RNA binding (Pumilio family; pf00806) and lipid-binding (e.g. StAR-related lipid-transfer; pf01852) (Table 5). Other, but more species-specific domains were still related to the above-mentioned functions (Additional file 4). Interestingly, a subset of Pfam were clearly related to innate immunity responses (e.g. Spätzle: pf16077 and CUB: pf00431) with some domains more specifically involved in coagulation (haemolectins/mucins; e.g. C8: pf08742; von Willebrand factor: pf00094; Mucin2_WxxW: pf13330) and, possibly, in the formation of the viscous (oily) haemolymph-exudate (Table 5; Additional file 4: Table S13-S14).

Genes putatively involved in CA biosynthesis

We searched in both transcriptomes for unigenes related to the MVA pathway and trans-farnesol branch (“terpenoid backbone biosynthesis”, KEGG map 009000), as well as for enzymes involved in the JH biosynthesis (“insect hormone biosynthesis”, KEGG map 00981). In these pathways, we selected genes previously reported as related with CA content in blister beetles (i.e. in Epicauta and/or Hycleus, whose sequence data were available); further selection was driven by protein sequence homology with the ones from Epicauta, which allowed recognising six transcripts with high similarity (70–93%). Those pathways are summarized in Fig. 3A, in which the six genes selected are highlighted. Intensity of expression, as revealed by FPKM counts in transcriptome analyses, is reported in Fig. 3B (a complete list of transcripts emerging from BLAST searches on LT and MV transcriptomes can be found in Additional file 5: Table S15). MVA enzymes expressed in LT and MV transcriptomes included acetyl-CoA C-acetyltransferase (atoB;, hydroxyl-methylglutaryl-CoA synthase (HMGCS;, hydroxyl-methylglutaryl-CoA reductase (HMGCR;, mevalonate kinase (mvaK1;, phosphomevalonate kinase (mvaK2; and diphosphomevalonate decarboxylase (MVD; Transcripts encoding these enzymes all showed high identity (~ 60–80%) in both species with Tribolium protein putative orthologs (except for mvaK1 showing a 44% of similarity, and, thus, not further considered), but none was significantly upregulated in MRT or haemolymph in respect to the whole body. Among these genes, HMGCR - a rate-limiting enzyme in MVA pathway and crucial for CA synthesis [34] - was downregulated in MRT (Fig. 3B). We also searched for transcripts related to the farnesyl diphosphate synthase (FDS;,, a key enzyme in farnesol metabolism, but also intervening in the formation of the geranyl diphosphate, a precursor of sesquiterpenoids, triterpenoids and monoterpenoids. This gene, moderately abundant in the whole body of both species (̴30–45 FPKM), was significantly downregulated in MRT and haemolymph of LT and MV (Fig. 3B). Among genes belonging to the branched chain of trans-farnesol (a potential precursor of CA), we identified NADP+-dependent farnesol dehydrogenase (FOHSDR:, protein-S-isoprenylcysteine O-methyltransferase (ICMT:, STE24 endopeptidase (STE24: and protein farnesyltransferase subunit beta (FNTB: All transcripts putatively encoding these genes were expressed but downregulated in MRT or haemolymph respect to whole body (Additional file 5: Table S15; Fig. 3B). Since a metabolite of JH could constitute a precursor of CA, we also investigated genes involved in JH biosynthetic pathway. Transcriptomic data showed that transcripts of two key enzymes responsible for the de novo synthesis of JH, i.e. methyl farnesoate epoxidase (MFE:; and JH acid O-methyltransferase JH-III synthase (JHAMT:, were expressed at very low levels in both species (Fig. 3B). On the contrary, the juvenile hormone epoxide hydrolase (JHEH), involved in the degradation of JH, was abundant in whole bodies (̴ 48–97 FPKM), though significantly downregulated in MRT of both species (Fig. 3B). Moreover, searching for members of the Cyp450 family, we identified unigenes showing 75% similarity with CYP4BM1, whose high expression level was demonstrated to be related to CA production [19]. These transcripts were significantly downregulated in MRT of both species as compared to the whole body (Additional file 5: Table S15) (as well as all transcripts showing a similarity with CYP4C7, a sesquiterpenoid omega-hydroxylase degrading JH III [44]).

Fig. 3
figure 3

Pathway and expression levels (Heat-Map and RT-qPCR) of genes putatively involved in CA biosynthesis in LT and MV. A) The diagram shows a recap of the pathways involved in the synthesis and degradation of JH, highlighting the enzymes whose transcripts were here selected and examined. Green ellipses represent the enzymes of the upstream mevalonate and terpenoid backbone pathways, the red ones indicate the enzymes specific for the downstream JH pathway. B) Colours in the Heat-Map were assigned according to percentiles, with yellow referring to the lowest FPKM values (below 30th percentile) and dark blue to the FPKM highest values (above 95th percentile). Identity percentages of tBLASTn using Epicauta protein sequences (as in GenBank Acc. N.) as a query are reported in the last column. B =  body, MRT =  male reproductive tract, H =  haemolymph. C) The graphs show RT-qPCR results expressed as fold change (2-ΔΔCt) for three enzymes of the JH pathway, i.e. JHAMT, JHEH and MFE in LT and MV, respectively

Finally, we did not detect specific transcripts associated to genes of the sesquiterpenoid, triterpenoid or monoterpenoid biosynthetic pathways. Also, we did not find transcripts possibly encoding enzymes catalysing cyclization reactions (as would be required for the transformation of a linear sesquiterpenoid into the tricyclo-decane structure of CA), such as the trans-isoprenyl diphosphate synthase involved in cyclization of 8-oxogeranial to build the iridoid ring scaffold in flea beetles (Coleoptera: Chrysomelidae: Alticini) [45]. Rather, we identified transcripts referable to two key enzymes involved in the formation of terpenoid-quinone biosynthesis, such as the geranylgeranyl pyrophosphate synthase ( and the 1,4-dihydroxy-2-naphthoate octaprenyltransferase (, which, however, were not differentially expressed in MRT or haemolymph of both species as well.

Identification of candidate genes related to CA detoxification and transport

The most conservative approach (logFC> 2, FDR < 0.001) identified 157 and 70 orthologous transcripts upregulated in both species in MRT and haemolymph, respectively. By relaxing the threshold (i.e., FDR < 0.01) the number of orthologous transcripts between the two species raised to 482 (MRT) and 283 (haemolymph), respectively. By applying the less conservative threshold, upregulated orthologous unigenes potentially involved in CA transport and/or detoxification were identified (Fig. 4A; Additional file 5: Table S16). Among them, two members belonging to ‘Major Facilitator Superfamily’, four to ‘ABC transporters’ and three to proteins involved in mobilization of small solutes (i.e. three sugar transporters, one sodium:solute symporter and one choline transporter) were identified in MRT, whereas three lipid-binding proteins were detected in the haemolymph (Fig. 4A). No unigenes coding for members of the most common families of detoxification enzymes (e.g. CypP450s, glutathione S-transferases) were identified.

Fig. 4
figure 4

Expression levels (Heat-Map and RT-qPCR) of genes putatively involved in sequestration mechanisms related to CA detoxification in LT and MV. A) Colours in the Heat-Map were assigned according to percentiles, with yellow referring to the lowest FPKM values (below 20th percentile) and dark blue to the FPKM highest values (above 95th percentile). Transcripts showing identities with those involved in salicin sequestration in Chrysomela populi (see Discussion) are marked with (*). B) The graphs show RT-qPCR results expressed as fold change (2-ΔΔCt) for transcripts showing similarities to the ones involved in salicin sequestration in Chrysomela populi, i.e. PCFT, SWEET1, SLC5A12 and SLC44A1 in LT and MV, respectively. C) The graphs show RT-qPCR results expressed as fold change (2-ΔΔCt) for other transcripts possibly involved in CA transport, i.e. PLT4, l(2)03659, ABCF2, MDR49 and ABCB6 in LT and MV, respectively

Reverse transcription-quantitative PCR (RT-qPCR) validation

Abundance of transcripts emerging from transcriptomic analyses has been confirmed by RT-qPCR by examining 12 transcripts differentially modulated in MRT and total body. We first examined transcripts coding for both JHAMT and MFE, two key enzymes of the latest steps of JH biosynthetic pathway, and for JHEH, a JH-degrading enzyme putatively involved in the production of a CA precursor. All three unigenes did not show any upregulation in MRT (Fig. 3C). Then, we analysed 4 transcripts coding for proteins whose possible homologs in Chrysomela populi Linnaeus, 1758 (Coleoptera, Chrysomelidae) are involved in transport of a plant-derived glycoside (salicin) from haemolymph to specific defensive glands (see Discussion): (i) proton-coupled folate transporter (PCFT), (ii) sugar transporter (SWEET1), (iii) sodium-coupled monocarboxylate transporter 2 (SLC5A12) and (iv) choline transporter-like protein 1 (SLC44A1) (Fig. 4A). Two of them, PCFT and SLC5A12 showed a strong differential expression in MRT compared to whole body (Fig. 4B). Finally, 5 transcripts encoding other putative transporters, chosen for ligand similarity, were analysed for their possible involvement in CA sequestration: i) probable polyol transporter 4 (PLT4), ii) probable multidrug resistance-associated protein lethal (2)03659 [l (2)03659], iii) ATP-binding cassette sub-family F member 2 (ABCF2), iv) multidrug resistance protein homolog 49 (MDR49), v) ATP-binding cassette sub-family B member 6 mitochondrial (ABCB6) (Fig. 4A). All of them exhibited a higher expression in MRT compared to whole body, particularly in MV (Fig. 4C).

The validation analysis was performed using independent, triplicate RNA samples extracted from whole bodies and MRT. Correlation analyses indicate statistically significant (R2 = 0.5878; p < 0.0001, for LT; R2 = 0.2518; p < 0.0286, for MV) linear relationships between RNA-seq and RT-qPCR results in both species (Fig. 5). This result strongly supports the transcriptional abundance profiles revealed by RNA-seq analyses, in terms of both absolute transcript abundance and relative fold-change differences of transcripts among tissue and samples.

Fig. 5
figure 5

RT-qPCR validation. Correlation between transcriptional abundance of 11 genes (in LT) and 10 genes (in MV) in both whole body and dissected MRT samples, as revealed by RT-qPCR and RNA-seq. Level of abundance is defined as the ratio between each sample value over the group median (mean FPKM and mean ΔCt for RNA-seq and RT-qPCR data, respectively) in both RT-qPCR and RNA-seq approaches. For both techniques, statistical evaluation throughout Pearson test was performed and relative results are reported in the figure insets


The present study produced a de novo assembly of male transcriptomes of LT and MV, two blister beetle species (Meloidae) belonging to the Meloinae subfamily. A high-throughput characterization of the upregulated portion of transcriptome in MRT and haemolymph of CA-exuding blister beetles is here provided for the first time. The two target tissues are essential in understanding the biological defence mechanisms in Meloidae due to their close interaction with CA; in fact, the MRT corresponds to a specialized body compartment where CA is accumulated in large quantities (and from which it is released upon mating), whereas haemolymph represents a system of internal vehiculation and outward emission of CA under stressful conditions. The most relevant aspects emerging from our comparative analysis are discussed below.

High quality of LT and MV de novo transcriptome assembly and conservation of gene expression patterns

A high-quality de novo transcriptome assembly, supported by the high percentage of (BUSCO) completeness (̴ 97%) in terms of the expected gene content, was produced for LT and MV (Table 1). The overall size of the assembled transcriptomes (̴ 200 Mb), as well as the number of retrieved transcripts (̴180,000 on average), were similar in the two target Meloinae species (Table 1). These results likely mirror a comparable amplitude in the two species of the gene repertoire active in males, as it is expected in (relatively) phylogenetically close taxa having analogous physiological, ecological, and behavioural traits [6, 7]. It is also worth noticing that the most frequent protein domains and functional categories to whom the annotated transcripts were ascribed were nearly the same for both beetles (Fig. 1) and partially overlapped with those recovered in Hycleus [36]. These were basically related to nucleic acid synthesis and regulation, signal transduction, cytoskeleton-mediated transport, and metabolite transfer, and might be coherent with the activation of metabolic processes of sexually mature male insects. Such a convergence may also depend on the fact that males were in both cases collected from field after their emerging period, maintained in laboratory at steady rearing conditions and kept nourished on their natural host plants. The observed resemblance, however, remains partly speculative, since, as in other non-model insects [46], nearly 50–65% of recovered transcripts were unannotated. This however likely constitutes a valuable portion of transcriptomes to disclose possible underlying differences, and, above all, to identify proteins conferring peculiar adaptations to blister beetles. Remarkably, we also observed that the number of transcripts specifically upregulated in MRT and haemolymph of both LT and MV was comparable (Fig. 2). This outcome, again, might suggest a conservation of transcriptomic responses in these two tissues to provide fundamental functions in both species, almost certainly including those related to CA metabolism.

Gene expression and role of the male reproductive system in CA biosynthesis

Overall, the transcriptome of MRT of both species encompassed a large set of protein domains whose enrichment was expected based on the fundamental role that these organs play in insect reproduction. Dynein, kinesin, microtubule-associated domains were the most abundant protein families, plausibly involved in mobilization of sperms, but also of secretory vesicles (and of their associated lipid and protein-based content) that were observed to replenish MAG [47]. Many other domains were associated to germ-cell proliferation or, more in general, to those developmental processes ongoing in reproductive organs of male insects, usually mediated by phosphorylation, nucleic acid regulation, ubiquitination and induction of many signalling cascades (as in T. castaneum) [48]. Among the latter, we observed many kinases, as well as a group of neuralized-like ubiquitin ligases known to play roles in Notch pathway-mediated cell fate decisions during development [49]. Furthermore, we also recorded carboxypeptidases that, along with other peptidases (e.g. aminopeptidases and trypsins in LT), are known to play a role in various aspects of male reproduction [50, 51].

Besides their reproductive functions, MRT, and particularly MAG, store large quantities of CA [13, 16,17,18,19]. However, MAG were recently disregarded as the site of CA production, as the enzymes involved in the earlier steps of CA biosynthesis were highly expressed in other organs [19, 52]. In fact, HGMCR - a key enzyme of the MVA pathway that directly alter the amount of downstream CA biogenesis in blister beetles [34] - showed a 60-fold higher expression in the CA-replenished fat body of E. sibirica than in other tissues [19]. Additionally, FDS, a key enzyme of the isoprenoid pathway synthesizing the farnesyl diphosphate (FPP), a precursor of various essential isoprenoid end-products, was hitherto found upregulated in the alimentary canal of males of H. cichorii [52]. Coherently with these previous observations, HMGCR and FDS (as well as unigenes encoding enzymes of the farnesol/farnesal-related branch of the terpenoid synthetic pathway; Additional file 5: Table S15) were both downregulated in MRT of LT and MV (Fig. 3B; Additional file 5: Table S15). Remarkably, among all, FDS appeared highly expressed in whole bodies of both species (and particularly in MV; see Fig. 3 and Additional file 5: Table S15). Hence, as already suggested [19, 52], FDS represents a plausible candidate gene for the synthesis of CA (along with other terpenoids) in blister beetles.

It has also been hypothesized that CA may originate from a JH-derived metabolite [35, 40]. In fact, silencing of two enzymes of the JH pathway, i.e. MFE and JHEH (but not JHAMT) caused a significant decrease in CA concentration in 3-to-7 days post-mated males of E. sibirica [35]. Hence, the JH III acid diol (a degradation product of JH III; Fig. 3A) was hypothesized as a potential precursor of CA [35]. JHs are sesquiterpenoids normally synthesized by corpora allata and released into the haemolymph where they regulate larval growth and metamorphosis [53]. In adults, JHs stimulate vitellogenesis and egg production in females [54] and orchestrate both the courtship behaviour and maturation of accessory glands in males [55, 56]. Since MRT of blister beetles can be reached by JH through the haemolymph, a degradation step of this hormone favouring a downstream CA synthesis could theoretically take place in this organ. Transcripts encoding two enzymes involved in the de novo synthesis of JH, i.e. MFE and JHAMT, were overall scarcely expressed (FPKM < 1) in LT and MV (Fig. 3B). This could imply that the JH titre at time of collection on field of (haemolymph-exuding) LT and MV was already adequate to ensure both the onset of reproductive maturation of males and CA renewal in their bodies. On the contrary, the mRNA of the JH-degrading enzyme JHEH was abundantly expressed in whole bodies of males of both species (Fig. 3B). The active conversion of JH into JH acid diol promoted by JHEH might serve both to suppress JH signalling and response [57] and supply JH-related metabolites to sustain CA biosynthesis. Since both our RNA-seq and RT-qPCR analyses confirmed that JHEH was significantly downregulated in MRT of both species (Fig. 3B-C), this might indicate that the JH-degrading route leading to CA was not expressed very highly in MRT at the time the samples were captured and dissected (i.e. CA may not be synthesized during the actual defence response, but accumulated before emission) and/or that it may more likely take place elsewhere in blister beetle body. Both options are plausible as in males of E. sibirica the highest transcript level of JHEH was observed 5th day after mating - in a phase of high CA productivity - and more significantly in the fat body [35].

Despite our data are coherent with previous findings downsizing the importance of MAG in the biosynthesis of CA, these might still constitutively express a key enzyme converting a linear terpenoid precursor into the typical tricyclo-decane structure of CA in one of the very latest steps of the biosynthetic pathway. Although we failed to detect transcripts pertaining to IDS-related terpene synthases - a family of enzymes known to perform terpene cyclization in insects [45] – some other overlooked biochemical processes converting a final metabolite delivered through the haemolymph, could eventually take place in MAG and contribute to accumulate CA in this organ. To conclude, it is worth to remark that tailored genomic analyses on fat bodies will be necessary to possibly identify the missing enzymes of the endogenous biochemical pathway. Anyhow, the here-provided transcriptomes of MRT still could represent a valuable background to uncover some still-unexplored enzymatic mechanisms related to CA de novo biosynthesis in blister beetles.

Possible mechanisms preventing self-intoxication from CA in blister beetles

CA is a potent toxin binding to protein phosphatase 2A, a metabolic constituent of all eukaryote cells [58]. Then, regardless of where CA is synthesized, mechanisms to cope with its toxicity should exist in the accumulating MAG, but, in general, in tissues and organs that interact with this terpene. In insects, toxicokinetic detoxification of lipophilic xenobiotics is typically achieved in two phases: i) Phase I, consisting in oxidation-reduction reactions converting the xenobiotic into a polar compound and generally performed by Cyp450s, carboxylesterases, alcohol/aldehyde dehydrogenases, hydroxylases and peroxidases, and ii) Phase II, consisting in the subsequent conjugation of the (polar) degradation product with an endogenous compound (sugars, amino acids, or glutathione) usually accomplished by glucosyl-transferases, acyl-transferases and glutathione S-transferases (GST) [59,60,61]. Logically, it would not be beneficial for blister beetles to activate Phase I and Phase II detoxification strategies in CA storing organs, because that would alter or inactivate a de novo synthetized and advantageous terpene. Consistently, although Cyp450s and other Phase I and II enzyme-related domains were abundant in whole transcriptomes of both LT and MV (Fig. 1A) - as expected in phytophagous insects protecting themselves against the detrimental effects of secondary metabolites of plants [62] - these were not significantly over-represented in MRT (Additional file 4: Table S11-S12). Then, other self-protective strategies are more likely than direct detoxification, such as sequestration of toxic compounds in intracellular compartments (Phase III), an additional mechanism which several insects evolved to tolerate xenobiotics [63,64,65]. Indeed, the observation of numerous replenished vesicles inside the cells of the glandular epithelia of MAG and vasa deferentia of M. proscarabeus [47] could support a compartmentalization strategy of CA.

Phase III is operated in many organisms by membrane transport proteins belonging to ATP-binding cassette (ABC) transporters or solute carrier proteins of the Major Facilitator Superfamily (MFS) [66, 67]. ABC transporters are amongst the largest carrier families present in many phyla [68] and known to be involved in detoxification in insects [69], whereas only recently MFS have been demonstrated to mediate xenobiotic resistance in arthropods [70,71,72]. Again, in phytophagous insects, proteins of both families can prevent the interference of chemical plant defences with physiological processes [46, 73]. The inspection of LT and MV orthologous unigenes upregulated in MRT allowed detecting some conserved transcripts encoding ABC and MFS proteins, as well as members of other protein families involved in the mobilization of small solutes (Fig. 4A; Additional file 5: Table S16). These transporters might function in the loading of MAG or other secretory reproductive tissues with seminal fluid components (e.g. prostaglandins, lipids, peptides, hormones) [74, 75]. Along with transportation of seminal fluid substances, some of these proteins could have also evolved the ability in blister beetles to sequester and concentrate CA (or noxious CA precursors) into the vesicles of MAG to ensure protection of MRT. More likely, to prevent damages, ABC or MFS proteins could shuttle into MAG a chemically deactivated form of CA carried by the haemolymph, such as, for example, an endogenously produced glycosylated intermediate. Alternatively, the inward transportation of a diet-acquired glycosyl donor from haemolymph could serve to deactivate a noxious CA precursor in MRT. Noteworthy, both the complete absence of a cuticular partitioning and the labyrinthine lacunar system observed in the third MAG pair of Meloe proscarabeus [47] might support the absorption and vehiculation of substances from the haemolymph into blister beetle MRT.

Among LT and MV orthologous transporters upregulated in MRT, the proton-coupled folate transporter, sugar transporter SWEET1, sodium-coupled monocarboxylate and choline transporter-like protein 1 (Fig. 4A-B) echoed the panel of proteins governing the sequestration of the plant-derived phenolic glycoside salicin in larvae of Chrysomela populi [76]. By feeding on plants, larvae of this leaf beetle sequester salicin that reaches a specialized closed and chitin-coated defensive glandular reservoir via the haemolymph. In larval defensive glands the acquired glycoside is then enzymatically converted into a biologically active form, the salicylaldehyde, which is emitted as a deterrent [77,78,79]. Actually, the enzymatic conversion of a glycosidic precursor gathered directly from plants is considered a recent novelty evolved in leaf beetles to release toxins in a safe manner. In more ancestral leaf beetle groups, such as in Phaedon cochleariae Fabricius, 1791 and Gastrophysa viridula De Geer, 1775, monoterpenoids (iridoids) are de novo synthetized [80,81,82] through the conversion of an autonomously-produced 8-hydroxygeraniol-8-O-β-d-glucoside into the iridoids in the glandular reservoir, where cleavage and toxin release occur [81]. Likewise, to readily avoid the hazardous effects of an endogenously produced toxic terpene, glycosylation of a CA precursor should occur autonomously in blister beetles, without having to rely on sequestering compounds from diet. Therefore we hypothesize that, in order to circumvent the auto-intoxicative effects of CA, blister beetles might shuttle autonomously-synthetized CA glycosylated precursors from haemolymph into MAG. Remarkably, the enrichment in both species of glycosyl transferase domains in the transcriptomic repertoire of haemolymph (i.e. pf03360 in LT and Pf04666 in MV, Additional file 4: Table S13-S14) could support the transfer of a saccharide moiety to a circulating CA precursor.

Five other unigenes of ABC and MFS transporters which might play a role in sequestration mechanisms were identified using a more ‘relaxed’ threshold (i.e. FDR < 0.01) (Fig. 4A,C). Specifically, three of them (PLT4, l(2)03659, ABCF2) were significantly upregulated in MRT as compared to whole body (Fig. 4C) and, then, deserve further attention to understand their role in detoxification.

Whatever the biochemical mechanism to prevent damages in the storing organs, the auto haemorrhage response requires a rapid release of CA in the haemolymph to reach the intersegmental membranes of blister beetle appendages and be discharged externally. Then, once the behavioural response is triggered, a mode of transport in the haemolymph capable of mitigating the toxicity of a released and free-circulating CA in beetle tissues should also be expected. Among the proposed carriers, lipocalins have been regarded as potential CA-binding proteins able to simultaneously carry out roles related to both transport and detoxification in blister beetle haemolymph [21, 22]. Insect lipocalins bind a variety of lipophilic (endogenous or exogenous) compounds and are involved in many functions, among which olfaction, pheromone transport and metabolite (e.g. retinoic acid) binding [83]. Among the three LT and MV orthologous lipid-binding proteins upregulated in haemolymph, one was related to a phorbol ester/diacylglycerol binding protein unc-13 (whose possible role is discussed in the next paragraph), but two were, as previously hypothesized [21], associated to lipocalins (Fig. 4A). We specifically identified: (i) a possible homolog of apolipoprotein D (ApoD) (Fig. 4A), a lipocalin associated with high-density lipoproteins playing fundamental roles in lipid metabolism or in binding for progesterone and arachidonic acid in humans [84] and, (ii) a low-density lipoprotein receptor (MALRD1) (particularly abundant in LT haemolymph) able to bind to ligands (such as lipocalins) and internalize them through receptor-mediated endocytosis [85]. Among these two, ApoD was previously found transcriptionally induced in T. castaneum in response to the injection of crude lipopolysaccharide [86] or upon septic injury [87], thus suggesting a role of this protein in both immunity and stress responses in tenebrionid beetles. Hence, under stressful conditions, ApoD could be strongly upregulated in haemolymph of blister beetles. Once there, this protein could bind the released CA and deliver it outward during the autohaemorrhaging exudation. Since substrate binding abilities of lipocalins are broad [83], more tailored approaches are needed to confirm this hypothesis.

Haemolymph-specific transcriptomes and possible molecular adaptation to reflex-bleeding

Blister beetles, if stressed, exude viscous-oily droplets of haemolymph from the autohaemorrhaging tissues, and this represents the manner they externally emit CA.

The autohaemorrhaging behaviour requires efficient coagulation and integument repairing mechanisms to minimize haemolymph loss [5]. Moreover, wounds in the intersegmental membranes provoked by increased haemolymph pressure are exposed to the attack of pathogens, which challenge the immunity system.

In insects the haemostatic system serves to stop bleeding, prevent microbial entry and favour wound sealing and healing. The first step in clotting formation in arthropods is achieved by the release of structural clot components following haemocyte degranulation. Among the most represented motifs in LT and MV haemolymph proteins (Table 5), we identified Von Willebrand Factors (vWF), C8 and trypsin inhibitor-like cysteine rich domains, typical of structural clotting proteins of arthropods [88, 89]. We also retrieved several transcripts encoding haemocytins - multidomain insect humoral lectins homologous to the mammalian vWF [90, 91] - owning vWF, C8 and trypsin inhibitor-like cysteine rich domains (Table 5) and that likely represent the main structural glycoproteins involved in coagulation of blister beetle haemolymph. We also found many transcripts related to calcium-dependent transglutaminases (i.e. haemocyte protein-glutamine gamma-glutamyltransferase) and prophenoloxidases (PPO), two protein classes known to cross-link the structural components in insects to harden the clot [88, 89]. The PPO-activating system constitutes an important component of insect defence mechanisms [92]: in fact, PPO released by haemocytes via a proteolytic cascade promotes the oxidation of phenolic molecules to produce melanin around invading pathogens and wounds [88, 89, 93]. We also observed other factors possibly engaged in coagulation and innate immunity, among which Spätzle and Laminin. Spätzle is a structural homolog of coagulogen, a clotting protein and functional equivalent of fibrinogen from horseshoe crab [94], but also a key regulator of the Toll pathway, leading to the expression of genes involved in immune defence to gram-positive bacteria and fungi [95]. Laminins are extracellular matrix glycoproteins of the basal lamina [96] but were demonstrated to interact with invading parasites and to play a dual role in immunity by both maintaining the basal levels of complement in the haemolymph and promoting the production of complement components through the interaction with LRIM1, as demonstrated in Anopheles gambiae Giles, 1902 [97].

Filamin, fascin, mucin and vinculin were among the most abundant domains in LT and MV haemolymph (Table 5; Additional file 4: Table S13-S14). Since these domains are typical of actin-bundling proteins producing gelation factors [98, 99], they could be responsible for the emission of the exuded haemolymph in the form of a mutable (liquid to solid) viscous-elastic oily substance in blister beetles. However, we cannot exclude that such factors might somehow be involved in tissue morphogenesis and repair following the autohaemorrhaging damage [100]. These repairing functions could also be performed by thymosines, whose domains were highly represented in the haemolymph of the two analysed species (Table 5). Thymosin β4, in particular, is associated with tissue repair and cell migration in vertebrates [101], albeit still scarcely characterized in insects [102]. Thymosines are also known to regulate haemopoietic stem cell proliferation and differentiation in crustaceans [103]. Indeed, haematopoiesis is expected to be enhanced after injury-induced haemolymph loss and nonself challenge [104] and an increased level of mitosis is also likely to occur within the hematopoietic tissue by immune stimulation [105]. Coherently, we found many domains related to mitosis (Table 5), such as mini-chromosome maintenance (MCM) proteins involved in chromosome replication.

We also detected some over-represented protein domains related to the diacylglycerol (DAG) /phosphatidic acid (PA) signalling pathway [106], i.e. diacylglycerol kinases (RalGDS/AF-6 and C1 domains), pleckstrin homology (PH) and PDZ-related proteins (Table 5; Additional file 4: Table S13-S14). The high representation of diacylglycerol kinases (and the upregulation of inositol 1,4,5 triphosphate phosphatase-related transcripts; Additional file 2: Table S4, S6) could indicate an active conversion of DAG to PA [107], the latter being the main substrate of two major constituents of cell membranes in eukaryotes, i.e. phosphatidylcholine and phosphatidylethanolamine [108]. PA is also known to activate the phosphoinositide signalling pathway [109] regulating many cell activities through the direct interaction of phosphoinositides with membrane proteins or by membrane recruitment of cytosolic proteins containing e.g. PH and PDZ domains able to bind phosphoinositides [110, 111]. Phosphoinositide signalling triggers cell proliferation and survival, but also induces cytoskeletal changes and actin remodelling for vesicle trafficking, membrane dynamics (and ruffling), cell division/cytokinesis and migration [112]. Hence, among all possible roles of DAG/PA signalling pathway, the increase of PA and the activation of phosphoinositide-induced processes in our target species could respond to the need of repairing cells and remodelling tissues damaged after autohaemorrhaging. These processes could also be sustained by the activity of multiple Rap/Rho GAP and Rho GEF-containing proteins (Table 5; Additional file 4: Table S13-S14), known to regulate cytoskeletal rearrangements necessary for cell-shape change, cell adhesion and migration [113]. More tailored investigations, however, are needed to clarify if the activation of the above-mentioned signalling pathways could have any role in blister beetle autohaemorrhaging response.

Finally, we observed several proteins exhibiting multiple C2 domains (Table 5). These Ca2+-dependent cysteine-rich modules, binding phospholipids and targeting cell membranes, are involved in several signal transduction cascades or membrane trafficking [114]. In LT and MV haemolymph C2 domains were mostly represented by transcripts annotated as Multiple C2 domain/Transmembrane region Proteins (i.e. MCTP, containing three consecutive C2 and a peroxin-like domain), synaptotagmin-1, phorbol ester/diacylglycerol binding protein unc-13 homolog B and Protein Kinase C. The exact homology and function of these transcripts need to be verified, since the presence of multiple-C2 domain proteins, synaptotagmin-1 and unc-13 homologs in insect haemolymph is rather unusual. In fact, in invertebrates these molecules are expressed in synaptic vesicles of the nervous system, where they regulate the secretion of neurotransmitters [115, 116]. Noteworthily, MCTPs were found expressed in the accessory cells of the olfactory organs of Drosophila and supposed to be secreted into the sensillum lymph that surrounds the olfactory receptor neuron dendrites [117]. If further analyses will confirm the homology of the retrieved transcripts with proteins involved in neuronal functions, then circulating MCTP, synaptotagmin-1 and unc-13 could neuromodulate the response to autohaemorrhaging ruptures starting from peripheral sensory organs placed in legs and antennae where these damages more often occur. These findings, if corroborated, could inform about the existence of a ‘neuronal alarm system’ allowing blister beetles to rapidly activate repairing functions and minimize haemolymph loss.


This study provided a high quality de novo assembly of LT and MV male transcriptomes, a valuable genetic resource to explore some of the still enigmatic aspects of blister beetle biology. By specifically producing transcriptomes of two tissues steadily in contact with CA, i.e. the MRT and haemolymph in two species, we were able to draw some hypotheses on the mechanisms employed by these insects to guarantee a safe storage, circulation and autohaemorrhaging emission of CA. We were also able to identify for the first-time a panel of haemolymph factors expressed in the reflex-bled exudate of blister beetles. Finally, we contributed to characterize the expression levels in MRT and haemolymph of genes so-far known to be involved in the autogenous production of CA in two species never investigated thus far.

The precise mechanisms allowing CA to be accumulated in MRTs and released into the haemolymph without causing damages to internal tissues of blister beetles are still elusive. In this regard, we here propose that sequestration through ABC, MFS or solute transporters of a purported glycosylated CA precursor could avoid auto-intoxication after accumulation in MAG of this terpene. We also identified some abundantly expressed lipocalins that could vehiculate and mitigate the reactivity of a freely circulating CA in the haemolymph during the autohaemorrhaging response. Future tailored studies aimed at verifying these hypotheses and/or definitively identifying factors granting self-protection in blister beetles would shed light on one of the most intriguing aspects of their biology. Nonetheless, unveiling these mechanisms would be of great help to design innovative drug-delivery systems for a safe therapeutic application of CA in medicine. As for gene upregulation in haemolymph, we detected many factors involved in coagulation and integument repairing mechanisms in blister beetles. Their over-representation possibly reflects the need to limit the fluid loss due to the reflex-bleeding behaviour, a possible common adaptation in autohaemorrhaging insects [5] which is worth to be further analysed in an evolutionary perspective. Finally, respect to CA biosynthesis, our results are coherent with recent literature excluding a leading role of MRT in CA production at both its upstream biosynthetic steps (i.e. MVA pathway) and those related to JH catabolism. Yet, in MAG, overlooked biochemical mechanisms converting a final terpenoid intermediate into the tricyclo-decane skeleton of CA cannot be excluded. Future research efforts shall be focused to elucidate the final phases of the CA biosynthetic pathway and candidate organs where these latest enzymatic steps may occur.


Beetle sampling

Males of LT (tot. n. = 45, all labeled: Lazio, Roma, Tolfa, Rio Fiume, 42°03′36″N 11°56′50″E, 75 m) and of MV (tot. n. = 69; 56 labeled: Lazio, Roma, Tolfa, Piantangeli - Grasceta dei Cavallari, 42°10′52″N 11°56′30″E, 440 m; 13 labeled: Abruzzo, Pescara, Popoli, R.N. Sorgenti del Pescara, 42°10′41″N 13°48′36″E, 350 m) were collected in Central Italy in June–July 2019 by hand-picking the specimens while resting on flower stems (mainly of Dipsacaceae and Asteraceae) during the warmer hours of the day. Specimens were identified on field by the taxonomist of our research team (MAB) using taxonomic keys (Bologna 1991). Conspecific individuals were sexed by examining the last ventrite shape (deeply emarginated in males) and immediately placed in fauna-boxes with either a paper or soil bottom and provided with sand. All specimens were then transported to laboratory and kept alive for 1–2 days or, at most, a week, and fed daily with fresh flowers (e.g. Knautia sp.) and fruits (e.g. apple slices).

Collection of tissues

The ‘whole body’ (or ‘total transcriptome’) reference set was constructed by immediately storing in TRIzol Reagent (Thermo Fisher Scientific, Wilmington, DE, USA) 9 males of both LT and MV (18 specimens in total). Biological replicates were obtained by sub-dividing the stored males in 3 pools, each composed by 3 individuals per species. Haemolymph from the remaining 36 LT and 60 MV was gathered after cutting the apical antennomers and applying a light pressure on the abdomen to facilitate the release of droplets from the cut areas. In case of poor haemolymph emission, cuts were also performed in the tarsal part of one or more legs. Haemolymph droplets were collected in ice-kept 1.5 ml Eppendorf vials filled with 1 ml TRIzol Reagent. Each of the three biological replicates of haemolymph was assembled from 12 LT and 20 MV, respectively. After haemolymph extraction, the MRT from 15 (out of 36) LT and 15 (out of 60) MV were dissected, and 3 biological replicates composed by 5 of each were generated. Each specimen was put in a sterile Petri dish and quickly dissected with scissors and tweezers under a stereomicroscope. Ventrites were removed to expose the male internal genitalia and a drop of RNase-Free water was added to facilitate the isolation of MAG from the contiguous internal organs. After their dissection, the three pairs of MAG, testicles and vasa deferentia were gently pulled out with tweezers and put in 1.5 ml vials (Eppendorf) kept on ice with 1 ml TRIzol.

RNA extraction, quantification and sequencing

Total RNA from tissues was isolated using a TRIzol-based procedure (Thermo Fisher Scientific, Wilmington, DE, USA). Total bodies were grounded in a sterile ceramic mortar using liquid nitrogen to obtain a fine powder, and then homogenized in 1 ml of TRIzol Reagent. Pooled tissues of haemolymph and MRT were lysed directly in TRIzol by pipetting the lysate several times until complete homogenization. RNA concentration and purity were determined by measuring absorbance using NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). 1 μg of total RNA was run on a 1% denaturing gel to verify RNA integrity. For each sample (‘pool’, see above) 10 μg (from total body or MRT) or 1–4 μg (from haemolymph) of high-quality RNA (RIN value > 8) was sent to IGA Technology Services s.r.l. (Udine, IT) for mRNA-seq stranded library preparation, validation, and sequencing, resulting in 9 libraries for each species. A 2 × 150 bp paired-end sequencing was performed using a NovaSeq 6000 System (Illumina, San Diego, California, USA) with a depth of about 30 million (haemolymph and MRT) or 60 million (total body) of reads per sample. Data generated in the present study are available in the Sequence Read Archive (SRA) database of NCBI (http://www. under bioproject number PRJNA674987.

De novo transcriptome assembly, abundance estimation and differential expression analysis

Read quality was assessed by FastQC software v0.11.4 ( and read quality trimming was performed by Trimmomatic software (v0.32) [118]. The whole quality-trimmed read dataset (including reads from all tissues and body samples) for each species, concatenated into two paired FASTQ files, was de novo assembled using the Trinity software (release 2.3.2) [119, 120] with default parameters and --SS_lib_type RF, −-jaccard_clip, −-normalize_reads, −-min_contig_length 400 flags set on the ADA Server at the Department of Biology, University of Naples Federico II (24 cores, 256 GB of memory) [121]. The quality of the assembled transcriptomes was assessed by BUSCO (v3.0.1) pipeline [122]. Transcript-level quantification for each sample was performed using RSEM software [123] and Bowtie aligner [124], as implemented in the Trinity software package. Quality control between RNA-seq replicates was performed using the PtR Trinity perl script (v0.32); very good correlation scores (R2 > 80) were obtained for all replicates in both species, except for the third replicate from total body tissue of LT that was excluded from subsequent analyses. DGE analysis was performed using edgeR software [125], which uses a negative binomial model for differential expression analysis, with cut-off values of FDR (False Discovery Rate) < 0.001 and FC (Fold Change) > 2. FDR = p-value adjusted for multiple testing with the Benjamini-Hochberg procedure.

Functional annotations

The functional categories present in the two assembled transcriptomes were investigated and summarized using Annocript pipeline ( [126] and UniProtKB reference database. The longest transcript of each Trinity gene cluster was employed as representative for the annotation. We performed the following similarity searches: BLASTX against TrEMBL/UniRef and SwissProt (parameters: e-value 1E-5, threshold 18, wordsize 4), RPSBLAST against CDD profiles (parameters: e-value 1E-5), BLASTn against Rfam and rRNAs (parameters: e-value 1E-5). GO (Gene Ontology) and PFAM annotations were obtained, for each transcript, using default Annocript pipeline parameters. The enrichment analysis was then performed on the GO and PFAM terms of DE transcripts, identified by edgeR in each pairwise tissue comparisons and for each species, using Annocript and the Fisher Exact Test (adjusted p-value < 0.01) in [R] package (

Manual BLAST searches for genes putatively involved in CA biosynthesis

Manual searches were performed to detect unigenes that coded enzymes supposed to be involved in CA regulation in blister beetles and mapped in the ‘terpenoid backbone’ (map 00900) and ‘insect hormone’ (map 00901) biosynthetic pathways of Kyoto Encyclopedia of Genes and Genomes (KEGG; tBLASTn searches were run on LT and MV transcriptomes by using Tribolium (Tenebrionoidea: Tenebrionidae) or preferentially (when available) Epicauta or Hycleus (Tenebrionoidea: Meloidae) protein sequences deposited in GenBank (and reported in detail in Additional file 5: Table S15) as a query. Among transcripts showing the highest matching scores (and lowest e-values), only those aligning with more than a half of residues of the query length and showing an identity percentage > 60% were considered. To confirm these results, manual BLASTn searches were also performed using the same matching threshold against LT and MV transcriptomes for the subset of genes whose mRNA sequences were available in Genbank for Epicauta and/or Hycleus. Transcripts with the highest scores and identity percentages for the searched gene (thus, representing the most plausible closest orthologs) and with FPKM> 1 for at least one of the three entries i.e. body (B), haemolymph (H) and male reproductive tract (MRT), were retained; if all variants of the “best orthologs” scored FPKM< 1, these were all reported.

Comparative analysis of LT and MV orthologs to identify genes involved in CA deactivation and/or transport

Based on the assumption that proteins responsible for CA detoxification in the accumulating organs should be conserved across species, we optimized the selection of candidate genes by searching for putative orthologous transcripts upregulated (FDR < 0.01) in MRT versus body samples of both LT and MV. Similarly, we searched for orthologous transcripts upregulated in haemolymph versus body samples of both species to identify possible transport-related proteins. The two MRT upregulated transcript lists, and the two haemolymph upregulated transcript lists were pairwise compared using a Best Reciprocal BLASTn Hit approach with a custom Perl script and a coverage cutoff of 30% and an e-value cut-off of 0.01. To detect potential CA detoxifying enzymes or transporters, we selected common upregulated transcripts in MRT and haemolymph containing functional protein domains typical of: i) enzymes participating to the modification/degradation of apolar xenobiotics (e.g. CypP450s, carboxylesterases, hydroxylases, peroxidases, GSTs, glycosyl transferases and acyltransferase), ii) transporters involved in sequestration of terpenes/drugs and/or xenobiotics in insects (e.g. ABC-transporters, Multidrug-resistance and MFS proteins [66, 127] and iii) small solute carriers. Additionally, we searched for common upregulated transcripts in haemolymph encoding proteins previously hypothesized to bind and mobilize CA, such as lipophorins, apolipoprotein, lipocalins and General Odorant Binding Proteins [21, 22]. After manual inspection, we selected candidate orthologous transcripts from a broader list provided by the ‘Best Reciprocal Hit’ output including 483 and 283 unsorted pairs of transcripts from gonads and haemolymph, respectively. In doing so, all “unknown transcripts” which failed to receive a functional domain by the automatic annotation procedure were further scanned using InterProScan [128] and HHpred [129].

Validation of transcripts by RT-qPCR

Total RNA was extracted from whole bodies and dissected MRT to obtain three novel replicates for RT-qPCR analysis. 1 μg of total RNA extracted from MRT and whole body samples was retro-transcribed into cDNA by SuperScript™ III Reverse Transcriptase (Thermo Fisher Scientific, Wilmington, DE, USA), following the manufacturer’s instructions. The RT-qPCR was performed using the SsoAdvanced™ Universal SYBR® Green Supermix (BIO-RAD Laboratories Inc., Hercules, CA, USA), according to manufacturer’s instructions. Amplifications were conducted for: i) β-Actin (β-ACT), Ribosomal Protein S7 (RP S7), as reference genes; ii) Juvenile hormone acid O-methyltransferase (JHAMT), Methyl farnesoate epoxidase (MFE), Juvenile hormone epoxide hydrolase (JHEH) of the JH pathway; iii) Proton-coupled folate transporter (PCFT), Sugar transporter Sweet 1 (SWEET1), Sodium monocarboxylate cotransporter 1 (SMCT1), Choline transporter-like protein 1 (SLC44A1), as potential orthologs of salicin sequestring proteins of C. populi); iv) Polyol transporter 4 (Polt4), Probable multidrug resistance-associated protein lethal (2)03659 (l(2)03659), Multidrug resistance protein homolog 49 (MDR49), ATP-binding cassette sub-family F member 2 (ABCF2), ATP-binding cassette sub-family B member 6 (ABCB6), as other (MFS/ABC) transporters of interest. Specific primers pairs were designed (Additional file 6: Table S17). All reactions were performed in triplicate in the Corbett Rotor-Gene 6000 (Qiagen, Hilden, Germany) and relative quantification was carried out with the ΔΔCT method [130] using the abundance of RPS7 mRNA as endogenous housekeeping control. The relative transcription levels as obtained by RT-qPCR analyses were therefore compared with abundance levels as detected by RNA-seq (FPKM, fragments per kilobase of exon model per million reads mapped - the same values employed to produce the Heat Maps in Fig. 3-4). Values from replicate experiments were averaged. Finally, to obtain values suitable for statistical comparisons, we calculated (for each gene) a fold-change (FC) value as the ratio of abundance of each transcript (ΔCt and FPKM for RT-qPCR and RNA-seq, respectively) over the group median. These values (plotted after conversion in log2 numbers) were used to evaluate the correlation between RNA-seq and RT-qPCR methods, applying statistical evaluation throughout the Pearson test (by using the Prism GraphPad software).

Availability of data and materials

The short-read DNA sequences have been deposited in the Sequence Read Archive (SRA) database of NCBI (http://www. under bioproject number PRJNA674987.



ATP-binding cassette


Alkaline phosphatase


Apolipoprotein D


acetyl-CoA C-acetyltransferase


Basic local alignment search tool


Benchmarking set of universal single-copy orthologues




Cytochrome P450




Differentially expressed genes


False discovery rate


Farnesyl diphosphate synthase


protein farnesyltransferase subunit beta


NADP+-dependent farnesol dehydrogenase


Fragments per kilobase million


Farnesyl diphosphate


Gene Ontology


Glutathione S-transferase


3-hydroxy-3-methylglutary-CoA reductase


hydroxyl-methylglutaryl-CoA synthase


protein-S-isoprenylcysteine O-methyltransferase


Kyoto encyclopedia of genes and genomes


Juvenile hormone


JH acid O-methyltransferase JH-III synthase


Juvenile hormone epoxide hydrolase


Lydus trimaculatus


Male accessory glands


MAM and LDL-receptor class A domain-containing protein


Mini-chromosome maintenance


Mylabris variabilis




Mevalonate kinase.


phosphomevalonate kinase


Methyl farnesoate epoxidase


Major Facilitator Superfamily


Multiple C2 domain and Transmembrane region Proteins


male reproductive tract


diphosphomevalonate decarboxylase


Phosphatidic acid


Protein Family


Pleckstrin homology




Reverse transcription-quantitative PCR


STE24 endopeptidase


UniProt reference clusters


von Willebrand Factor


  1. Hollande CH. L’autohémorrhée ou le rejet du sang chez les insectes (toxicologie du sang). Archs Anal Microsc. 1911;13:171–318.

    Google Scholar 

  2. Holloway GJ, de Jong PW, Brakefield PM, de Vos H. Chemical defence in ladybird beetles (Coccinellidae). 1. Distribution of coccinelline and individual variation in defence in 7-spot ladybirds (Coccinella septempunctata). Chemoecology. 1991;2(1):7–14.

    Article  CAS  Google Scholar 

  3. Karystinou A, Thomas APM, Roy HE. Presence of haemocyte-like cells in coccinellid reflex blood. Physiol Entomol. 2004;29(1):94–6.

    Article  Google Scholar 

  4. Knapp M, Řeřicha M, Židlická D. Physiological costs of chemical defence: repeated reflex bleeding weakens the immune system and postpones reproduction in a ladybird beetle. Sci Rep. 2020;10(1):9266.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Stocks I. Reflex bleeding (autohemorrhage). In: Capinera JL, editor. Encyclopedia of. Entomology, Dordrecht: Springer; 2008.

    Chapter  Google Scholar 

  6. Bologna MA. Coleoptera Meloidae. Fauna D'italia. XXVIII. Calderini: Bologna; 1991.

    Google Scholar 

  7. Bologna MA, D'Inzillo B, Cervelli M, Oliverio M, Mariottini P. Molecular phylogenetic studies of the Mylabrini blister beetles (Coleoptera, Meloidae). Mol Phyl Evol. 2005;37(1):306–11.

    Article  CAS  Google Scholar 

  8. Bologna MA, Oliverio M, Pitzalis M, Mariottini P. Phylogeny and evolutionary history of the blister beetles (Coleoptera, Meloidae). Mol Phyl Evol. 2008;48(2):679–93.

    Article  CAS  Google Scholar 

  9. Bravo C, Mas-Peinado P, Bautista LM, Blanco G, Alonso JC, García-París M. Cantharidin is conserved across phylogeographic lineages and present in both morphs of Iberian Berberomeloe blister beetles (Coleoptera, Meloidae). Zool J Linnean Soc. 2017;180(4):790–804.

    Article  Google Scholar 

  10. Gisondi S, Gasperi T, Roma E, Tomai P, Gentili A, Vignoli L, et al. Cantharidin content in two Mediterranean species of blister beetles, Lydus trimaculatus and Mylabris variabilis (Coleoptera: Meloidae). Entomol Sci. 2019;22(3):258–63.

    Article  Google Scholar 

  11. Capinera JL, Gardner DR, Stermitz FR. Cantharidin levels in blister beetles (Coleoptera: Meloidae) associated with alfalfa in Colorado. J Econ Entomol. 1985;78(5):1052–5.

    Article  CAS  Google Scholar 

  12. Blodgett SL, Carrel JE, Higgins RA. Cantharidin content of blister beetles (Coleoptera: Meloidae) collected from Kansas alfalfa and implications for inducing cantharidiasis. Environ Entomol. 1991;20(3):776–80.

    Article  Google Scholar 

  13. Carrel JE, McCairel MH, Slagle AJ, Doom JP, Brill J, McCormick JP. Cantharidin production in a blister beetle. Experientia. 1993;49(2):171–4.

    Article  CAS  PubMed  Google Scholar 

  14. Nikbakhtzadeh MR, Tirgari S. Cantharidin component of Iranian blister beetles (col: Meloidae) and their differences between Iranian and exotic species. Iran J Public Health. 2002;31:113–7.

    CAS  Google Scholar 

  15. Mebs D, Pogoda W, Schneider M, Kauert G. Cantharidin and demethylcantharidin (palasonin) content of blister beetles (Coleoptera: Meloidae) from southern Africa. Toxicon. 2009;53(4):466–8.

    Article  CAS  PubMed  Google Scholar 

  16. Sierra JR, Woggon WD, Schmid H. Transfer of cantharidin during copulation from the adult male to the female Lytta vesicatoria (Spanish fly). Experientia. 1976;32(2):142–4.

    Article  CAS  Google Scholar 

  17. Nikbakhtzadeh MR, Hemp C, Ebrahimi B. Further evidence for the role of cantharidin in the mating behaviour of blister beetles (Coleoptera: Meloidae). Integr Biosci. 2007;11(2):141–6.

    Article  Google Scholar 

  18. Jiang M, Lü SM, Zhang Y. The potential organ involved in cantharidin biosynthesis in Epicauta chinensis Laporte (Coleoptera: Meloidae). J Insect Sci. 2017;17(2):2.

    Article  CAS  Google Scholar 

  19. Jiang M, Lü SM, Qi ZY, Zhang YL. Characterized cantharidin distribution and related gene expression patterns in tissues of blister beetles, Epicauta chinensis. Insect Sci. 2019;26(2):240–50.

    Article  CAS  PubMed  Google Scholar 

  20. Selander RB. Sexual behavior in blister beetles (Coleoptera: Meloidae): I. The Genus Pyrota. The Can Entomol. 1964;96(8):1037–82.

    Article  Google Scholar 

  21. McCormick JP, Carrel JE. Cantharidin biosynthesis and function in meloid beetles. In: Prestwich GD, Blomquist G, editors. Pheromone biochemistry. Reno: Academic Press; 1987. p. 307–50.

    Chapter  Google Scholar 

  22. Dettner K. Inter-and intraspecific transfer of toxic insect compound cantharidin. In: Dettner K, Bauer G, Völkl W, editors. Vertical food web interactions. Berlin: Springer; 1997. p. 115–45.

    Chapter  Google Scholar 

  23. Nikbakhtzadeh MR, Vahedi M, Vatandoost H, Mehdinia A. Origin, transfer and distribution of cantharidin-related compounds in the blister beetle Hycleus scabiosae. J Venom Anim Toxins Trop Dis. 2012;18(1):88–96.

    Article  CAS  Google Scholar 

  24. Eisner T, Rossini C, González A, Iyengar VK, Siegler MV, Smedley SR. Paternal investment in egg defence. In: Hilker M, Meiners T, editors. Chemoecology of insect eggs and egg deposition. Berlin: Blackwell Verlag GmbH; 2002. p. 91–116.

    Chapter  Google Scholar 

  25. Beauregard H. Les insectes vesicants. Paris: F. Alcan; 1890.

    Google Scholar 

  26. Wang GS. Medical uses of Mylabris in ancient China and recent studies. J Ethnopharmacol. 1989;26(2):147–62.

    Article  CAS  PubMed  Google Scholar 

  27. Moed L, Shwayder TA, Chang MW. Cantharidin revisited: a blistering defense of an ancient medicine. Arch Dermatol. 2001;137(10):1357–60.

    Article  CAS  PubMed  Google Scholar 

  28. Vakharia PP, Chopra R, Silverberg NB, Silverberg JI. Efficacy and safety of topical cantharidin treatment for Molluscum contagiosum and warts: a systematic review. Am J Clin Dermatol. 2018;19(6):791–803.

    Article  PubMed  Google Scholar 

  29. Liu D, Chen Z. The effects of cantharidin and cantharidin derivates on tumour cells. Anti Cancer Agents Med Chem. 2009;9(4):392–6.

    Article  CAS  Google Scholar 

  30. Verma AK, Prasad SB. Bioactive component, cantharidin from Mylabris cichorii and its antitumor activity against Ehrlich ascites carcinoma. Cell Biol Toxicol. 2012;28(3):133–47.

    Article  CAS  PubMed  Google Scholar 

  31. Deng LP, Dong J, Cai H, Wang W. Cantharidin as an antitumor agent: a retrospective review. Curr Med Chem. 2013;20(2):159–66.

    Article  CAS  PubMed  Google Scholar 

  32. Puerto Galvis CE, Vargas Mendez LY. Kouznetsov VV Cantharidin-based small molecules as potential therapeutic agents. Chem Biol Drug Des. 2013;82(5):477–99.

    Article  CAS  PubMed  Google Scholar 

  33. Naz F, Wu Y, Zhang N, Yang Z, Yu C. Anticancer attributes of cantharidin: involved molecular mechanisms and pathways. Molecules. 2020;25(14):3279.

    Article  CAS  PubMed Central  Google Scholar 

  34. Lü S, Jiang M, Huo T, Li X, Zhang Y. 3-hydroxy-3-methyl glutaryl coenzyme a reductase: an essential actor in the biosynthesis of cantharidin in the blister beetle Epicauta chinensis Laporte. Insect Mol Biol. 2016;25(1):58–71.

    Article  CAS  PubMed  Google Scholar 

  35. Jiang M, Lü S, Zhang Y. Characterization of juvenile hormone related genes regulating cantharidin biosynthesis in Epicauta chinensis. Sci Rep. 2017;7(1):1–11.

    Article  CAS  Google Scholar 

  36. Huang Y, Wang Z, Zha S, Wang Y, Jiang W, Liao Y, et al. De novo transcriptome and expression profile analysis to reveal genes and pathways potentially involved in cantharidin biosynthesis in the blister beetle Mylabris cichorii. PLoS ONE. 2016;11(1):e0146953.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Carrel JE, Doom JP, McCormick JP. Cantharidin biosynthesis in a blister beetle: inhibition by 6-fluoromevalonate causes chemical disarmament. Experientia. 1986;42(7):853–4.

    Article  CAS  PubMed  Google Scholar 

  38. Peter MG, Woggon WD, Schmid H. Identification of farnesol as an intermediate in the biosynthesis of cantharidin from mevalonolactone. Helv Chim Acta. 1977;60(8):2756–62.

    Article  CAS  PubMed  Google Scholar 

  39. Peter MG, Waespe HR, Woggon WD, Schmid H. Incorporation experiments with (3H and 14C) doubly labelled farnesols into cantharidin. Helv Chim Acta. 1977;60:1762–72.

    Google Scholar 

  40. Woggon WD, Hauffe SA, Schmid H. Biosynthesis of cantharidin: evidence for the specific incorporation of C-4 and C-11′ of farnesol. J Chem Soc Chem Commun. 1983;6(6):272–4.

    Article  Google Scholar 

  41. McCormick JP, Carrel JE, Doom JP. Origin of oxygen atoms in cantharidin biosynthesized by beetles. J Am Chem Soc. 1986;108(25):8071–4.

    Article  CAS  Google Scholar 

  42. Blum MS. Chemical Defence of arthropods. New York: Academic Press; 1981.

    Google Scholar 

  43. Nicolson SW. Water replenishment following reflex bleeding in the blister beetle Decapotoma lunata Pallas (Coleoptera: Meloidae). African Entomol. 1994;2:21–3.

    Google Scholar 

  44. Sutherland TD, Unnithan GC, Feyereisen R. Terpenoid omega-hydroxylase (CYP4C7) messenger RNA levels in the corpora allata: a marker for ovarian control of juvenile hormone synthesis in Diploptera punctata. J Insect Physiol. 2000;46(8):1219–27.

    Article  CAS  PubMed  Google Scholar 

  45. Beran F, Rahfeld P, Luck K, Nagel R, Vogel H, Wielsch N, et al. Novel family of terpene synthases evolved from trans-isoprenyl diphosphate synthases in a flea beetle. Proc Natl Acad Sci U S A. 2016;113(11):2922–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Oppenheim SJ, Baker RH, Simon S, DeSalle R. We can't all be supermodels: the value of comparative transcriptomics to the study of non-model insects. Insect Mol Biol. 2015;24(2):139–54.

    Article  CAS  PubMed  Google Scholar 

  47. Muzzi M, Di Giulio A, Mancini E, Fratini E, Cervelli M, Gasperi T, et al. The male reproductive accessory glands of the blister beetle Meloe proscarabaeus Linnaeus, 1758 (Coleoptera: Meloidae): anatomy and ultrastructure of the cantharidin-storing organs. Arthropod Struct Dev. 2020;59:100980.

    Article  PubMed  Google Scholar 

  48. Khan SA, Eggleston H, Myles KM, Adelman ZN. Differentially and Co-expressed Genes in Embryo, Germ-Line and Somatic Tissues of Tribolium castaneum. G3. 2019;9:2363–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Yeh E, Dermer M, Commisso C, Zhou L, McGlade CJ, Boulianne GL. Neuralized functions as an E3 ubiquitin ligase during Drosophila development. Curr Biol. 2001;11:1675–9.

    Article  CAS  PubMed  Google Scholar 

  50. Sonenshine DE, Bissinger BW, Egekwu N, Donohue KV, Khalil SM, Roe M. First transcriptome of the testis-vas deferens-male accessory gland and proteome of the spermatophore from Dermacentor variabilis (Acari: Ixodidae). PLoS ONE. 2011;6(9):e24711.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Gotoh A, Shigenobu S, Yamaguchi K, Kobayashi S, Ito F, Tsuji K. Transcriptome profiling of the spermatheca identifies genes potentially involved in the long-term sperm storage of ant queens. Sci Rep. 2017;7(1):5972.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zha S, Yin Y, Wang Y, Huang Y, Li Y, Wang Z. Cloning and functional analysis of farnesyl pyrophosphate synthase (FPPS) gene from Mylabris cichorii. Biotechnol Appl Biochem. 2017;64(5):667–76.

    Article  CAS  PubMed  Google Scholar 

  53. Riddiford LM. Cellular and molecular actions of juvenile hormone. I. General considerations and pre-metamorphic actions. Adv Insect Physiol. 1994;24:213–27.

    Article  CAS  Google Scholar 

  54. Wyatt GR, Davey KG. Cellular and molecular actions of juvenile hormone. II. Roles of juvenile hormone in adult insects. Adv Insect Physiol. 1996;26:2–155.

    Google Scholar 

  55. Wijesekera TP, Saurabh S, Dauwalder B. Juvenile hormone is required in adult males for Drosophila courtship. PLoS One. 2016;11(3):e0151912.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Konopova B, Jindra M. Juvenile hormone resistance gene Methoprene-tolerant controls entry into metamorphosis in the beetle Tribolium castaneum. Proc Natl Acad Sci U S A. 2007;104(25):10488–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Khlebodarova TM, Gruntenko NE, Grenback LG, Sukhanova MZ, Mazurov MM, Rauschenbach IY, et al. A comparative analysis of juvenile hormone metabolizing enzymes in two species of Drosophila during development. Insect Biochem Mol Biol. 1996;26(8-9):829–35.

    Article  CAS  PubMed  Google Scholar 

  58. Li YM, Casida JE. Cantharidin-binding protein: identification as protein phosphatase 2A. Proc Natl Acad Sci U S A. 1992;89(24):11867–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Kennedy CJ, Tierney KB. Xenobiotic protection/resistance mechanisms in organisms. In: Laws E, editor. Environmental Toxicology. New York: Springer; 2012.

    Chapter  Google Scholar 

  60. Khan RA, Rashid M, Wang D, Zhang YL. Toxicology and biochemical basis of cantharidin effects on Helicoverpa armigera (hub.) (Lepidoptera: Noctuidae). Pak J Zool. 2013;45:769–77.

    Google Scholar 

  61. Rashid M, Khan RA, Zhang YL. Over-expression of cytochrome P450s in Helicoverpa armigera in response to bio-insecticide, cantharidin. Int J Agriculture Biol. 2013;15:993–7.

    CAS  Google Scholar 

  62. Scott JG, Wen Z. Cytochromes P450 of insects: the tip of the iceberg. Pest Manag Sci. 2001;57(10):958–67.

    Article  CAS  PubMed  Google Scholar 

  63. Sorensen JS, Dearing MD. Efflux transporters as a novel herbivore countermechanism to plant chemical defenses. J Chem Ecol. 2006;32(6):1181–96.

    Article  CAS  PubMed  Google Scholar 

  64. Opitz SE, Müller C. Plant chemistry and insect sequestration. Chemoecology. 2009;19(3):117–54.

    Article  CAS  Google Scholar 

  65. Heckel DG. Insect detoxification and sequestration strategies. Voelckel C, Jander G, Insect-Plant Interactions. Section I. Biochemistry of Insect-Plant Interactions. New York: Wiley; 2014. 77–114, Insect Detoxification and Sequestration Strategies,

  66. Merzendorfer H. ABC transporters and their role in protecting insects from pesticides and their metabolites. In: Cohen E, editor. Advances in insect Physioloy, target receptors in the control of insect pests: part II, vol. 46. Cambridge: Academic Press; 2014. p. 1–72.

    Chapter  Google Scholar 

  67. Dermauw W, Van Leeuwen T. The ABC gene family in arthropods: comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Biol. 2014;45:89–110.

    Article  CAS  PubMed  Google Scholar 

  68. Holland IB. ABC transporters, mechanisms and biology: an overview. Essays Biochem. 2011;50(1):1–17.

    Article  CAS  PubMed  Google Scholar 

  69. Broehan G, Kroeger T, Lorenzen M, Merzendorfer H. Functional analysis of the ATP-binding cassette (ABC) transporter gene family of Tribolium castaneum. BMC Genomics. 2013;14(1):6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Dermauw W, Wybouw N, Rombauts S, Menten B, Vontas J, Grbić M, et al. A link between host plant adaptation and pesticide resistance in the polyphagous spider mite Tetranychus urticae. Proc Natl Acad Sci U S A. 2013;110(2):E113–22.

    Article  PubMed  Google Scholar 

  71. de la Paz C-MM, Wheat CW, Vogel H, Söderlind L, Janz N, Nylin S. Mechanisms of macroevolution: polyphagous plasticity in butterfly larvae revealed by RNA-Seq. Mol Ecol. 2013;22(19):4884–95.

    Article  CAS  Google Scholar 

  72. Snoeck S, Greenhalgh R, Tirry L, Clark RM, Leeuwen TV, Dermauw W. The effect of insecticide synergist treatment on genome-wide gene expression in a polyphagous pest. Sci Rep. 2017;7(1):13440.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Birnbaum SSL, Rinker DC, Gerardo NM, Abbot P. Transcriptional profile and differential fitness in a specialist milkweed insect across host plants varying in toxicity. Mol Ecol. 2017;26(23):6742–61.

    Article  CAS  PubMed  Google Scholar 

  74. Hull JJ, Chaney K, Geib SM, Fabrick JA, Brent CS, Walsh D, et al. Transcriptome-based identification of ABC transporters in the western tarnished plant bug Lygus hesperus. PLoS ONE. 2014;9(11):e113046.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Rangel J, Shepherd TF, Gonzalez AN, Hillhouse A, Konganti K, Ing NH. Transcriptomic analysis of the honeybee (Apis mellifera) queen spermathecae reveals genes that may be involved in sperm storage after mating. PLoS ONE. 2021;16(1):e0244648.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Schmidt L, Wielsch N, Wang D, Boland W, Burse A. Tissue-specific profiling of membrane proteins in the salicin sequestering juveniles of the herbivorous leaf beetle, Chrysomela populi. Insect Biochem Mol Biol. 2019;109:81–91.

    Article  CAS  PubMed  Google Scholar 

  77. Pavan M. Studi sugli antibiotici e insetticidi di origine animale. I. Sul principio attivo della larva di Melasoma populi L. (Col. Chrysomelidae). Arch Zool Ital. 1953;38:157–84.

    CAS  Google Scholar 

  78. Pasteels JM, Rowell-Rahier M, Braekman JC, Dupont A. Salicin from host plant as precursor of salicylaldehyde in defensive secretion of chrysomeline larvae. Physiol Entomol. 1983;8(3):307–14.

    Article  CAS  Google Scholar 

  79. Burse A, Boland W. RNAi based functional analysis of biosynthetic enzymes and transport proteins involved in the chemical defense of juvenile leaf beetles. In: Hoffmann KH, editor. Insect molecular biology and ecology. Boca Raton: CRC Press; 2014. p. 351–75.

    Google Scholar 

  80. Feld BK, Pasteels JM, Boland W. Phaedon cochleariae and Gastrophysa viridula (Coleoptera: Chrysomelidae) produce defensive iridoid monoterpenes de novo and are able to sequester glycosidically bound terpenoid precursors. Chemoecology. 2001;11(4):191–8.

    Article  CAS  Google Scholar 

  81. Burse A, Schmidt A, Frick S, Kuhn J, Gershenzon J, Boland W. Iridoid biosynthesis in Chrysomelina larvae: fat body produces early terpenoid precursors. Insect Biochem Mol Biol. 2007;37(3):255–65.

    Article  CAS  PubMed  Google Scholar 

  82. Strauss AS, Peters S, Boland W, Burse A. ABC transporter functions as a pacemaker for sequestration of plant glucosides in leaf beetles. eLife. 2013;2:e01096.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Flower DR. The lipocalin protein family: structure and function. Biochem J. 1996;318(1):1–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Rassart E, Bedirian A, Do Carmo S, Guinard O, Sirois J, Terrisse L, et al. Apolipoprotein D. Biochim Biophys Acta. 2000;1482(1-2):185–98.

    Article  CAS  PubMed  Google Scholar 

  85. Dantuma NP, Potters M, De Winther MP, Tensen CP, Kooiman FP, Bogerd J, et al. An insect homolog of the vertebrate very low-density lipoprotein receptor mediates endocytosis of lipophorins. J Lipid Res. 1999;40(5):973–8.

    Article  CAS  PubMed  Google Scholar 

  86. Altincicek B, Knorr E, Vilcinskas A. Beetle immunity: identification of immune-inducible genes from the model insect Tribolium castaneum. Dev Comp Immunol. 2008;32(5):585–95.

    Article  CAS  PubMed  Google Scholar 

  87. Altincicek B, Elashry A, Guz N, Grundler FM, Vilcinskas A, Dehne HW. Next generation sequencing based transcriptome analysis of septic-injury responsive genes in the beetle Tribolium castaneum. PLoS ONE. 2013;8(1):e52004.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Theopold U, Li D, Fabbri M, Scherfer C, Schmidt O. The coagulation of insect hemolymph. Cell Molec Life Sci. 2002;59(2):363–72.

    Article  CAS  PubMed  Google Scholar 

  89. Theopold U, Schmidt O, Söderhäll K, Dushay MS. Coagulation in arthropods: defence, wound closure and healing. Trends Immunol. 2004;25(6):289–94.

    Article  CAS  PubMed  Google Scholar 

  90. Kotani E, Yamakawa M, Iwamoto S, Tashiro M, Mori H, Sumida M, et al. Cloning and expression of the gene of hemocytin, an insect humoral lectin which is homologous with the mammalian von Willebrand factor. Biochim Biophys Acta. 1995;1260(3):245–58.

    Article  PubMed  Google Scholar 

  91. Goto A, Kumagai T, Kumagai C, Hirose J, Narita H, Mori H, et al. (2001). A Drosophila haemocyte-specific protein, hemolectin, similar to human von Willebrand factor. Biochem J. 2001;359(1):99–108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Eleftherianos I, Revenis C. Role and importance of phenoloxidase in insect hemostasis. J Innate Immun. 2011;3(1):28–33.

    Article  CAS  PubMed  Google Scholar 

  93. Cerenius L, Kawabata SI, Lee BL, Nonaka M, Söderhäll K. Proteolytic cascades and their involvement in invertebrate immunity. Trends Biochem Sci. 2010;35(10):575–83.

    Article  CAS  PubMed  Google Scholar 

  94. Bergner A, Oganessyan V, Muta T, Iwanaga S, Typke D, Huber R, et al. Crystal structure of a coagulogen, the clotting protein from horseshoe crab: a structural homologue of nerve growth factor. EMBO J. 1996;15(24):6789–97.

    Article  PubMed  PubMed Central  Google Scholar 

  95. Valanne S, Wang JH, Rämet M. The Drosophila toll signaling pathway. J Immunol. 2011;186(2):649–56.

    Article  CAS  PubMed  Google Scholar 

  96. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. The extracellular matrix of animals. In: molecular biology of the cell, 4th edition. New York: Garland Science; 2002. Available at

    Google Scholar 

  97. Lombardo F, Christophides GK. Novel factors of Anopheles gambiae haemocyte immune response to plasmodium berghei infection. Parasit Vectors. 2016;9(1):78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Bryan J, Edwards R, Matsudaira P, Otto J, Wulfkuhle J. Fascin, an echinoid actin-bundling protein, is a homolog of the Drosophila singed gene product. Proc Natl Acad Sci U S A. 1993;90(19):9115–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Sokol NS, Cooley L. Drosophila filamin encoded by the cheerio locus is a component of ovarian ring canals. Curr Biol. 1999;9(21):1221–30.

    Article  CAS  PubMed  Google Scholar 

  100. Te Boekhorst V, Preziosi L, Friedl P. Plasticity of cell migration in vivo and in silico. Annu Rev Cell Dev Biol. 2016;32(1):491–526.

    Article  CAS  Google Scholar 

  101. Huff T, Otto AM, Müller CS, Meier M, Hannappel E. Thymosin β4 is released from human blood platelets and attached by factor XIIIa (transglutaminase) to fibrin and collagen. FASEB J. 2002;16(7):691–6.

    Article  CAS  PubMed  Google Scholar 

  102. Jing J, Sun X, Zhou C, Zhang Y, Shen Y, Zeng X, et al. Cloning, expression and effects of P. americana thymosin on wound healing. Int J Mol Sci. 2019;20(19):4932.

    Article  CAS  PubMed Central  Google Scholar 

  103. Saelee N, Noonin C, Nupan B, Junkunlo K, Phongdara A, Lin X, et al. β-thymosins and hemocyte homeostasis in a crustacean. PloS ONE. 2013;8(4):e60974.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. Smith VJ, Accorsi A, Malagoli D. Hematopoiesis and hemocytes in pancrustacean and molluscan models. In: Malagoli D, editor. The evolution of the immune system. Cambridge: Academic Press; 2016. p. 1–28.

    Chapter  Google Scholar 

  105. Hillyer JF. Insect immunology and hematopoiesis. Dev Comp Immunol. 2016;58:102–18.

    Article  CAS  PubMed  Google Scholar 

  106. Prestwich GD. Phosphoinositide signaling: from affinity probes to pharmaceutical targets. Chem Biol. 2004;11(5):619–37.

    Article  CAS  PubMed  Google Scholar 

  107. Xie S, Naslavsky N, Caplan S. Diacylglycerol kinases in membrane trafficking. Cell Logist. 2015;5(2):e1078431.

    Article  PubMed  PubMed Central  Google Scholar 

  108. Dong X, Li X, Li Q, Jia H, Zhang H. The inducible blockage of RNAi reveals a role for polyunsaturated fatty acids in the regulation of dsRNA-endocytic capacity in Bactrocera dorsalis. Sci Rep. 2017;7(1):1–11.

    Article  CAS  Google Scholar 

  109. Lehmann M. Diverse roles of phosphatidate phosphatases in insect development and metabolism. Insect Biochem Mol Biol. 2020.

  110. Kutateladze TG. Translation of the phosphoinositide code by PI effectors. Nat Chem Biol. 2010;6(7):507–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Gallardo R, Ivarsson Y, Schymkowitz J, Rousseau F, Zimmermann P. Structural diversity of PDZ–lipid interactions. Chem Biochem. 2010;11(4):456–67.

    Article  CAS  Google Scholar 

  112. Balla T. Phosphoinositides: tiny lipids with giant impact on cell regulation. Physiol Rev. 2013;93(3):1019–137.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  113. Van Buul JD, Geerts D, Huveneers S. Rho GAPs and GEFs: controlling switches in endothelial cell adhesion. Cell Adhes Migr. 2014;8(2):108–24.

    Article  Google Scholar 

  114. Rizo J, Südhof TC. C2-domains, structure and function of a universal Ca2+−binding domain. J Biol Chem. 1998;273(26):15879–82.

    Article  CAS  PubMed  Google Scholar 

  115. Aravamudan B, Fergestad T, Davis WS, Rodesch CK, Broadie K. Drosophila UNC-13 is essential for synaptic transmission. Nat Neurosci. 1999;2(11):965–71.

    Article  CAS  PubMed  Google Scholar 

  116. Dubuque SH, Schachtner J, Nighorn AJ, Menon KP, Zinn K, Tolbert LP. Immunolocalization of synaptotagmin for the study of synapses in the developing antennal lobe of Manduca sexta. J Comp Neurol. 2001;441(4):277–87.

    Article  CAS  PubMed  Google Scholar 

  117. Tunstall NE, Herr A, de Bruyne M, Warr CG. A screen for genes expressed in the olfactory organs of Drosophila melanogaster identifies genes involved in olfactory behaviour. PLoS ONE. 2012;7(4):e35641.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  118. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  119. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  120. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.

    Article  CAS  PubMed  Google Scholar 

  121. Meccariello A, Salvemini M, Primo P, Hall B, Koskinioti P, Dalíková M, et al. Maleness-on-the-Y (MoY) orchestrates male sex determination in major agricultural fruit fly pests. Science. 2019;365(6460):1457–60.

    Article  CAS  PubMed  Google Scholar 

  122. Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and Phylogenomics. Mol Biol Evol. 2018;35(3):543–8.

    Article  CAS  PubMed  Google Scholar 

  123. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011;12(1):323.

    Article  CAS  Google Scholar 

  124. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  125. Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  126. Musacchia F, Basu S, Petrosino G, Salvemini M, Sanges R. Annocript: a flexible pipeline for the annotation of transcriptomes able to identify putative long noncoding RNAs. Bioinformatics. 2015;31(13):2199–201.

    Article  CAS  PubMed  Google Scholar 

  127. Wu C, Chakrabarty S, Jin M, Liu K, Xiao Y. Insect ATP-binding cassette (ABC) transporters: roles in xenobiotic detoxification and Bt insecticidal activity. Int J Mol Sci. 2019;20(11):2829.

    Article  CAS  PubMed Central  Google Scholar 

  128. Quevillon E, Silventoinen V, Pillai S, Harte N, Mulder N, Apweiler R, et al. InterProScan: protein domains identifier. Nucleic Acids Res. 2005;1(Web Server):33–W120.

    Article  CAS  Google Scholar 

  129. Söding J, Biegert A, Lupas AN. The HHpred interactive server for protein homology detection and structure prediction. Nucleic Acids Res. 2005;33(Web Server):244–W248.

    Article  CAS  Google Scholar 

  130. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-[Delta][Delta] CT method. Methods. 2001;25(4):402–8.

    Article  CAS  PubMed  Google Scholar 

Download references


Thanks to all students and colleagues helping collecting specimens of blister beetles on field, and R. Poloni (Modena) for providing individuals of Mylabris variabilis from Abruzzo (Italy). We wish to thank A. Riccieri (Roma Tre University) for providing starting material to set-up RT-qPCR and M.V. Modica for her precious advice on data analysis workflow. We are sincerely grateful to Prof. W.D. Woggon (University of Basel) for sharing data and opinion on cantharidin biosynthesis and blister beetle defensive behaviour. We are also grateful to three anonymous reviewers, whose suggestions greatly improved the manuscript.


This work was carried out in the frame of the project “CanBBio - Cantaridina: dalla Biodiversità alle Biotecnologie” (prot. n° 85–2017-14967) financed by Regione Lazio (Progetti Gruppi di ricerca – Conoscenza e Cooperazione per un Nuovo Modello di Sviluppo, L.R. Lazio 13/08) and co-funded by MIUR-Italy (Departments of Excellence, L. 232/2016, art.1: 314–337, awarded to Dept. of Science, University of Roma Tre, Rome Italy for 2018–2022).

Author information

Authors and Affiliations



EF carried out field collections, tissue dissections, molecular lab procedures (RNA extraction, RT-qPCR), blast searches on selected gene families and participated to drafting tables, figures and text of the manuscript; MS performed bioinformatic analyses (i.e. raw data processing, transcriptome de novo assembly, differential expression analyses) and contributed to drafting the manuscript; FL helped in setting up the molecular lab workflow, performed statistical analyses for RT-qPCR validation, helped in writing and reviewed the manuscript; MMu helped to carry out field collections, reared the specimens, dissected tissues for RNA extraction and reviewed the manuscript; MMo helped in field collections, performed blast searches on selected gene families and participated to drafting tables; SG contributed to set-up lab procedures for haemolymph collection and reviewed the manuscript; TP and VD helped in discussing data on detoxification strategies and reviewed the manuscript; TG and ER helped in discussing data on cantharidin biosynthesis and reviewed the manuscript; PM helped in setting up molecular lab procedures and reviewed the manuscript; ADG participated to field collections, supervised tissue dissections and reviewed the manuscript; MAB supervised and coordinated field collections, performed specimen identification, helped in conceiving the study, writing and reviewing the manuscript; MC supervised and coordinated molecular lab activities (RNA extraction, RT-qPCR), participated to field collections and drafted the manuscript; EM conceived, designed and coordinated the study, participated to field collections, performed statistical analyses, prepared tables and figures and wrote the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Manuela Cervelli or Emiliano Mancini.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The author declares no bias nor conflicts of interest in this research.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1 Table S1-S2.

Annotation and sequence data of LT and MV longest transcript for each gene cluster.

Additional file 2 Table S3-S6.

Annotation of upregulated transcripts in MRT and haemolymph (as compared to whole body) of LT and MV (FDR < 0.001 –logFC> 2).

Additional file 3 Table S7-S10.

Complete list of significantly (adj. p < 0.05) over-represented Gene Ontology (GO) terms in differentially expressed transcripts of MRT and haemolymph of LT and MV.

Additional file 4 Table S11-S14.

Complete list of significantly (adj. p < 0.05) upregulated protein families (Pfam) in MRT and haemolymph of LT and MV.

Additional file 5 Table S15-S16.

Expression data (FPKM) and complete list of transcripts putatively involved in CA biosynthesis pathways and sequestration/self-detoxification mechanisms in LT and MV.

Additional file 6 Table S17.

List of primers used for RT-qPCR validation.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fratini, E., Salvemini, M., Lombardo, F. et al. Unraveling the role of male reproductive tract and haemolymph in cantharidin-exuding Lydus trimaculatus and Mylabris variabilis (Coleoptera: Meloidae): a comparative transcriptomics approach. BMC Genomics 22, 808 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: