Identification of putative Type-I sex pheromone biosynthesis-related genes expressed in the female pheromone gland of Streltzoviella insularis.

Species-specific sex pheromones play key roles in moth sexual communication. Although the general pathway of Type-I sex pheromone biosynthesis is well established, only a handful of genes encoding enzymes involved in this pathway have been characterized. Streltzoviella insularis is a destructive wood-boring pest of many street trees in China, and the female sex pheromone of this species comprises a blend of (Z)-3-tetradecenyl acetate, (E)-3-tetradecenyl acetate, and (Z)-5-dodecenyl acetate. This organism therefore provides an excellent model for research on the diversity of genes and molecular mechanisms involved in pheromone production. Herein, we assembled the pheromone gland transcriptome of S. insularis by next-generation sequencing and identified 74 genes encoding candidate key enzymes involved in the fatty acid biosynthesis, β-oxidation, and functional group modification. In addition, tissue expression patterns further showed that an acetyl-CoA carboxylase and two desaturases were highly expressed in the pheromone glands compared with the other tissues, indicating possible roles in S. insularis sex pheromone biosynthesis. Finally, we proposed putative S. insularis biosynthetic pathways for sex pheromone components and highlighted candidate genes. Our findings lay a solid foundation for understanding the molecular mechanisms underpinning S. insularis sex pheromone biosynthesis, and provide potential targets for disrupting chemical communication that could assist the development of novel pest control methods.

Introduction

Lepidoptera sex pheromones, which are usually secreted by female moths to attract conspecific males, play a key role in sexual communication, and are used as a monitoring and trapping tool in integrated pest management programs [1-3]. In general, moth sex pheromones are composed of two or more components in a unique ratio, and are classified into four types (Type-I, Type-II, Type-III, and Type-0) according to their site of production, chemical structure, and biosynthetic features [4]. Type-I sex pheromones are alcohols and their derivatives (acetates and aldehydes) with long straight chains (C10–C18) which are used by most moths [1, 5]. Type-II sex pheromones are composed of C17–C23 hydrocarbons with two or three double bonds at the three, six, or nine positions, or their corresponding epoxide derivatives [1, 5]. Compared with Type-I and Type-II sex pheromones, Type-III sex pheromones with one or more methyl branches possess distinct biosynthetic features, and these components include C17–C23 saturated and unsaturated hydrocarbons, as well as functionalized hydrocarbons [5]. Type-0 sex pheromones, short-chain secondary alcohols or ketones similar to some general plant volatile compounds, are utilized by the oldest non-ditrysian lineages of Lepidoptera species and are thought to represent the ancestral type of sex pheromone [5-7]. Moth sex pheromones, particularly Type-I, are mainly biosynthesized in and released from the sex pheromone gland (PG) located at the inter-segmental membrane between the eighth and ninth abdominal segments [4, 8].

The general biosynthesis pathway for Type-I sex pheromones in moths is well established; they are synthesized de novo through modified fatty acid biosynthesis pathways, and several enzymatic reactions are indispensable, including desaturation, oxidation, reduction, and acetylation [1, 4, 9–12]. All carbon atoms of the fatty acid are derived from acetyl-CoA, acetyl-CoA carboxylase (ACC) converts acetyl-CoA into the fatty acid precursor malonyl-CoA [13], and fatty acid synthetase (FAS) produces palmitic acid (C16) or stearic acid (C18) using acetyl-CoA and malonyl-CoA as substrate and NADPH as reducing agent [14-15]. Double bonds are introduced into the acyl chain at specific positions by desaturases (DESs), of which seven (Δ5 [16], Δ6 [17], Δ9 [18], Δ10 [19], Δ11 [20], Δ12 [11], and Δ14 [21]) have been identified in Lepidoptera species based on signature motifs. For instance, Δ9-desaturases have been divided into two groups: one with a substrate chain length preference of C16 >C18 (NPVE motif), and the other with a substrate chain length preference of C18 >C16 (KPSE motif) [22]. Subsequently, the unsaturated fatty acid is subjected to chain-shortening by β-oxidation, generating sex pheromone precursors of specific chain length [23], and the carbonyl carbon is modified to form an oxygenated functional group, such as an aldehyde, alcohol, or acetate ester, and these modifications involve some key biosynthesis enzymes; fatty acyl-CoA reductase (FAR) converts these acyl chains into fatty alcohols that act as actual sex pheromone components in various moths [24-26], but most fatty alcohols are either oxidized into the corresponding aldehyde by dehydrogenases [27-28] or esterified to form acetate esters by acetyltransferase (ATF) [29-31], resulting in the final functional groups.

Streltzoviella insularis (Staudinger) (Lepidoptera: Cossidae) is a destructive wood-boring pest and occurs in many provinces and cities in China. It mainly damages various street trees, such as Fraxinus americana, Ginkgo biloba, Sophora spp., and Ulmus spp., causing great economic losses to urban forestry [32-34]. The female sex pheromone of S. insularis is a blend of (Z)-3-tetradecenyl acetate (Z3-14:OAc), (E)-3-tetradecenyl acetate (E3-14:OAc), and (Z)-5-dodecenyl acetate (Z5-12:OAc) [35-36], and these acetate esters are typical of Type-I sex pheromones. These different components indicate the involvement of different desaturases, β-oxidases, and reductases during sex pheromone production. Characterization of the genes encoding putative key enzymes involved in this process may not only help to elucidate the sex pheromone biosynthesis pathway in S. insularis, but may also provide potential targets for disrupting sexual communication for pest control purposes. Hence, in the present study, we first constructed a transcriptome library of S. insularis PGs and identified a series of genes that might be involved in sex pheromone biosynthesis. Tissue expression patterns and phylogenetic analysis were performed to postulate the potential functions of the identified genes. Based on the results, we propose putative biosynthetic pathways for the sex pheromone components in S. insularis.

