Genome-wide identification and analysis of the thiolase family in insects

Data resources

In this study, 20 representative species were selected from Lepidoptera, Hymenoptera, Hemiptera, Diptera, Coleoptera, Phthiraptera, and Isoptera (Table 1). The annotated genes and genomes of B. mori were retrieved from SilkDB v3.0 (https://silkdb.bioinfotoolkits.net). The sequence information of Danaus plexippus and Heliconius melpomene were downloaded from http://metazoa.ensembl.org/. Manduca sexta was downloaded from https://i5k.nal.usda.gov/. The other sequences were retrieved from GenBank (https://www.ncbi.nlm.nih.gov/), including Papilio xuthus, Plutella xylostella, Culex quinquefasciatus, Anopheles gambiae, Drosophila melanogaster, Anoplophora glabripennis, Nicrophorus vespilloides, Tribolium castaneum, Halyomorpha halys, Acyrthosiphon pisum, Cimex lectularius, Diachasma alloeum, Apis mellifera, Bombus impatiens, Pediculus humanus, and Zootermopsis nevadensis.

Identification of insect thiolase genes

The known thiolase sequences of Homo sapiens and Mycobacterium tuberculosis were retrieved from GenBank (Table S1; Anbazhagan et al., 2014) and used as queries to perform BLASTP search (E- value <0.01) against the protein database of predicted genes in each species. Hidden Markov Model (HMM) files of Thiolase_N (PF00108) and Thiolase_C (PF02803) domains were downloaded from Pfam database (http://pfam.xfam.org/), which were used to screen the protein database of each species with hmmsearch in HMMER 3.0 (E- value <0.01). Based on BLASTP and hmmsearch analyses, the candidate thiolase genes were identified and subsequently checked by conserved domain search (CD-Search) in NCBI and hmmscan against Pfam database (E- value <1e−5). The candidate sequences that have Thiolase_N and/or Thiolase_C domains were recognized as thiolases. Those identified thiolases were used as new queries to perform BLASTP search against the protein database of each species until no more novel loci can be found. All the validated thiolase genes were used for further analysis.

Phylogenetic analysis

The protein sequences of thiolases from 20 insects, H. sapiens and M. tuberculosis were aligned using MUSCLE (Edgar, 2004). Positions that had a high percentage of gaps (>70%) were trimmed. The handled alignment of protein sequences was used for checking the most suitable model of evolution by ProtTest 3.2 (Darriba et al., 2011). Maximum-likelihood (ML) trees were reconstructed using RAxML version 8.2.12 (Stamatakis, 2014) with the most suitable model (PROTGAMMAVTF) and 500 bootstrap replicates. FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/) was used for plotting the final phylogenetic tree. The clustering and classification of the thiolase sequences in the ML tree were done using known functional properties of H. sapiens and M. tuberculosis (Anbazhagan et al., 2014).

Chromosome distribution, gene structure, and syntenic analysis

To localize the thiolase genes on chromosomes, B. mori, H. melpomene, D. melanogaster, T. castaneum, and A. mellifera were selected because their genome sequences have been assembled into chromosomes. Based on the GFF (General Feature Format) file of each species, every thiolase gene was mapped to the corresponding chromosomes. Using protein sequences of the thiolases, the precise exon/intron structures were generated through BLAT search against the genome sequences with Scipio server (https://www.webscipio.org/). The synteny events between two species were detected by Multiple Collinearity Scan toolkit (MCScanX) with the default parameters (Wang et al., 2012). The syntenic map of B. mori and H. melpomene was constructed with family_circle_plotter.java in MCScanX software.

Molecular modeling of protein structure

The three-dimensional (3D) structure prediction of insect thiolases was conducted using the homology modeling method. Structures of T1-, T2-, CT-, AB-, TFE-, and SCP2-thiolases were predicted on-line at the SWISS-MODEL Interactive Workspace (Arnold et al., 2006). The known protein that has the highest sequence similarity to the thiolase to be analyzed is used for homology modeling. The predicted models of monomer and multimer were visualized in Swiss-PdbViewer 4.1.0 (Guex & Peitsch, 1997). To understand the 3D structural similarities among the insect thiolases, all the other structures were compared with BmorT2 using magic fit algorithm in Swiss-PdbViewer, respectively. The root-mean-square distance (RMSD) values were calculated to express the structural similarity. The lower value of RMSD means higher similarity between two structures (Carugo & Pongor, 2001).

Reverse transcription-polymerase chain reaction (RT-PCR)

The various tissues on day 3 of fifth-instar larvae were dissected in the silkworm. The sex pheromone glands (PGs) from five individuals were used as one sample at each developmental stage. All the samples were preserved in RNAlater (Ambion, 98 Austin, USA) and stored at −80 °C for RNA isolation. Total RNA was extracted using Trizol reagent (Invitrogen, USA). The first strand of cDNA was synthesized by M-MLV reverse transcriptase following the manufacturer’s instructions (Promega, USA). RT-PCR primers were listed in Table S2. The silkworm RpL3 gene was used as an internal control for relative quantitative analysis of RT-PCR. PCRs were performed with the following cycling parameters: 95 °C for 3 min (min), followed by 25 cycles of 30 s (s) at 95  °C, 30 s annealing (temperatures listed in Table S2) and 30 s extension (72 °C), and a final extension at 72 °C for 10 min. The amplification products were monitored on 1.5% agarose gels.

Genome-wide identification and phylogeny of insect thiolase proteins

To identify thiolases in insects, human and M. tuberculosis thiolase protein sequences were used as queries to perform homologous searches in whole genomes. In total, 137 thiolase genes were identified in 20 insects from seven orders, and the gene numbers of the species were ranged from 4 to 15 (Table 1; Table S1). The thiolase protein sequences were used to reconstruct the maximum-likelihood phylogenetic tree (Fig. 1). Based on the nomenclature rules in humans and M. tuberculosis (Anbazhagan et al., 2014), each thiolase was named in insects. It was indicated that insect thiolases were grouped 7 classes, namely CT, T1, T2, TFE, SCP2 (type-1), TFEL (type-2), and AB (Table 1; Fig. 1). Relatively, the classification of insect thiolases was more similar to that of the human than M. tuberculosis. Unlike humans, out of 20 insect species, only P. xuthus has gene members in AB class. Interestingly, 5 out of 6 Lepidopteran insects have no CT-thiolase. Furthermore, TFEL (type-2) class was only detected in C. lectularius (Hemipter) and bacterium M. tuberculosis.

To understand the evolutionary mode, gene gain and loss of thiolases were analyzed. It was indicated that most of the gain and loss events were occurred in a certain species (Fig. 2). Especially, P. xylostella (Lepidoptera), A. pisum (Hemipter), and P. xuthus (Lepidoptera) showed more duplications after or during the formation of the species, resulting in a total number of 15, 12, and 10 genes, respectively. Except for the gene duplication of a single species lineage, the common ancestor of D. alloeum, A. mellifera and B. impatiens showed 2 duplications (Fig. 2). This phenomenon was also noted in the clade of Hemipteran H. halys, C. lectularius, and A. pisum. For those recent duplication genes, they were often phylogenetically closely related to its ancestral genes (Fig. 1). Conversely, A. mellifera (Hymenoptera), H. halys (Hemipter), and Z. nevadensis (Isoptera) presented 3, 3, and 2 gene losses during speciation, resulting in fewer genes in these species. Generally, gene gain and loss rates are important for understanding the role of natural selection and adaptation in shaping gene family sizes. For the species with more gene expansion, whether these duplicated genes play roles in adapting to special habitats deserves further study.

Gene structures of insect thiolase genes

A comparative analysis of exon-intron structures was conducted for the 137 insect thiolase genes (Figs. 3A and 3B; Fig. S2A). The insect thiolase genes have a different number of introns ranging from 0 to 22. It was indicated that only 12 genes have no intron, and 17 genes have only one intron, accounting for 21.17% in total. The intronless genes were distributed in T2 (5), SCP2 (type-1) (2), CT (2), AB (2), and TFEL (1). Previous studies revealed that introns can delay regulatory response and are selected against in genes whose transcripts need to be adjusted quickly to meet environmental challenges (Jeffares, Penkett & Bahler, 2008). The intronless thiolase genes and the genes contained fewer introns might play important roles in survival for environmental changes. In addition, the intron number and exon/intron structures of thiolase genes are very different, even the orthologous genes of different species in the same class have a large differentiation. It was suggested that the differentiation of the intron number may result in the diversification of thiolase gene structures in insects.

Chromosome distribution and gene synteny

In order to explore the chromosomal distribution of thiolase genes, five representative species were analyzed. It was indicated that most of the thiolase genes were randomly distributed on different chromosomes (Figs. S3A–S3E), for example, six thiolase genes were scattered on four chromosomes in D. melanogaster (Fig. S3A), which is similar to thiolase genes in the human (Anbazhagan et al., 2014). However, different members of a certain class are often tandem distribution, such as DmelSCP2 (type-1)-1 and DmelSCP2 (type-1)-2 in D. melanogaster (Fig. S3A) and BmorT1-1, BmorT1-2, and BmorT1-3 in B. mori (Fig. S3C). Meanwhile, we also detected the distribution of the 10 T1-thiolase genes in P. xylostella and 7 CT-thiolase genes in A. pisum, which were distributed on several small unassembled scaffolds. Whether they are distributed in tandem, we still need to wait for the scaffold sequences to be integrated into the corresponding chromosome in future. In general, tandem duplication might be the main mechanism for enlarging thiolase family in insects.

The syntenic relationships of thiolase genes were investigated among B. mori, H. melpomene, D. melanogaster, T. castaneum, and A. mellifera because their genome sequences have been assembled into chromosome levels. The results indicated that only four genes exhibited the syntenic relationships between B. mori and H. melpomene, that is, BmorT2 and HmelT2, BmorTFE and HmelTFE (Fig. S3F). Interestingly, except for the tandemly duplicated genes, amount of thiolase genes often present orthologous relationships among insect, human, and M. tuberculosis (Fig. 1). It was suggested that thiolase family is an ancient gene family. Even in insects, the age of thiolase gene differentiation is relatively long. Thus, the discovery of fewer syntenic genes implies that thiolases might mainly locate in some non-conserved genomic blocks (The Heliconius Genome Consortium, 2012).

Subcellular localization of thiolase proteins

Subcellular localization refers to the specific location of a certain protein or the expression product of a certain gene in the cell. Protein subcellular localization is closely related to protein functions (Pereto, Lopez-Garcia & Moreira, 2005; Wang et al., 2014). Only when the protein is positioned correctly can it perform normal biological functions. In this study, subcellular localization of all the 137 insect thiolase proteins was predicted by PSORT II server (https://www.genscript.com/tools/psort), which were cytosolic, mitochondrial, or peroxisomal enzymes (Fig. 3C; Fig. S2B). Generally, most of the TFE- and T2-thiolase proteins were located in the mitochondrion, T1- and CT-thiolases were cytosolic, and SCP2-thiolases were peroxisomal proteins. Previous studies suggested that the mitochondrial and peroxisomal thiolase proteins were mainly involved in the fatty acid β-oxidation pathway (Pereto, Lopez-Garcia & Moreira, 2005), and cytosolic localization was related to the biosynthesis of acetoacetyl-CoA (Kursula et al., 2005). However, in a certain class of thiolase, there are always a few exceptions to the cellular location in some species (Fig. 3C; Fig. S2B), which suggested that its function might have diverged during evolution.

Conserved domain characteristics and catalytic residues

To identify the potential domains of insect thiolase proteins (Table S1), it was performed hmmscan analysis in Pfam database. The results indicated that all the thiolases contained Thiolase_N and Thiolase_C domains (Fig. 4A). In addition to the thiolase domains, SCP2-thioloase (type-1) has a typical sterol carrier protein 2 (SCP2) domain at C terminal. Unexpectedly, some of the members in TFE, CT, and T2 classes contained a ketoacyl-synt (beta-ketoacyl synthase) domain within Thiolase_N (Table S1). We carefully checked the alignments of hmmscan search. It was found that the E value was around the threshold 1e−5, and only about 50 amino acids can be aligned, which are much shorter than 250 amino acids of the ketoacyl-synt domain (Pfam ID, PF00109). Thus, thiolases may not contain the real ketoacyl-synt domain, and just show certain similarities with it (Huang et al., 1998). Therefore, based on the domain characteristics, all the insect thiolase encoding genes were classified as 2 groups (Fig. 4A).

The conserved sequence blocks of the 20 insects, humans, and M. tuberculosis were analyzed (Fig. S4 ; Fig. 4B). The CxS-motif is the most important sequence fingerprint in the N-terminal domain, which provides the nucleophilic cysteine (Zeng & Li, 2004; Mazet et al., 2011). Except for some incomplete sequences, almost all the thiolases contained the cysteine residue (Fig. S4 ; Fig. 4B). The histidine of the GHP-motif contributes to the oxyanion hole of the thioester oxygen (Merilainen et al., 2009). It was indicated that GHP-motif was highly conserved in insects, humans, and M. tuberculosis (Fig. S4 ; Fig. 4B). The cysteine of CxGGGxG-motif provides the catalytic residue of the active sites. Except for SCP2-thiolases, the catalytic cysteine was retained in almost all the other thiolases (Fig. S4). In addition, the asparagine side chain of the NEAF-motif interacts with important catalytic water (Mazet et al., 2011). However, NEAF-motif was replaced by HDCF-motif in all of the SCP2-thiolases (Fig. S4 ; Fig. 4B). Based on the comparison of the sequence fingerprints, it was indicated that the catalytic mechanisms of the insect thiolases might be similar to that of thiolases from mammals and bacteria.

Molecular modeling of insect thiolases

In recent years, the crystal structures of some thiolases have been gradually resolved in bacteria, fish, and mammals (Harijan et al., 2013; Kim et al., 2015; Xia et al., 2019). The high sequence similarities (>60%) may help to build more accurate 3D structures for the insect thiolases (Arnold et al., 2006). Based on homology modeling using SWISS-MODEL Interactive Workspace, we found that thiolase sequences within a class were very conserved among different organisms. For instance, BmorTFE-thiolase and BmorSCP2-thiolase (type-1) shared 67.58% and 61.63% identities with its corresponding modeling templates from human (PDB ID: 6dv2.1.A) and zebrafish (6hrv.2.A), respectively. In this study, the modeling structures of some representative thiolases were presented (Figs. 5A–5F). The structural similarities of the monomeric forms were detected with magic fit in Swiss-PdbViewer (Guex & Peitsch, 1997). The RMSD values were ranged from 0.25 Å to 0.91 Å among BmorT2, BmorT1-1, DmelCT, and PxutAB-1, while BmorTFE and BmorSCP2 (type-1) shared from 1.12 Å to 2.01 Å with the others (Table S3). It was indicated that T1, T2, CT, and AB classes share more similar 3D structures than TFE-thiolase and SCP2-thiolase (type-1) (Figs. 5A–5F). This phenomenon is widespread in both humans and M. tuberculosis (Harijan et al., 2013; Anbazhagan et al., 2014). For the quaternary structure, different thiolases also have certain differences. For example, BmorT1-1 and BmorSCP2 (type-1) were homo-tetramer and homo-dimer, respectively (Figs. 5G and 5F). The results of 3D structural modeling showed that different classes of thiolase genes still present some extent divergence in tertiary or quaternary structures.

SCP2-thiolase (type-1) was widely distributed in insects, mammals, and bacteria (Table 1). One single structural gene referred to as the sterol carrier protein x (SCPx) gene encodes a full-length protein comprised of 3-oxoacyl-CoA thiolase (known as SCP2-thiolase) and sterol carrier protein 2 (Seedorf et al., 1994; Gallegos et al., 2001). The C-terminal SCP2-domain containing the peroxisomal targeting signal is needed for the targeting of full-length SCPx into the peroxisomes. The SCP2-thiolase and SCP2 protein are produced from SCPx via proteolytic cleavage by peroxisomal proteases (Seedorf et al., 1994). Based on the homology modeling, the tertiary and quaternary structures of mature SCP2-thiolase (type-1) protein were presented in Figs. 5F and 5H, respectively. For insect SCP2-thiolases, the canonical CxGGGxG-motif is also absent, and the NEAF-motif has been replaced by HDCF-motif (Fig. 4B). The previous studies indicated that HDCF-motif might provide the catalytic cysteine in bacteria, mammals, and fish (Harijan et al., 2013; Kiema et al., 2019). Based on the structural modeling, the cysteine of HDCF-motif is very close to the other two catalytic sites in protein spatial conformation (Fig. 5F). Therefore, the catalytic cysteine of the insect SCP2-thiolases might be not provided by CxGGGxG-motif but HDCF-motif.

Expression profile and potential functional diversity

To understand the potential functional diversity of the insect thiolases, the silkworm, B. mori, was used as a model organism to perform expression profile analysis in the various tissues and sex pheromone glands (PGs) at different developmental stages. In the silkworm, genome-wide microarray with 22,987 oligonucleotides was designed and surveyed the gene expression profiles in multiple tissues on day 3 of the fifth-instar larvae (Xia et al., 2007). 5 out of 6 thiolase genes were found its corresponding probes (Fig. 6A). The microarray data indicated that BmorSCP2 (type-1), BmorT2, BmorTFE, and BmorT1-1 have expression signals at least one of the 9 tissues. Relatively, BmorT2 and BmorT1-1 showed ubiquitous expressions. Meanwhile, the expression profiles of the four genes were similar between females and males, respectively.

To validate the expression profiles of the silkworm thiolase genes, the mixed male and female tissues were used to perform RT-PCR validation on day 3 of the fifth-instar larvae (Fig. 6B). In total, 5 out of the 6 thiolase genes presented expression evidence. Relatively, BmorT1-2 and BmorTFE showed predominant expressions in hemocyte and head, respectively (Fig. 6B), while BmorT1-1 was widely expressed in various tissues. In addition, sex pheromone glands of different developmental stages were used to detect the expressions of thiolase genes (Fig. 6C). Four thiolase genes presented expression signals in the silkworm PGs. Relatively, BmorT1-1 showed the highest expression on day 8 of pupae. BmorTFE, BmorT2, and BmorSCP2 (type-1) presented expressions at all the developmental stages. Interestingly, the expression levels of all three genes were declined in the mated female PGs (Fig. 6C). These expression analyses might help us understand the functional divergence of the thiolase genes in the silkworm.