Odoriferous Defensive Stink Gland Transcriptome to Identify Novel Genes Necessary for Quinone Synthesis in the Red Flour Beetle,

Summary

Chemical defense is one of the most important traits, which endow insects the ability to conquer a most diverse set of ecological environments. Chemical secretions are used for defense against anything from vertebrate or invertebrate predators to prokaryotic or eukaryotic parasites or food competitors. Tenebrionid beetles are especially prolific in this category, producing several varieties of substituted benzoquinone compounds. In order to get a better understanding of the genetic and molecular basis of defensive secretions, we performed RNA sequencing in a newly emerging insect model, the red flour beetle Tribolium castaneum (Coleoptera: Tenebrionidae). To detect genes that are highly and specifically expressed in the odoriferous gland tissues that secret defensive chemical compounds, we compared them to a control tissue, the anterior abdomen. 511 genes were identified in different subtraction groups. Of these, 77 genes were functionally analyzed by RNA interference (RNAi) to recognize induced gland alterations morphologically or changes in gland volatiles by gas chromatography-mass spectrometry. 29 genes (38%) presented strong visible phenotypes, while 67 genes (87%) showed alterations of at least one gland content. Three of these genes showing quinone-less (ql) phenotypes – Tcas-ql VTGl; Tcas-ql ARSB; Tcas-ql MRP – were isolated, molecularly characterized, their expression identified in both types of the secretory glandular cells, and their function determined by quantification of all main components after RNAi. In addition, microbe inhibition assays revealed that a quinone-free status is unable to impede bacterial or fungal growth. Phylogenetic analyses of these three genes indicate that they have evolved independently and specifically for chemical defense in beetles.

Introduction

Insects are among the most diverse group of animals on the planet and amazingly include more than a million described species, which is more than half of all known living organisms [1], [2]. Moreover, they have conquered almost every environment on earth. From a series of distinctive attributes that orchestrate together to endow them with the ability to live in a wide range of ecological environments, chemical defense is one of the most important traits [3]. Many chemical secretions have repellent or irritant properties [4], [5]. Tenebrionid beetles are especially prolific by producing several various substituted benzoquinone compounds [5]–[9]. Tribolium beetles (Coleoptera: Tenebrionidae) have dragged attentions of researchers to their particular secretions, since it was noted that their flour medium turns pink over time due to the secretion of a gaseous substance from adults [10], which also deleteriously affects the viscous and elastic properties of dough made from such infested flour [11].

Tribolium beetles possess for the purpose of chemical defense two pairs of specialized secretory organs – one in the prothorax and one in the posterior abdomen – termed odoriferous or stink glands [12]. The glands located in the prothorax are called prothoracic, thoracic, or anterior glands, while the other pair in the abdomen is referred to as abdominal, posterior, or pygidial glands. The fine structure of these glands revealed two types of secretory units composed of two slightly different types of cells with particular vesicular organelles (cell type 1 and cell type 2), tubules, reservoir, ducts and muscles [12]–[14]. At least four members of the genus Tribolium (T. anaphe, T. castaneum, T. confusum, and T. destructor) use the glands to produce the quinone derivatives 2-methoxybenzoquinone, ethyl-1,4-benzoquinone (EBQ), and methyl-1,4-benzoquinone (MBQ) [9], [12], [14]–[18]. However, only the latter two substances were detected in T. confusum [4], [6]. Besides benzoquinone derivatives, hydrocarbons were also reported as major secretion components. T. confusum secrets 1-pentadecene [19], 1,6-pentadecadiene and smaller amounts of 1-hexadecene, 1,6-hexadecadiene, hexadecatriene, 1-heptadecene, 1,8-heptadecadiene and heptadecatriene [20], [21]. In T. castaneum, 1-pentadecene and 1,6-pentadecadiene were identified, with the former as the main component [9], [18], plus two still unidentified hydrocarbons were found [18].

The red flour beetle, T. castaneum has been developed into a highly sophisticated genetic model organism [22] with plenty of genetic and genomic tools: reverse genetics based on systemic RNA interference [23], [24], forward genetics based on insertional mutagenesis [25], [26], transgene-based mis-expression systems [27], [28], as well as a fully annotated genome sequence [29]. Moreover, several mutants with odoriferous gland phenotypes, such as melanotic stink glands (msg, with both pairs of glands melanized) [30], tar (only prothoracic glands are darkly pigmented), and box (A^box, similar to tar, but only the abdominal glands are affected) [31] are known. For quinone biosynthesis, it has been reported in another quinone-producing tenebrionid beetle Eleodes longicollis that alkylated benzoquinones are formed by acetate condensation whereas the p-benzoquinones are generated from preformed aromatic rings of amino acids (tyrosine and phenylalanine) [5], [32], [33]. In the glandular secretory cells, p-quinones are in a form of phenolic ß-glucoside contained in the more apical part, which were then transferred to inner part of the gland and form active quinones by a series of enzymatic reactions [14]. The alkenes are biosynthesized from fatty acids by oxidative decarboxylation [21]. However, no data are available on the genes involved in these processes. In addition, arthropods have evolved mechanisms to reduce the autointoxicative effects of a variety of toxic compounds produced by themselves for defense purpose [5]. For example, a paradoxsomatid millipede, Oxidus gracilis, which produces toxic phenol from tyrosine, possesses the ability to rapidly convert the exogenous phenol to tyrosine by tyrosine phenol lyase [34]. Benzoquinones are highly reactive, unstable and also toxic. Besides the purpose of chemical defense, Tribolium beetles and other insects use them also as tanning agent and to sclerotize cuticles [35]–[37], which requires perfect handling and detoxication systems. Thus understanding of the mechanisms involved in the autodetoxication of the defensive compounds might provide inspirations to manage this cosmopolitan pest and potentially other coleopteran pests. Obviously, tenebrionids are protected from their own toxic secretions by cuticular linings both internally and externally [5]. Tribolium beetles have the ability to partition the secretion away from the somatic cells, firstly by producing the secretions in the cuticle-lined organelles [14] and then keeping them in storage sacs (reservoirs) that are formed from invaginations of the cuticle [12]. The newly emerged Tribolium adults lack the defensive secretions, implying the need for building up an adequate self-protective barrier [8]. Consequently, if such self-protection systems could be broken, the pests will harm themselves. Here, we present the first transcriptome data on beetle stink glands, which allowed us to identify three genes necessary for the production of defensive quinones.