Materials and methods

Ethics statement

S. insularis is not on the List of Endangered and Protected Animals in China. The Beijing Municipal Bureau of Landscape and Forestry issued a permit for field collection.

Sample collection

S. insularis individuals were collected from Fraxinus americana at Beijing Forestry University North Road, Haidian District, Beijing, China, in May 2017. Damaged trunks were chopped down, taken to the laboratory, and larvae inside trunks were fed on the phloem and xylem of the host under natural environmental conditions. Adult moths were sexed after emergence according to the genitalia. The pheromone gland and associated ovipositor valves, as well as parts of the terminal abdominal segments (together abbreviated as PG) were dissected from 1-day-old and 2-day-old female adults during the scotophase, which is reported to be the calling period of this moth [35, 37]. In addition, antennae and legs were also collected at the same time, immediately placed in RNAlater (Ambion, Austin, TX, USA), and stored at -80°C.

RNA extraction

Total RNA was extracted from 15 PGs (seven PGs from 1-day-old females and eight PGs from 2-day-old females) using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions, with three biological replicates. RNA purity was evaluated with a NanoDrop 2000 instrument (Thermo, Waltham, MA, USA), and RNA concentration was measured using a Qubit RNA Assay Kit and a Qubit 2.0 Fluorimeter (Life Technologies, CA, USA). RNA integrity was determined by an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA), and RNA degradation and contamination were monitored by 1% agarose gels to ensure the quality of the RNA samples for subsequent transcriptome sequencing.

cDNA library construction and Illumina sequencing

cDNA library construction and Illumina sequencing of samples were performed at Shanghai Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China). According to the TruSeq RNA Sample Preparation Guide V2 (Illumina), mRNA was purified from total RNA using Oligo (dT) magnetic beads, then fragmented by adding fragmentation buffer. Random hexamer primers were used to synthesize first-strand cDNA, followed by synthesis of the second strand using dNTPs, RNaseH, and DNA polymerase I. All remaining overhangs were converted into blunt ends via polymerase. After end-repair, poly-A tailing, and ligation of adapters, 150–200 bp cDNA fragments were purified using an AMPure XP system (Beckman Coulter, Beverly, MA, USA), and 3μl USER Enzyme (NEB, USA) was incubated with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by incubation at 95°C for 5 min, prior to PCR amplification. PCR products were purified using an AMPure XP system, and library quality was assessed on the Agilent Bioanalyzer 2100 system. Finally, S. insularis cDNA libraries were sequenced on an Illumina Hiseq 4000 platform, and paired-end reads were generated.

Sequence assembly and functional annotation

To obtain the clean reads, the raw reads were processed to remove low-quality reads and adapter sequences. Then, GC Content, Q20 and Q30 were used to assess the sequencing quality. The qualified reads assembly was carried out with the short reads assembling program-Trinity [38]. The largest alternative splicing variants in the Trinity results were called unigenes. The annotation of unigenes was performed by the National Center for Biotechnology Information (NCBI) BLASTx searches against the non-redundant (Nr) protein database, with a cut-off E-value of 10−5. Unigenes were also annotated using other protein databases including Gene Ontology (GO) [39], Clusters of Orthologous Groups of proteins (COG) [40], and Kyoto Encyclopedia of Genes and Genomes (KEGG) [41]. The longest open reading frame (ORF) for each unigene was determined by the NCBI ORF Finder tool (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). Fragments per kilobase of exon per million mapped reads (FPKM) values were calculated by RSEM (RNA-Seq by Expectation-Maximization) with default parameters represented gene expression in S. insularis PG tissue [42].

Identification of putative genes involved in sex pheromone biosynthesis

Putative unigenes involved in sex pheromone biosynthesis of S. insularis were confirmed by analysis with the BLASTx program. All candidate pheromone biosynthesis-activating neuropeptide receptor (PBANR), acetyl-CoA carboxylase (ACC), fatty acid synthase (FAS), desaturase (DES), acyl-CoA oxidase (ACO), acyl-CoA dehydrogenase (ACD), enoyl-CoA hydratase (ECH), L-3-hydroxyacyl-CoA dehydrogenase (HCD), 3-ketoacyl-CoA thiolase (KAT), fatty acyl-CoA reductase (FAR), alcohol dehydrogenase (AD), aldehyde reductase (AR) and acetyltransferase (ATF) genes were manually checked by tBLASTn in NCBI online.

Sequence and phylogenetic analyses

Amino acid sequences of candidate desaturases were aligned with those of other insect species using ClustalW by MEGA (Version 5.0) [43]. Phylogenetic tree construction was performed using the neighbor-joining method as implemented in MEGA (Version 5.0) with a p-distance model and pairwise deletion of gaps. Bootstrap support of tree branches was assessed by re-sampling amino acid positions 1000 times [44]. Phylogenetic trees were colored and arranged using FigTree (Version 1.4.2) [45].

Expression analysis by quantitative real-time PCR (RT-qPCR)

Expression patterns of putative ACC and DES genes in different tissues (antennae, legs, and PGs) were analyzed by RT-qPCR using a Bio-Rad CFX96 PCR System (Hercules, CA, USA). Total RNA was extracted from 25 antennae, 10 legs, and 15 PGs of female moths following the method described above, and transcribed into cDNA using a PrimeScript RT Reagent Kit with gDNA Eraser (No. RR047A; TaKaRa, Shiga, Japan). Gene-specific primers were designed using Primer 3 Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) and are listed in S1 Table. The S. insularis actin gene served as an internal reference gene. Each RT-qPCR mixture was composed of 12.5 μl of TB Green Premix Ex Taq II (Tli RNaseH Plus; No. RR820A; TaKaRa), 1 μl of forward primer (10 μM), 1 μl of reverse primer (10 μM), 2 μl of cDNA, and 8.5 μl of sterilized H2O. RT-qPCR cycling parameters were as follows: 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 30 s, and 65°C to 95°C in increments of 0.5°C for 5 s to generate melting curves. To check reproducibility, each reaction for each tissue was performed with three biological replicates and three technical replicates. Negative controls without template were included in each experiment. Relative expression levels were calculated according to the comparative 2-ΔΔCt method (the amplification efficiency was close to 100% for 12 genes) [46]. Leg samples were used for calibration, and actin was used for calculating and normalizing target gene expression, and correcting for sample to sample variation. Data in the form of means ± standard error (SE) from different samples were subjected to one-way nested analysis of variance, followed by Tukey’s honestly significant difference tests, implemented in SPSS Statistics 22.0 (IBM, Chicago, IL, USA).

Results and discussion

Illumina sequencing and unigene assembly

We constructed cDNA libraries utilizing mRNAs from S. insularis PG tissue samples as template with an Illumina Hiseq 4000 platform, and included three biological replicates. A total of 63,881,910, 54,395,274, and 58,219,720 raw reads were obtained from each library. After removing low-quality reads and adaptors, we finally acquired 60,708,992, 51,561,536, and 55,208,486 clean reads, respectively (Table 1). Subsequently, assembly of all clean reads together resulted in 30,307 unigenes with an N50 value of 2072 bp, an average length of 1385 bp, and a longest length of 26,771 bp. Raw reads have been deposited in the NCBI SRA database under accession number SRP179142.

Table 1

Summary of sequencing results.

	Raw data			Clean data
	Repeat 1	Repeat 2	Repeat 3	Repeat 1	Repeat 2	Repeat 3
Read number	63,881,910	54,395,274	58,219,720	60,708,992	51,561,536	55,208,486
Base number	9,646,168,410	8,213,686,374	8,791,177,720	9,000,405,945	7,640,975,358	8,189,668,630
Q20 (%)	97.11	96.94	97.03	98.27	98.18	98.24
Q30 (%)	93.08	92.73	92.98	94.92	94.71	94.92
GC (%)	46.69	47.02	44.4	46.57	46.87	44.26

Homology searching and functional annotation

Among the 30,307 unigenes, 16,304 (53.80%) were successfully matched using the BLASTx homology search (cut-off E-value of 10−5) to entries in the NCBI Nr protein database. The best matches were obtained for Danaus plexippus sequences (30.62%), followed by Bombyx mori (25.94%), Papilio xuthus (2.54%), and Acyrthosiphon pisum (1.63%), as shown in Fig 1.

Fig 1

Species distribution based on homology searches of S. insularis unigenes against the NCBI Nr protein database.

GO annotation was used to classify the unigenes into three functional groups (molecular function, cellular component, and biological process) according to the GO categories. Of 30,307 unigenes identified in S. insularis, 8053 (26.57%) were annotated. As shown in Fig 2, 20,072 unigenes were assigned to the ‘molecular function’ category, and ‘binding’ (4141 unigenes, 43.14%) and ‘catalytic activity’ (3695 unigenes, 38.49%) were the most highly represented terms in this category. A total of 12,115 unigenes were assigned to GO terms in the ‘cellular component’ category, and ‘cell part’ (2409 unigenes, 19.88%) and ‘cell’ (2409 unigenes, 19.88%) were the most abundant terms. A further 20,072 unigenes were assigned to GO terms in the ‘biological process’ category, and the main terms were ‘cellular process’ (4329 unigenes, 21.56%) and ‘single-organism process’ (3326 unigenes, 16.57%). In addition, all unigenes were searched against the COG database for functional prediction and classification, and the results showed that 3865 unigenes (12.75%) could be assigned to 25 specific categories (Fig 3); ‘signal transduction mechanisms’ (567 unigenes, 14.67%) was the largest group, and ‘cell motility’ (5 unigenes, 0.13%) was the smallest group. Furthermore, KEGG annotation was used to divide unigenes into five KEGG pathways (cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems; Fig 4). Most unigenes were assigned to the ‘processes’ branch, and ‘global and overview maps’ (1251 unigenes, 28.07%) was the most highly represented term.

Fig 2

GO classification of S. insularis unigenes.