Results

Stink gland transcriptome sequencing

The odoriferous defensive stink glands were dissected and identified by their special morphological structures, vesicular organelles [13], [14] (Figure S1). mRNA sequencing was performed in six stink gland samples and one control sample (anterior abdomen, where no odoriferous glands are located) on a next generation sequencing (NGS) platform. The abbreviations for the samples were, s1: sample 1, anterior abdomen from wild-type; s2: prothoracic glands from tar mutant; s3: wild-type male prothoracic glands; s4: wild-type female prothoracic glands; s5: wild-type male abdominal glands; s6: wild-type female abdominal glands. After sequencing, 27.8 to 29.7 million reads were obtained from each sample. About 50% of them were successfully mapped to the reference database, i.e. mRNAs of the Tribolium official gene set (OGS) in Beetlebase [22], [38]. And the average depths/coverages were in the range of 23.5 to 27.8. Moreover, ratios of covered region in reference varied from 52.5% to 68.0%. Detailed statistics are presented in Table 1. The mapping of the different transcriptomics samples as individual tracks to the newest Tribolium genome assembly is provided by the public iBeetle Genome Browser (http://bioinf.uni-greifswald.de/gb2/gbrowse/tcas).

mRNA-seq library subtractions

The constructed mRNA-seq libraries are presented in Dataset S1, which shows the reads and coverage (depth) of each gene in all the samples, respectively. Table 1 indicates that the relative total reads of all samples were quite close to each other. Therefore, the read number represented actual expression levels of all the Tribolium genes in various samples. In order to screen out the differentially or specifically expressed genes in odoriferous glands, ten different groups of subtractions were performed and in total 511 genes were identified (Figure 1; for detailed list and the subtractive conditions, see Dataset S2). There were 62 genes in Group 1, standing for the gland-specific genes, which were highly expressed in all wild-type gland samples [s3, s4, s5, s6] but not in the anterior abdomen control [s1]. Group 2 presented 23 thoracic gland-specific genes, which had higher reads in wild-type thoracic glands [s3 and s4] than in wild-type abdominal glands [s5 and s6 (and s1 control)]. Group 3 had 40 abdominal gland-specific genes [s5 and s6 against s3, s4 (and s1)]. For sex related subtractions, Group 4 offered zero male stink gland-specific genes [s3 and s5 against s4, s6 (and s1)], but Group 5 had four female stink gland-specific genes [s4 and s6 against s3 and s5]. There were zero male thoracic gland-specific genes in Group 6 [s3 against s4 and s1], yet 299 male abdominal gland-specific genes in Group 7 [s5 against s6 and s1]. Meanwhile, there were also zero female thoracic gland-specific genes in Group 8 [s4 against s3 and s1], but 39 female abdominal gland-specific genes in Group 9 [s6 against s5 and s1]. And Group 10 presented 44 genes which were either up or down regulated in prothoracic glands of tar mutants. The high number of genes identified in Group 7 is probably due to the contamination of the abdominal gland sample by male accessory glands, which are hard to separate from the abdominal glands.

Gene ontology annotation

Gene ontology (GO) annotation allows meta-analyses of gene populations and associates the targeted genes to specific terms with hierarchical vocabularies describing three independent ontologies: biological process, molecular function, and cellular component [39]. Analyses were performed with 1451 genes abundant in the control, 1206 genes abundant in wild-type glands, and the 511 genes from subtraction Group 1 to Group 10 (The genes are listed in Dataset S2). Results (Figure 2A) showed that many genes were classified to metabolic and cellular processes in the GO term of biological process, and to catalytic activity and binding in molecular function. For the cellular component, most genes belonged to cell, macromolecular complex and organelle. These implied the existence of strong metabolisms in both gland and anterior abdomen. Moreover, similar trends were observed in separated analyses in the different subtraction groups (Figure 2B). Detailed GO results are presented in Dataset S3.

Transcriptomic exploration of candidate genes for quinone synthesis

Glucosidases, phenol oxidases, and peroxidases have been considered to be involved in the production of quinones in the odoriferous glands [14] and were annotated in the Tribolium genome [29]. In our stink gland transcriptome analysis, we have now explored these candidate genes for expression at the gland transcriptome level (Figure 3, details in Dataset S4). In total, 19 glucosidase, 14 phenol oxidase, and 18 peroxidase encoding genes were identified through blast searches and conserved domain confirmation. Transcriptomic explorations revealed that at least four glucosidase (TC000223, TC000537, TC002741, and TC011354), five phenol oxidase (TC000821, TC005376, TC006769, TC010489, and TC10490) and four peroxidase (TC010355, TC010362, TC012328 and TC014929) genes have increased reads in the gland samples, which confirms the importance of these three types of enzymes in defensive secretions and verifies the reliability of our transcriptome data.

Functional analysis of the most highly and gland-specifically expressed genes

In order to find novel gene functions involved in quinone synthesis, we functionally analyzed 77 genes from transcriptomic subtraction groups 1, 2, and 10 that were at least 64× higher expressed in the glands compared to the control tissue. RNAi of these genes resulted in various abnormal visible phenotypes (Figure 4). Additionally, gas chromatography–mass spectrometry (GC-MS) measurements revealed the alterations of different chemical components in both pairs of glands (an example of the chromatogram is depicted in Figure S2). The main components identified are listed in Table 2. Based on the extents of alterations of the chemicals, phenotypes were classified into six strengths (Figure 5A). 29 genes (38%, strength 1–3) showed strong changes, i.e. more than 75% reduction of at least one component, which were mostly accompanied with visible phenotypes (Figure 5B). In total, 67 of 77 genes (87%) showed alterations of at least one secreted chemical. Detailed descriptions on the phenotypic changes of all the 77 genes can be found in Dataset S5. In addition, gland cellular morphology was explored in all the 29 genes with strong phenotypes, but no visible abnormalities were observed in the secretory cells (data not shown).

Quantification of volatile gland contents

Previous research has revealed the amount of different glandular components only on the whole beetle level [8], [9], [16], [18], [40]–[43]. In order to elucidate the chemical compositions of the volatiles in the different pairs of glands and the extent of reduction after gene knock-downs, three genes with strong quinone-less phenotypes were chosen from the 77 tested genes to quantify different glandular components. Wild-type and Enhanced green fluorescent protein (EGFP) double stranded RNA (dsRNA)-injected beetles were used as controls. Figure 6 shows the complete losses of all quinones in both pairs of glands in both females and males from the knock-down beetles of these three quinone-less genes (except for one GT39 dsRNA injected male out of sixteen injected males). In comparison, the alkenes were reduced to different extents. Statistical analyses revealed significant differences between wild-type and the knock-downs (Dataset S6), EGFP dsRNA injection surprisingly caused a few significant differences in alkenes compared to wild-type but not in quinones. Interestingly, all the alkenes in prothoracic glands of GT63 knock-downs were not statistically different from the wild-type, while only heptadecadiene and heptadecene in the abdominal glands showed the same trend. Sex differences were also analyzed (Dataset S6), which showed that most chemicals had no significant differences between males and females, except for all the alkenes in abdominal glands of GT39 knock-downs and in thoracic glands of GT62. In wild-type, only heptadecadiene showed a significant difference between different sexes while the other chemicals did not.

Additionally, the amount of all the main components in wild-type A10 beetles is presented in Table 3 (details in Dataset S6). The prothoracic glands possess about 40% of either quinones or alkenes of all the stored secretions in the whole beetle, while abdominal glands have about 60%. But the molar ratios of quinones to alkenes are almost the same in both pairs of glands (thr, 2.60; abd, 2.70–2.74). And the molar ratios of MBQ to EBQ vary from 0.77 to 0.88 in different gland and sex levels. The only major dissimilarity between those two glands is the composition of distinct alkenes. The prothoracic glands have higher portions of 1-heptadecene (17ene) and 1,8-heptadecadiene (17diene), especially the former, but a lower portion of 1-pentadecene (15ene) (15ene: 17diene: 17ene = ∼60%: 28%: 12% in thr, ∼88%: 4%: 8% in abd).

Full length cDNAs of the three quinone-less genes, GT39, GT62 and GT63, were cloned and the sequences submitted to GenBank (accession numbers JX569829, JX569830, and JX569831, respectively). Based on the phenotypes and their homology (see Dataset S7, S8, and S9 for respective protein sequences), GT39 has been designated as Tcas-quinone-less vitellogenin-like (Tcas-ql VTGl), GT62 as Tcas-quinone-less arylsulfatase b (Tcas-ql ARSB) and GT63 as Tcas-quinone-less multi-drug resistance protein (Tcas-ql MRP).

Phylogeny of the three newly identified quinone-less proteins

After quantification of the gland volatiles in these three gene knock-downs, their phylogeny was explored. In the phylogenetic tree (Figure 7A), Tcas-ql VTGl (GT39) was classified together with eight other Tribolium homologs, and 11 homologs were grouped in a branch close by, including proteins from Strongylocentrotus purpuratus, Danio rerio, Gallus gallus, Homo sapiens, Mus musculus, and Nasonia vitripennis (see Dataset S7 for all the sequences). Tcas-ql ARSB (GT62) was grouped together with two other Tribolium homologs (Figure 7B, see Dataset S8 for all the sequences), which were closest related to three Nasonia proteins, similarly to Tcas-ql MRP (GT63, Figure 7C, see Dataset S9 for all the sequences). Since all three quinone-less gene encoded proteins had several Tribolium homologs, we checked their homologs expression levels in the gland tanscriptome. The homologs were linked with corresponding GLEAN predictions and explored in the transcriptomic libraries (Dataset S10). It was shown (Figure 8) that no other gene was as highly expressed in the gland samples except for the three identified quinone-less genes, although some genes contained more than 16 times (2 to the power of 4, GI:91085475, Figure 8C) or 4 times (2 to the power of 2, GI:189236319, Figure 8B) higher reads in some cases. In conclusion, this indicates that all three quinone-less genes most probably have evolved independently and specifically for the quinone-based chemical defensive system in Tribolium.

Expression patterns of the quinone-less genes in gland tissue

The expression patterns of Tcas-ql VTGl (GT39), Tcas-ql ARSB (GT62) and Tcas-ql MRP (GT63) were explored in dissected odoriferous glands by gland whole mount fluorescent in situ hybridization (GWMFISH). All three genes showed strong expressions in both gland cell type 1 and cell type 2 [13], [14] of both pairs of glands (Figure 9), which confirms their involvement in Tribolium defensive secretion.

Microbe inhibition assays

Microbe inhibition tests were carried out to identify the effect of losing quinones. It had been shown that beetle chemical secretions – especially the benzoquinones – can inhibit the growth of several microbes common to flour with artificial MBQ having the same effect [44], [45]. In our experiments, a fungus, Aspergillus niger, and a gram positive bacterium, Arthrobacter globiformis, were used to test the strength of the chemical defense mediated by beetle glands. As a control, we used RNAi against GT12 (detailed analysis on this gene will be presented elsewhere), representing a gene causing an alkene-less phenotype at knock-down, which is accompanied by only a slight reduction of quinones (see Dataset S5). The results (Figure 10) showed that wild-type beetle secretions could inhibit microbe growth in a certain area, but loss of quinones led to undetectable inhibitions (GT63), while an alkene-less state had only reduced inhibition effects (GT12).

Phenol oxidase activity assays

Since chemical defense systems are responsible for defending the host from infection, we wanted to test, whether also a part of the innate immune system is affected directly or indirectly by the quinone-less phenotype. Thus, after RNAi-mediated knock-down of the quinone-less genes, phenol oxidase (PO) activities were measured as a general index of the melanization innate immune responses in invertebrates [46]. Compared to wild-type beetles, the three quinone-less gene knock-downs had significantly reduced levels of PO-activity both in females (Figure 11) and males (Figure S3), while buffer or EGFP dsRNA injected animals did not show significant changes. These data indicate that the extra-corporal chemical defense may be linked with the function of a part of the innate immune system in Tribolium.

Discussion

Transcriptome sequencing was performed in odoriferous gland samples from Tribolium castaneum to identify novel gene functions in defensive quinone synthesis. In a next generation sequencing approach, 27.8–29.7 million reads were obtained, of which only 50 to 61% reads were successfully mapped to the Tribolium genome. Part of the low mapping ratio is probably explainable by the strain differences used for genome sequencing (inbred Georgia 2 strain) [29] and the San Bernardino strain we used for transcriptomics together with the strict mapping parameters. However, this did neither affect the comparison and subtraction between different samples, nor the library subtraction results, since all our wild-type samples have the same strain background and all samples were treated identically. Only the tar mutation is in different genetic background derived from sooty [47] and Maxillopedia-Dachs3 strains (mxp^Dachs3) [31]. However, the mapping ratio for this sample was similar to the others. Another reason for the low mapping ratio might be that the official gene set was from a combination of different automated gene prediction programs and annotation pipelines [29], which makes it possible that there were still unidentified genes in the Tribolium genome.

The comparisons we performed were between different tissues in wild-type and tar mutant. During the library subtractions, we chose a general cut-off of fold change at 64 times, which is much higher compared to many microarray analyses (mostly two times), but the number of the genes we got was reasonably high and seemed suitable to start to work with. Only subtraction Group 7 (male abdominal gland-specific genes), with 299 genes (Figure 1), showed an unusual high number, which might, however, be explained by the fact that the male accessory glands are hard to dissect away from the abdominal glands and many of the genes in Group 7 might actually be male accessory gland-specifically expressed genes. As these were not a topic of this study, these potentially interesting genes remain for future analysis.

GO annotations revealed that the glands have quite active metabolisms with many catalytic and binding related genes being highly expressed. The GO annotation rate of the genes, that had coverage of more than 50 in all the glands, was 70.1%, while the control had 78.8% (Dataset S3). This suggests that there are more orphan genes expressed in Tribolium odoriferous glands than in the control. In addition, only 53.6% of all the 511 genes from the subtractions were annotated, and surprisingly, Group 7 (male specific glands genes) had an annotation rate of only 42.5%, suggesting an even higher number of genes with unknown functions. We identified also some glucosidases, phenol oxidases and peroxidases highly transcribed in the glands, that are candidate enzymes to be involved in quinone biosynthesis. In conclusion, our transcriptome data have reliably detected candidate genes involved in quinone biosynthetic mechanisms of chemical defense in the red flour beetle.

Exploring the functions of a first batch of highly gland-specifically expressed genes by RNAi and GC-MS to potentially identify novel gene functions in quinone biosynthesis, 67 of 77 genes (87%) showed alterations of at least one secreted chemical, which not only confirmed their importance for semiochemical synthesis, but also signified the effectiveness of our transcriptome screening. Surprisingly, some genes with very high reads in the glands showed no big changes at the chemical level. For example, GT26 (TC007317), GT35 (TC010551) and GT41 (TC011337) from Group 1 had more than 441, 514 and 690 times higher reads respectively in all wild-type gland samples than control (Dataset S1), however, their knock-downs showed less than 75% reductions, or even neglectable changes (Dataset S5). We suggest that these genes might be involved in other biological processes indirectly related to chemical secretion. Additionally, in Group 10, the GT23 gene (TC006131) that had more than 126 times enriched reads in s2-tthr (tar prothoracic gland sample) compared to s3-mthr and s4-fthr (male and female prothoracic gland samples; Dataset S1), showed no changes in prothoracic glands, but slightly increased amounts of quinones and reduced amounts of alkenes in abdominal glands. Encoding the odorant binding protein 21 (OBP21, GI:270012767), GT23 might be involved in olfaction system. It is possible that the mutated tar somehow caused the mis-expression of this gene in a different type of tissue, or OBP21 belongs to the ubiquitous OBP type, such as encapsulin [48], which is probably involved in diverse physiological functions [49]. Moreover, GT23 showed about 10 times more reads in wild-type prothoracic glands than in the control sample. In addition, another OBP (GT76) and two chemosensory protein (GT30, GT77) encoding genes showed remarkable expression changes in the prothoracic glands of wild-type compared to tar mutants in the opposite direction. GT30, GT76, and GT77 are expressed at high levels in wild-type prothoracic glands, but their expression is strongly reduced in tar mutants. None of these four genes is expressed at significant levels in the abdominal glands, indicating a specific function for the anterior glands. However, no significant changes in volatile gland contents could be detected after RNAi knock-down of those genes.

During the morphological analyses, many abnormal glands were observed (Figure 4). However, besides the quantitative changes in the main glandular volatiles, GC-MS of the referred gene knock-downs showed no additional peak compared to wild-type, which indicates non-volatility of the accumulated brownish or black substances. The gland phenotypes could be explained as the knock-downs triggered the inhibition of chemical syntheses or the blocking of their transportation, or in some cases, the accumulation of intermediate substrates or unknown brown or black polymers. It was proposed for previously identified msg mutants that the black material is of high-molecular-weight and polymeric consisting of polymerized prematurely formed quinones due to the absence of the inhibitor in oxidation of hydroquinone [30], [50]. Moreover, the lack of specific amphipathic molecules – such as terminal alkenes – might cause the various phases (probably organic and aqueous) seen e.g. in Figure 4 G4, G5.

In the quantification part, the main glandular contents were quantified separately in both pairs of glands for the first time. Assuming that stage A10 at 32.5°C is equal to A12 at 30°C, which was predicted based on the Tribolium life parameter table [51], the MBQ and EBQ amounts we got (Dataset S6) in wild-type were 20–30% (females) and 40–70% (males) higher than the amounts detected by Unruh et al. [8]. This indicates that the dissection based extraction is much more accurate than the homogenization based method, since the latter may cause the loss of unstable chemicals during the crude preparation. Moreover, in our experiments, the males and females were not separated before harvesting, which is much closer to the natural conditions compared to the method Unruh et al. used [8]. In addition, more EBQ was detected in our tests (molar ratio of MBQ/EBQ: 0.81 in female, 0.87 in male), while the previously reported ratio was in the range of 0.59–0.61 [18], [42]. However, the weight ratio of quinones in the whole beetle (61.5%) was only a bit higher than 58.3% reported previously [18]. Interestingly, different hydrocarbon compositions were observed between the prothoracic and abdominal glands, which might reflect the dissimilar usage of their precursors, fatty acids [52], in distinct body parts and sexes. Furthermore, except for heptadecadiene in abdominal glands, all other components presented no significant differences between male and female at stage A10. Therefore, we propose that both sexes possess similar secretion levels in normal environment, and the higher level of benzoquinones in female observed before [8] was due to a different energy allocation when they produce no or less eggs as virgins, since reproduction (mating and egg production) could change the energy allocation and fitness in several other species [53]–[57].

After these first functional analyses of the complete set of highly gland-specifically expressed genes, three genes causing a quinone-less phenotype at knock-down were characterized further in our work. The quantitative data clearly showed the loss of quinones and the reduction of alkenes in both pairs of glands. Surprisingly, all the alkenes in Tcas-ql MRP (GT63) knock-downs showed no significant differences to wild-type in prothoracic glands but were significantly different in abdominal glands, while Tcas-ql VTGl (GT39) knock-downs caused significant differences to wild-type in almost all the comparisons except the three alkenes in the male abdominal glands. Tcas-ql ARSB (GT62) knock-downs showed significant changes to wild-type in all comparisons (Dataset S6). This implies that Tcas-ql MRP (GT63) has distinct functions in different glands, while Tcas-ql VTGl (GT39) and Tcas-ql ARSB (GT62) function diversely in both glands of the distinct sexes.

Potential molecular functions of these three genes in quinone synthesis are predicted on account of the homology and conserved domain analyses. Possessing a pancreatic lipase-like enzyme domain, Tcas-ql VTGl (GT39) is a vitellogenin like protein (31% identity). Vitellogenin is classified as a glycolipoprotein, having properties of a sugar, fat and protein, and belongs to a lipid transport protein family [58]. Some data in mealworm showed that a vitellogenin like protein (19.6% identity) could enhance the melanin synthetic process [59], in which o-quinones were produced in an intermediate step right before synthesizing eumelanin [60], [61]. Possibly, the benzoquinone production in odoriferous glands has similar pathways, in which Tcas-ql VTGl (GT39) regulates the related enzyme activity or reactions yet with more important roles, since almost no quinones were detected in its knock-down. Provided that common synthetic steps exist, the black material in msg and tar mutants might actually be composed of melanin or something alike.

Tcas-ql ARSB (GT62) is an arylsulfatase B (ARSB) protein (45% identity), which has a sulfatase domain responsible to hydrolyze sulfates in the body by breaking down large sugar molecules called glycosaminoglycans (GAGs). ARSB targets two GAGs in particular: dermatan sulfate and chondroitin sulfate. ARSB is located in lysosomes, compartments within cells that digest and recycle different types of molecules [62]. The deficiency of ARSB is the cause of mucopolysaccharidosis VI (MPS VI) which occurs in humans and cats, called also Maroteaux-Lamy syndrome. MPS VI is a progressive condition that causes many tissues and organs to enlarge and become inflamed or scarred, skeletal abnormalities are also common in this condition [63], [64]. Since lysosomes are the waste disposal system in the cell [65], a few possible functions of Tcas-ql ARSB (GT62) are proposed. Firstly, Tcas-ql ARSB (GT62) may have important roles in the detoxication of toxic substances in gland cells, whose knock-down leads to the initiation of an assumed feedback loop, then to the inhibitions of secretory chemical syntheses. Secondly, since different subcellular localization prediction tools gave distinct results (data not shown), Tcas-ql ARSB (GT62) may not be located in lysosomes but in the cytoplasm or elsewhere, functioning as an essential transporter for the intermediates involved in both quinone and alkene production, or as a key enzyme responsible for the activation of the newly translated transporters or other vital related proteins, or simply controlling the energy source of transportation, such as the pH gradient, or ion donators.

Tcas-ql MRP (GT63) belongs to ATP-binding cassette transporters, subfamily C (ABCC), which is also known as multidrug resistance-associated protein (MRP). Depending on ATP, the members of the MRP family can transport both hydrophobic uncharged molecules and water-soluble anionic compounds [66]. The latter includes the substrates conjugated with anions, such as glutathione, glucuronate or sulfate [67]. Considering our data, Tcas-ql MRP (GT63) might be an important transporter for the quinone precursors in all gland secretory cells, which also transports some alkene precursors in abdominal glands. And the transportations may occur from the hemolymph to glandular cells, or from the secreting cells to the vesicular compartments, where the toxicant-producing reactions are being segregated to [14].

Despite of these predictions on the functions of the three quinone-less gene encoded proteins based on their homologs, it should be noted that our phylogenic analyses show that they evolved independently and particularly for the chemical defense in the red flour beetle. Moreover, based on the GWMFISH results that confirm their expressions in the secretory units, all three genes play essential roles in producing the defensive quinones in the odoriferous glands of T. castaneum.

In chemical secretion of Tribolium and other quinone-producing arthropods, the quinones and the alkenes (hydrocarbons) have their respective roles. Quinones are toxic and quite reactive, therefore mainly responsible for the defense against pathogens, parasitoids, and predators [5]. Alkenes and mixtures of other organic solvents partition quinones to the hydrocarbon phase, help to spread them on the arthropod cuticle, and aid penetration of quinones through the cuticles of various enemies to improve the defensive effects [5], [68], [69]. Our inhibition assay results confirm these descriptions, since the quinone-less gene knock-downs (Tcas-ql MRP, GT63) showed no microbe inhibitions at all, whereas the alkene-less gene knock-down (GT12) presented only reduced inhibition effects.

In insects, innate immune responses include antimicrobial peptides, phagocytosis, nodulation, melanotic encapsulation and wound healing, which endow potential hosts the abilities to defend against pathogens and parasites [70]–[74]. The chemical defense system is also responsible for defending the host from infection by other organisms. Therefore, it was interesting to check how innate immunity was affected, when the chemical defense system is knocked-down. An oxidoreductase (PO) activity has been commonly assayed to provide a general index of melanization innate immune responses in invertebrates [46], since the activation of PO leads to the melanization reactions of invading pathogens, which is a major aspect of the innate immune system [75]. Therefore, PO activities were examined in quinone-less gene knock-downs compared to wild-type, buffer injected and EGFP dsRNA injected beetles. The obtained results indicate that the chemical defense may be linked with melanotic encapsulation innate immune responses. Thus the three quinone-less genes necessary for quinone biosynthesis in odoriferous glands might also be involved in the melanization cascade. In addition, the genomic annotation and expression analyses of the Tribolium POs showed that some POs are highly expressed in all gland samples but low in control, which implies that they might mainly function in the chemical defense system, and their preferred substrates are in the glands. However, whether they are involved in the melanin pathway is so far unknown. In case the POs are multifunctional, a linkage and crosstalk between these two systems would exist.

Materials and Methods

Transcriptome sequencing

Prothoracic and abdominal odoriferous glands were dissected separately from A10–A30 (reared at 32.5°C, 10–30 days after eclosion) adult beetles and stored in RNAlater solution (Ambion, Life Technologies GmbH, Darmstadt, Germany, Cat. No. AM7020) on ice. Males and females were separately prepared except for the prothoracic glands from tar mutants. About 500 beetles were used for each gland sample, while the anterior abdomen, where no glands are located, was taken as a control tissue. Then total RNA was extracted using RNAqueous-Micro Kit (Ambion, Life Technologies GmbH, Darmstadt, Germany, Cat. No. AM1931) and treated by DNase I. Transcriptome sequencing (mRNA-seq) was performed by Macrogen Inc. (Seoul, South Korea), on a next generation sequencing (NGS) platform (Illumina/Solexa Genome Analyzer IIx). After sequencing, reads (38 bp each) were mapped to the mRNAs of the official gene set (OGS) from Beetlebase 3.0 [22], [38] by Maq tool (http://maq.sourceforge.net/). The samples (s) were, s1: anterior abdomen; s2: prothoracic glands from tar mutant; s3: male prothoracic glands; s4: female prothoracic glands; s5: male abdominal glands; s6: female abdominal glands. Except of s2, all other tissues were from wild-type. Coverage (depth) is calculated as reads times 38 divided by specific length of gene transcript. The transcriptomic sequences derived from the different samples have also been mapped as individual tracks to the Tribolium genome (iBeetle Genome Browser: http://bioinf.uni-greifswald.de/gb2/gbrowse/tcas).

Gene ontology annotation and mRNA-seq library subtractions

The genes, which had coverage over 50 (about 2 times of the whole sequencing coverage), were regarded as abundant or richly expressed in either all wild-type gland samples or control. Their functionalities were explored by gene ontology (GO) annotation [39] using Blast2go [76], [77]. In order to screen gland-specific genes, statistical subtractions were carried out among different samples for various comparisons. In general, the cut-off for logarithm of fold change, with 2 as the base, was 6, which meant 64 times more reads in one sample than the other. The detailed subtraction conditions are presented in Dataset S2.

Transcriptomic exploration of candidate genes for quinone synthesis

Tribolium glucosidase, phenol oxidase and peroxidase were searched initially in protein database at the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/protein/). The obtained proteins were characterized based on conserved domains (CDD of NCBI, http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) and probed back to the publicly accessible Tribolium genome at Beetlebase with blastp algorithm in order to be linked with OGS, avoid redundancies, and identify the homologs which were not covered by the previous searches. The newly identified proteins were then analyzed in CDD for confirmation.

Functional analysis of the most highly and gland-specifically expressed genes

To evaluate the subtraction results, 77 candidate genes were chosen from the gland transcriptome screening and functional analysis was performed by using RNA interference (RNAi) [78], [79]. An online tool, E-RNAi [80] was used to design fragments for double stranded RNA (dsRNA) synthesis with no or lowest off-target effects. Primers were designed by Primer Premier 5 [81] and listed in Dataset S5. Animals were injected with dsRNAs at mid pupal [79] or larval L5–L6 stage [82], and were checked at A10 and A24 (32.5°C) for morphological phenotypes on prothoracic and abdominal glands. Furthermore, both pairs of glands were dissected carefully and intact from one male and one female beetle and smashed in 100 µl methanol (Merck Millipore SupraSolv, Merck KGaA, Darmstadt, Germany, Cat. No. 106011), which represents a more specific gland volatile analysis than previously performed by Howard 1987 (whole beetle [17]) or Villaverde et al. 2007 (headspace [9]). Then the samples were stored at −20°C and measured within 24 hours. One microliter was loaded by a split injector into an Agilent gas chromatograph coupled with a mass spectrometer (GC-MS) (Detailed parameters are described in Text S1). The areas of the signals in chromatograms were calculated using the software MSD ChemStation D.02.00.275 (Agilent Technologies, Santa Clara, USA) under auto-integration mode. Then the data were compared between each candidate gene knock-down and the control. The phenotypes were grouped according to strengths of the alterations of the major components. For the three genes with quinone-less phenotypes, second independent dsRNA fragments, which had no overlaps with the first fragments, were designed with the same tools and used to confirm the phenotypes.

Quantification of volatile gland contents

In order to quantify different volatile components in the secretion, the following chemicals were obtained from commercial sources: methyl-1,4-benzoquinone (MBQ) (abcr GmbH & Co. KG, Karlsruhe, Germany, Cat. No. AB208176), 2-methylhydroquinone (MHQ) (abcr GmbH & Co. KG, Karlsruhe, Germany, Cat. No. AB132029), ethyl quinol (also called ethyl-1,4-hydroquinone: EHQ) (abcr GmbH & Co. KG, Karlsruhe, Germany, Cat. No. AB148997), 1-pentadecene (Fluka, SIGMA-ALDRICH Chemie GmbH, Munich, Germany, Cat. No. 76560) and 1-heptadecene (Fluka, SIGMA-ALDRICH Chemie GmbH, Munich, Germany, Cat. No. 51665). Then authentic standard solution series were made and a five-point calibration was performed by GC-MS. Based on the standard curves, the areas of the abundances from GC-MS were transformed to masses. Ethyl-1,4-benzoquinone (EBQ) and heptadecadiene, which were commercially unavailable, were calculated as equivalents based on the standard curve of EHQ and 1-heptadecene respectively. Quantification was carried out in wild-type, buffer injected, a dsRNA EGFP injected control, and three quinone-less gene knock-downs. After pupal RNAi, glands were prepared from A10 beetles (15–30 animals each sex). It was proposed that more than 80% of the glandular quinones are benzoquinones [8]. Because the small amounts of hydroquinones detected are precursors of benzoquinones [24], [64], the quantities of hydroquinones and benzoquinones were summed up and treated as secreted quinones. After quantification, statistical analyses were performed with software JMP 9.0.2 (SAS Institute, 2010) using student t-test for sex comparisons and a nonparametric method (Mann–Whitney–Wilcoxon) for group comparisons.

Phylogeny of the three novel quinone-less proteins

Full length cDNAs obtained from RACE reactions (see Text S1) were analyzed by the online tool ORF Finder (Open Reading Frame finder, http://www.ncbi.nlm.nih.gov/gorf/orfig.cgi) and translated to proteins. The amino acid sequences were submitted to NCBI to find the homologs through blastp search in Reference Proteins Database and the first fifty sequences were chosen, in which the Tribolium homologs were blasted again in Beetlebase (http://beetlebase.org/) to check redundancies and find the corresponding OGS numbers (listed in Dataset S10) in order to analyze their expressions at the glandular transcriptome level. Then all proteins (listed in Dataset S7, S8, and S9) were aligned by using MAFFT [83] and analyzed with FastTree [84] using maximum likelihood methods to construct dendrograms, which were displayed, marked and computed based on the branching frequencies (cut-off was 60%) using MEGA5 [85].

Gland whole mount fluorescent in situ hybridization

The protocol for gland whole mount fluorescent in situ hybridization (GWMFISH) was based on previous methods [86]–[90] with a few modifications. Details are described in Text S1. The only big difference was using the HNPP Fluorescent Detection Set (Roche Applied Science, Roche Diagnostics Deutschland GmbH, Mannheim, Germany, Cat. No. 11758888001) to detect the labeled probes after incubation of Anti-digoxigenin (DIG)-alkaline phosphatase (AP) Fab fragments (Roche Applied Science, Roche Diagnostics Deutschland GmbH, Mannheim, Germany, Cat. No. 11093274910; from sheep). Glands were counterstained with Hoechst 33342 (SIGMA-ALDRICH Chemie GmbH, Munich, Germany, Cat. No. B2261) prior to mounting and embedding in Aqua-Poly/Mount (Polyscience, Niles, Illinois, USA, Cat. No. 18606). The stainings were observed and captured with confocal laser scanning microscope Zeiss LSM780. 3D (3-dimensional) construction was performed using software ZEN2011 (Carl Zeiss MicroImaging GmbH, Oberkochen, Germany). Contrast and brightness were adjusted using Adobe Photoshop CS2.

Microbe inhibition assays

A fungus, Aspergillus niger, and a gram positive bacterium, Arthrobacter globiformis [91], were used to test the strength of the chemical defense. The A. niger strain was an isolate from old beetle cultures (GJ, unpublished), which was determined by the German collection of microorganisms and cell cultures (DSMZ) as Aspergillus niger, a common soil fungus also growing e.g. on bread and other food, known as ‘Black mold’. A. globiformis (from DSMZ, strain DSM20124) was another basic soil microbe and believed to have no contacts with Tribolium in nature. The culture media are listed in Text S1. Microbe lawns were made in the petri dishes using the method of [44], [92], by the modification that we put dissected abdominal glands to the holes on the lawn poked by a sterile glass pipet (one pair of glands per hole and breaking of the reservoirs in the holes) instead of freezing beetles on the plates. The plates were incubated at 25°C for 72 h (A. niger) or 28°C for 48 h (A. globiformis) respectively, and then the inhibition zones were photographed with a digital camera. The areas of the inhibition zones were measured with freeware ImageJ 1.44p.

Phenol oxidase activity assays

After RNAi, A10 beetles were harvested and frozen individually at −80°C in 150 µl Bis-Tris buffer (0.1 M, pH 7.5, sterile filtered. Bis-Tris: Fluka, SIGMA-ALDRICH Chemie GmbH, Munich, Germany, Cat. No. 14880) for at least 24 hours. To the frozen samples, a sterile steel ball (Ø 3 mm) was added each, and samples were homogenized using a GenoGrinder tissue homogenizer for 30 seconds at a speed of 1000 strokes per minute. After grinding, samples were placed on ice immediately before centrifuging three times at 6200 rpm 4°C (Eppendorf centrifuge 5810R) for five minutes to remove beetle debris. After each centrifugation step the supernatant was transferred to a new tube on ice before being centrifuged again. For measuring actual PO activity, a flat bottom 96well plate was prepared on ice with 50 µl sterile deionized water and 50 µl Bis-Tris buffer. In each well 20 µl of an individual sample extract was pipetted, or 20 µl Bis-Tris buffer when the well was serving as a blank. As PO activates the transfer of DOPA to Dopamine in insects [93], we added 50 µl L-DOPA (3,4-Dihydroxy-L-phenylalanine, SIGMA-ALDRICH Chemie GmbH, Munich, Germany, Cat. No. D9628; 4 mg/ml in Bis-Tris buffer, sterile filtered) into each well on ice. As the addition of substrate starts the reaction, plates needed to go to the Eon Microplate Spectrophotometer (Biotek Instruments, Inc., Bad Friedrichshall, Germany) immediately. Plates were read at 490 nm and 37°C with readings every two minutes for 90 minutes. After correcting the self-darkening of the substrate by subtracting the blanks, PO activity was estimated as Vmax of the linear phase of the reaction on every individual sample well (also compare with previous data [94]).

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14