Fig 3

COG classification of S. insularis unigenes.

Fig 4

KEGG classification of S. insularis unigenes.

Pheromone biosynthesis-activating neuropeptide receptor (PBANR)

The biosynthesis of Type-I sex pheromones in female moths has been shown to be regulated by a C-terminally amidated 33 amino acid neuropeptide termed PBAN that is released from the subesophageal ganglion in the brain and transported through the hemolymph to the PG [47-48]. The binding of PBAN to its receptor in the PG cell membrane will induce the opening of Ca2+ channels causing the influx of extracellular Ca2+, which then initiates sex pheromone production [49-50]. PBANR, a G protein-coupled receptor (GPCR), was first cloned from the PG of Helicoverpa zea [51]. PBANR has since been identified in Bombyx mori [52] and other Lepidoptera species [49, 53]. PBANRs exist as PBANR multiple isoforms (PBANR-As, -A, -B, and -C) based on alternative splicing of the C-terminus [54]. The various isoforms play different functional roles in the ligand-induced internalization [55], a phase of GPCR feedback regulation and desensitization in diverse moth species [56-57]. Herein, we identified a single PBANR in the S. insularis PG transcriptome that is 84% identical to Mamestra brassicae PBANR isoform B (ARO85772.1) and is very low in abundance (0.56 FPKM; Table 2 and S1 Text). The number of PBANR-encoding genes in the S. insularis PG was in accordance with Plutella xylostella [25], Agrotis segetum [58], and Agrotis ipsilon [59]. In addition, previous studies identified three isoforms of PBANR in Ostrinia nubilalis [60] and Mamestra brassicae [61]. However, we did not discover other isoforms of PBANR in our transcriptomic data, which may be explained by lower expression levels in S. insularis.

Table 2

Putative sex pheromone biosynthesis-related genes identified in the S. insularis pheromone gland transcriptome.

Name	Gene length (bp)	ORF length (bp)	Intact ORF	FPKM value	Best BLASTX match
					Function	ACC number	Species	Score	E-value	Identity
PBANR
SinsPBANR	1538	1224	Yes	0.56	pheromone biosynthesis activating neuropeptide receptor isoform B	ARO85772.1	Mamestra brassicae	726	0	84%
ACC
SinsACC1	723	399	No	3.60	PREDICTED: acetyl-CoA carboxylase	XP_013185423.1	Amyelois transitella	164	2E-42	63%
SinsACC2	7616	7101	Yes	28.33	PREDICTED: acetyl-CoA carboxylase isoform X3	XP_013146614.1	Papilio polytes	8781	0	90%
FAS
SinsFAS1	312	273	No	0.28	fatty acid synthase	BAM19658.1	Papilio xuthus	153	3E-42	88%
SinsFAS2	8170	7173	Yes	90.41	fatty acid synthase	AGR49310.1	Agrotis ipsilon	3623	0	81%
SinsFAS3	301	87	No	1.00	PREDICTED: fatty acid synthase	XP_013141731.1	Papilio polytes	177	2E-49	85%
SinsFAS4	459	441	No	0.42	fatty acid synthase 1	AKD01760.1	Helicoverpa assulta	209	5E-65	62%
SinsFAS5	310	135	No	0.59	fatty acid synthase-like	XP_021208123.1	Bombyx mori	164	4E-47	69%
DES
SinsDES1	244	225	Yes	1.81	PREDICTED: acyl-CoA Delta(11) desaturase-like	XP_011561954.1	Plutella xylostella	134	2E-36	77%
SinsDES2	449	267	Yes	1.15	acyl-CoA Delta(11) desaturase-like	XP_026752209.1	Galleria mellonella	168	4E-48	80%
SinsDES3	864	825	No	0.71	acyl-CoA delta-11 desaturase	AAL16642.1	Argyrotaenia velutinana	394	5E-135	63%
SinsDES4	350	228	No	1.20	stearoyl-CoA desaturase 5-like	XP_026757907.1	Galleria mellonella	171	2E-49	75%
SinsDES5	305	192	No	0.86	stearoyl-CoA desaturase 5-like	XP_021195328.1	Helicoverpa armigera	186	5E-56	80%
SinsDES6	1278	1002	Yes	65.89	desaturase	ARD71185.1	Spodoptera exigua	496	1E-172	70%
SinsDES7	1200	996	Yes	4.82	acyl-CoA Delta(11) desaturase-like	XP_028166624.1	Ostrinia furnacalis	531	0	78%
SinsDES8	2500	1032	Yes	361.30	acyl-CoA Delta(11) desaturase	XP_028982113.1	Diachasma alloeum	345	6E-108	52%
SinsDES9	2947	1143	Yes	5.42	desaturase	AAQ74260.1	Spodoptera littoralis	590	0	74%
SinsDES10	7148	1962	Yes	3.36	acyl-CoA-delta9-3a-desaturase	ABX71810.1	Dendrolimus punctatus	628	0	87%
SinsDES11	911	393	Yes	1.15	putative C-5 sterol desaturase	KPJ05936.1	Papilio machaon	395	4E-134	81%
SinsDES12	483	111	Yes	1.47	fatty acyl desaturase	AHW98356.1	Cydia pomonella	98.6	2E-21	77%
SinsDES13	1476	984	Yes	22.89	desaturase	AIM40219.1	Cydia pomonella	581	0	85%
SinsDES14	1435	1128	Yes	0.41	desaturase	AIM40222.1	Cydia pomonella	638	0	80%
SinsDES15	1563	966	Yes	86.74	sphingolipid delta(4)-desaturase DES1	XP_004930794.1	Bombyx mori	612	0	89%
SinsDES16	1484	1017	Yes	1.66	desaturase	ARD71181.1	Spodoptera exigua	515	2E-179	72%
SinsDES17	2225	1062	Yes	426.07	acyl-CoA Delta(11) desaturase-like isoform X1	XP_021183600.1	Helicoverpa armigera	624	0	82%
ACO
SinsACO1	405	363	No	1.10	probable peroxisomal acyl-coenzyme A oxidase 1	XP_026758799.1	Galleria mellonella	215	1E-63	73%
SinsACO2	2480	2013	Yes	42.90	PREDICTED: probable peroxisomal acyl-coenzyme A oxidase 1	XP_013188704.1	Amyelois transitella	1166	0	85%
SinsACO3	2792	2070	Yes	1.92	peroxisomal acyl-coenzyme A oxidase 3	XP_022819471.1	Spodoptera litura	1181	0	80%
SinsACO4	3173	2097	Yes	11.84	peroxisomal acyl-CoA oxidase 3	AID66678.1	Agrotis segetum	1165	0	77%
SinsACO5	375	189	No	0.57	PREDICTED: probable peroxisomal acyl-coenzyme A oxidase 1	XP_014367103.1	Papilio machaon	236	5E-77	89%
SinsACO6	2104	1899	No	26.85	probable peroxisomal acyl-coenzyme A oxidase 1 isoform X1	XP_022821900.1	Spodoptera litura	964	0	73%
SinsACO7	1919	1893	No	9.16	PREDICTED: probable peroxisomal acyl-coenzyme A oxidase 1	XP_013149571.1	Papilio polytes	992	0	75%
SinsACO8	279	243	No	0.00	probable peroxisomal acyl-coenzyme A oxidase 1	XP_021195539.1	Helicoverpa armigera	181	5E-52	90%
ACD
SinsACD1	1214	768	Yes	19.22	3-hydroxyacyl-CoA dehydrogenase type-2	XP_026727946.1	Trichoplusia ni	464	1E-161	89%
SinsACD2	3886	1902	Yes	214.36	very long-chain-specific acyl-CoA dehydrogenase, mitochondrial isoform X1	XP_026737732.1	Trichoplusia ni	944	0	80%
SinsACD3	1056	774	Yes	3.34	3-hydroxyacyl-CoA dehydrogenase type-2-like isoform X1	XP_026761478.1	Galleria mellonella	429	7E-149	79%
SinsACD4	1252	933	Yes	105.77	hydroxyacyl-coenzyme A dehydrogenase, mitochondrial-like	XP_022822785.1	Spodoptera litura	581	0	89%
SinsACD5	1547	1266	Yes	145.78	short/branched-chain-specific acyl-CoA dehydrogenase, mitochondrial	XP_023946257.1	Bicyclus anynana	808	0	92%
SinsACD6	2320	1830	Yes	11.66	PREDICTED: acyl-CoA dehydrogenase family member 9, mitochondrial	XP_013192619.1	Amyelois transitella	902	0	69%
SinsACD7	3306	1236	Yes	12.58	short-chain-specific acyl-CoA dehydrogenase, mitochondrial-like isoform X1	XP_028162581.1	Ostrinia furnacalis	697	0	81%
SinsACD8	2324	1275	Yes	210.50	probable medium-chain-specific acyl-CoA dehydrogenase, mitochondrial	NP_001298861.1	Papilio xuthus	712	0	84%
SinsACD9	4692	1224	Yes	17.89	short-chain-specific acyl-CoA dehydrogenase, mitochondrial	XP_026489065.1	Vanessa tameamea	709	0	89%
ECH
SinsECH1	1207	990	Yes	10.81	PREDICTED: probable enoyl-CoA hydratase	XP_013137975.1	Papilio polytes	484	2E-168	82%
SinsECH2	1321	303	Yes	2.54	enoyl-CoA hydratase domain-containing protein 3, mitochondrial isoform X2	XP_022822616.1	Spodoptera litura	393	3E-85	82%
SinsECH3	1439	894	Yes	13.68	enoyl-CoA hydratase domain-containing protein 2, mitochondrial	XP_028167557.1	Ostrinia furnacalis	445	1E-152	79%
HAD
SinsHAD1	1214	768	Yes	19.22	3-hydroxyacyl-CoA dehydrogenase type-2	XP_026727946.1	Trichoplusia ni	464	1E-161	89%
SinsHAD2	1056	774	Yes	3.34	3-hydroxyacyl-CoA dehydrogenase type-2-like isoform X1	XP_026761478.1	Galleria mellonella	429	7E-149	79%
SinsHAD3	1252	933	Yes	105.77	hydroxyacyl-CoA dehydrogenase	AID66694.1	Agrotis segetum	575	0	87%
KAT
SinsKAT1	1395	1194	Yes	10.21	3-ketoacyl-CoA thiolase, mitochondrial-like	XP_028176321.1	Ostrinia furnacalis	491	5E-169	63%
FAR
SinsFAR1	1805	1545	No	5.57	PREDICTED: fatty acyl-CoA reductase 1-like	XP_013185409.1	Amyelois transitella	761	0	71%
SinsFAR2	1867	1692	No	1.62	fatty acyl reductase 5	ATJ44463.1	Helicoverpa armigera	816	0	73%
SinsFAR3	2335	1875	Yes	33.44	fatty acyl-CoA reductase 2	ADI82775.1	Ostrinia nubilalis	992	0	80%
SinsFAR4	457	354	No	0.85	fatty acyl reductase	ARD71192.1	Spodoptera exigua	193	3E-56	77%
SinsFAR5	967	723	No	0.56	fatty acyl-CoA reductase 1	XP_021197389.1	Helicoverpa armigera	360	4E-118	57%
SinsFAR6	2395	1494	Yes	476.06	fatty acyl reductase	AID66655.1	Agrotis segetum	441	1E-143	46%
SinsFAR7	2040	1575	Yes	15.48	putative fatty acyl-CoA reductase CG5065	XP_004925992.1	Bombyx mori	900	0	84%
SinsFAR8	2910	1578	Yes	0.60	putative fatty acyl-CoA reductase CG5065	XP_026483533.1	Vanessa tameamea	1019	0	92%
SinsFAR9	1820	1605	No	9.83	fatty acyl reductase	ARD71186.1	Spodoptera exigua	726	0	73%
SinsFAR10	1807	1560	Yes	2.70	fatty acyl-CoA reductase 1	XP_021197389.1	Helicoverpa armigera	783	0	73%
SinsFAR11	1875	1500	Yes	65.72	putative fatty acyl-CoA reductase CG5065	XP_028038252.1	Bombyx mandarina	792	0	73%
SinsFAR12	2448	1590	Yes	35.68	putative fatty acyl-CoA reductase CG5065	XP_022835056.1	Spodoptera litura	635	0	64%
SinsFAR13	4732	1533	Yes	24.19	putative fatty acyl-CoA reductase CG8306	XP_004930778.1	Bombyx mori	855	0	79%
AD
SinsAD1	1221	975	Yes	28.73	alcohol dehydrogenase	BAR64763.1	Ostrinia furnacalis	529	0	80%
SinsAD2	1746	813	Yes	13.93	alcohol dehydrogenase AD1	AII21999.1	Sesamia inferens	360	4E-118	66%
SinsAD3	2923	1131	Yes	35.49	alcohol dehydrogenase class-3	XP_021189392.1	Helicoverpa armigera	658	0	94%
SinsAD4	1395	1059	Yes	7.76	alcohol dehydrogenase	BAR64764.1	Ostrinia furnacalis	579	0	80%
SinsAD5	1209	750	Yes	37.20	alcohol dehydrogenase AD2	AKQ06148.1	Cydia pomonella	327	7E-108	71%
AR
SinsAR1	853	807	No	6.56	aldo-keto reductase AKR2E4-like	XP_028160456.1	Ostrinia furnacalis	377	2E-128	69%
SinsAR2	1125	1011	Yes	27.65	aldo-keto reductase AKR2E4-like	XP_028177948.1	Ostrinia furnacalis	506	5E-177	71%
SinsAR3	1738	1092	Yes	25.10	aldo-keto reductase AKR2E4-like	XP_022830935.1	Spodoptera litura	553	0	75%
SinsAR4	1252	1032	Yes	26.15	PREDICTED: aldo-keto reductase AKR2E4-like	XP_013136681.1	Papilio polytes	498	8E-174	70%
SinsAR5	1242	987	Yes	125.24	aldehyde reductase 7	ATJ44502.1	Helicoverpa armigera	507	1E-177	71%
ATF
SinsATF1	1775	1269	Yes	23.49	acetyl-CoA acetyltransferase, mitochondrial	XP_028157143.1	Ostrinia furnacalis	777	0	89%
SinsATF2	379	285	Yes	0.47	PREDICTED: acetyl-CoA acetyltransferase, mitochondrial isoform X2	XP_013192024.1	Amyelois transitella	177	7E-51	79%

Acetyl-CoA carboxylase (ACC)

The first step of saturated long-chain fatty acid biosynthesis is the ATP-dependent carboxylation of acetyl-CoA to malonyl-CoA catalyzed by ACC, a rate-limiting enzyme [13-14]. In the S. insularis PG transcriptome, we identified two ACCs with lengths of 723 and 7616 bp (Table 2 and S1 Text), similar to the numbers reported previously for other moth species (two in A. ipsilon [59], one in P. xylostella [25], and one in A. segetum [58]). SinsACC1 with an ORF of 399 bp encodes for an ACC with 63% amino acid identity with the ACC of Amyelois transitella (XP_013185423.1), and SinsACC2 has an intact ORF of 7101 bp that shares high amino acid identity (90%) with the ACC of Papilio polytes (XP_013146614.1). The RT-qPCR results (Fig 5) showed that SinsACC1 was more strongly expressed in the antennae than in the other tissues, whereas SinsACC2 was mainly expressed in the PG. However, both were present in low abundance (3.6 and 28.33 FPKM) in the S. insularis PG transcriptome. It was reported that the plastid-specific ACC is inhibited by herbicides that target the eukaryotic form of the enzyme in monocotyledonous plants [62-64]. Eliyahu et al. (2003) subsequently demonstrated that the herbicide diclofop inhibits PBAN-activated sex pheromone production in Helicoverpa zea, thereby implicating ACC plays a key regulatory role in fatty acid biosynthesis [65], which provides a basis for the development of a new pest control method based on disruption of sex pheromone production in females.

Fig 5

Expression profiles of putative ACCs and DESs in different S. insularis tissues.

A, antennae; L, legs; P, pheromone glands. Actin was used as an internal reference gene for normalizing target gene expression. Standard errors are represented by error bars, and different lowercase letters (a–c) above bars denote significant differences (p <0.05).

Fatty acid synthase (FAS)

FAS is the multifunctional protein that catalyzes acetyl-CoA, malonyl-CoA, and NADPH through-multienzyme complex that catalyzes the synthesis of long-chain fatty acids. Labeling studies demonstrated that palmitic acid (C16) and stearic acid (C18) are the FAS products in most moth PGs [15, 66–67]. Herein, we identified five FASs with lengths ranging from 301 bp to 8170 bp in the S. insularis PG transcriptome (Table 2 and S1 Text), These results are similar to those reported for other insects, with six and three FASs in A. segetum [58] and Sesamia inferens [68], respectively. Among the five FASs, only SinsFAS2 has an intact ORF. BLASTX results showed that FASs share high sequence similarity with Lepidoptera FASs in the NCBI Nr protein database (>60%). The FPKM analysis showed that SinsFAS2 displayed the highest expression level (90.41 FPKM) in the S. insularis PG.

Desaturase (DES)

Double bonds are introduced into the fatty acid chain at specific positions by a variety of desaturases [69]. Three putative sex pheromone compounds of S. insularis were identified as Z3-14:OAc, E3-14:OAc, and Z5-12:OAc, which are unsaturated fatty acids with acetate esters as the functional group. It is therefore reasonable to assume that the saturated fatty acid precursor of S. insularis sex pheromones is palmitic acid (C16), which is desaturated by Δ5-desaturase and Δ9-desaturase to form the precursors Z/E5-16:acyl-CoA and Z9-16:acyl-CoA in the production of two major (Z3-14:OAc and E3-14:OAc) and one minor (Z5-12:OAc) sex pheromone component, respectively (Figs 6 and 7). From the S. insularis PG transcriptome, we identified 17 putative DESs with lengths ranging from 244 to 7148 bp (Table 2 and S1 Text). The number of DESs identified in S. insularis was more than that in A. ipsilon [59], P. xylostella [25], and A. segetum [58]. Of these DESs, the identity of the best BLASX match in the NCBI NR database ranged from 52% to 89%. Notably, SinsDES15 identified in the S. insularis transcriptome shared the highest identity (89%), comparable with DES1 in Bombyx mori (XP_004930794.1). Of the 17 DESs, nine DES sequences were either less than 1000 bp, or no common sites were found for computing distances; thus, we only used the remaining eight S. insularis DES sequences to construct our phylogenetic tree (Fig 8). In the tree, SinsDES13 and SinsDES16 are clustered in the ‘Δ11-desaturases’ clade. The SinsDES17 sequence shares high sequence homology with ‘Δ9-desaturases’, and it clusters with other enzymes also possessing the NPVE motif. The remaining DESs clustered into the ‘other desaturases’ ortholog clade. The qRT-PCR results (Fig 5) revealed that SinsDES6 and SinsDES8 were highly expressed in S. insularis PG compared with the other tissues, suggesting that they may play roles in S. insularis sex pheromone production. The other five DESs (SinsDES1, SinsDES3, SinsDES4, SinsDES6, and SinsDES7) were expressed at significantly higher levels in antennae than in other tissues. All DESs except SinsDES8 and SinsDES17 were present at low abundance (from 0.41 to 86.74 FPKM) in the S. insularis PG transcriptome. DESs play important roles in the generation of structural diversity in Lepidopteran sex pheromone biosynthesis, owing to the evolution of diverse enzymatic properties [22]. Based on the most likely sex pheromone biosynthetic pathways in S. insularis, both the Δ5- and Δ9-desaturase are likely involved, but it is not clear which of the 17 desaturase genes identified in our study encode these enzymes. Further biochemical analyses of these desaturases are required to determine which ones are involved in pheromone biosynthesis.

Fig 6

Putative biosynthesis pathway of the sex pheromone components Z3-14:OAc and E3-14:OAc in S. insularis.

The saturated fatty acid precursor palmitic acid (16:0) is desaturated by Δ5-desaturase to form the precursor Z/E5-16:acyl-CoA in the production of two major pheromone components (Z3-14:OAc and E3-14:OAc).

Fig 7

Putative biosynthesis pathway of the sex pheromone component Z5-12:OAc in S. insularis.

The saturated fatty acid precursor palmitic acid (16:0) is desaturated by Δ9-desaturase to form the precursor Z9-16:acyl-CoA in the production of the minor pheromone component Z5-12:OAc.

Fig 8

Neighbor-joining phylogenetic tree of selected Lepidopteran DES enzymes.

The stability of nodes was assessed by bootstrap analysis with 1000 replicates, and only bootstrap values ≥0.5 are shown at the corresponding nodes. The scale bar represents 2.0 substitutions per site. S. insularis sequences are colored red.

β-oxidation enzymes

After a specific Δ5 or Δ9 double bond is introduced into palmitic acid to form a fatty acyl CoA precursor, the chain of the precursors is then shortened sequentially via a β-oxidation catabolic process to generate different shorter chain pheromone precursors (14C and 12C). Each cycle of β-oxidation involves four reactions: (1) acyl-CoA oxidases (ACOs, in peroxisomes) and acyl-CoA dehydrogenases (ACDs, in mitochondria) act on acyl-CoA to form E2-enoyl-CoA; (2) E2-enoyl-CoA is reversibly hydrated by enoyl-CoA hydratase (ECH) to form L-3-hydroxyacyl-CoA; (3) L-3-hydroxyacyl-CoA dehydrogenase (HAD) catalyzes the reversible dehydrogenation of L-3-hydroxyacyl-CoA to 3-ketoacyl-CoA; and (4) 3-ketoacyl-CoA is cleaved by 3-ketoacyl-CoA thiolase (KAT) [37, 70–72]. In the S. insularis PG transcriptome, we identified eight ACO genes, nine ACD genes, three ECH genes, three HAD genes, and one KAT gene (Table 2 and S1 Text). The derived protein sequences of these 24 transcripts share 63–92% amino acid identity with their homologs in other insects. All transcripts were present in low abundance (from 0 to 214.36 FPKM) in the S. insularis PG.

Fatty acyl-CoA reductase (FAR)

Chain-shortened fatty acyl CoA precursors are reduced to the corresponding alcohols by alcohol-generating FARs. Fatty alcohols can serve as sex pheromone components in many moths including Plutella xylostella [25]. Herein, we detected 13 transcripts homologous to putative FAR genes in the S. insularis PG transcriptome (Table 2 and S1 Text), similar to the number identified in other moth species (13 in A. ipsilon [59] and 10 in A. segetum [58]). Among them, SinsFAR6 was expressed at the highest level (476.06 FPKM). The FARs in S. insularis encode proteins shared 46–92% amino acid sequence identity with homologs in other Lepidoptera moths such as B. mori, Helicoverpa armigera, and Spodoptera exigua.

Alcohol dehydrogenase (AD)

Fatty alcohols can also be used as pheromone intermediates to produce corresponding aldehydes by ADs [73]. In the S. insularis PG, five homologous full-length AD genes were identified (Table 2 and S1 Text). The number of AD-encoding genes in S. insularis was in accordance with P. xylostella [25] and A. ipsilon [59]. Two ADs (SinsAD1 and SinsAD4) encode proteins that are homologous to ADs in Ostrinia furnacalis (BAR64763.1 and BAR64764.1) and share relatively high amino acid sequence identity (70%); SinsAD2 encodes a protein sharing 66% identity with Sesamia inferens AD1 (AII21999.1), SinsAD3 encodes a protein sharing 94% identity with the AD of Helicoverpa armigera (XP_021189392.1), and SinsAD5 encodes a protein sharing 71% identity with the AD of Cydia pomonella (AKQ06148.1). FPKM value analysis revealed low expression levels in the S. insularis PG for all five ADs (FPKM <50).

Aldehyde reductase (AR)

ARs are a group of the aldo-keto reductases that catalyze the reduction of fatty aldehydes to alcohols [74]. Whether ARs first produce aldehydes which are then converted to alcohols, or vice versa, is very difficult to distinguish in sex pheromone biosynthesis. Herein, we identified five AR genes in the S. insularis PG transcriptome, and four included intact ORFs (Table 2 and S1 Text). The number of ARs identified in S. insularis was less than that in A. ipsilon [59] and P. xylostella [25]. The deduced protein sequences of these five genes share high amino acid sequence identity (>60%) with their homologs in other Lepidoptera species, and all were expressed at low levels (from 6.56 to 125.24 FPKM) in the S. insularis PG.

Acetyltransferase (ATF)

ATF catalyzes the conversion of fatty alcohols to acetate esters, and this is the final enzyme in the pheromone biosynthetic pathway of the S. insularis. Previous studies showed that ATF is found almost exclusively in the PG, and is active during the photophase and all adult stages [75-76]. ATF is microsomal and exhibits specificity for the Z isomer of 12-, 14-, and 16-carbon monounsaturated fatty alcohol substrates [29–30, 75–76]. However, the enzyme has not been identified at the gene level in any moth so far [58]. In the present study, we identified two transcripts predicted to encode ATFs in the S. insularis PG (Table 2 and S1 Text). The number of ATF-encoding genes in the S. insularis PG was in accordance with P. xylostella [25]. The BLASTX results revealed 89% and 79% amino acid sequence identity shared with putative ATFs of Ostrinia furnacalis and Amyelois transitella (XP_028157143.1 and XP_013192024.1), respectively. Both ATF transcripts were present at low abundance (23.49 and 0.47 FPKM) in the S. insularis PG.

Primers used for RT-qPCR analysis of ACCs and DESs in S. insularis.

(DOCX)

Click here for additional data file.

Nucleic acid sequences of all putative sex pheromone biosynthesis-related genes identified in the S. insularis pheromone gland transcriptome.

(DOCX)

Click here for additional data file.