Transcriptome sequencing and alignment
To explore the transcript expression patterns of S. litura at different developmental stages, we collected 15 S. litura samples from three developmental stages, larvae, chrysalis, and adults. Transcriptome sequencing of the cDNA library was performed on the Illumina HiSeq 2000 platform. A total of 446,991,085 raw reads and 134.1 Gb data were obtained from all samples. After quality control of the transcripts, 436,543,284 (97.66%) clean reads were obtained. We then mapped the clean reads to the reference genome and the results showed that the mapping rate was 81.31% ~ 94.96% (Supplementary Table S1).
Transcript expression analysis
We found that 23,020 transcripts were expressed, which accounted for 88.96% of the known transcripts. Next, DEseq2 was used to standardize the raw expression matrix and then perform the differential transcripts analysis between groups. The principal component analysis (PCA) demonstrated that the larvae, chrysalis, and adult groups of S. litura were clearly separated. At the same time, there was also a clear separation between the female and male adults (Figure S1).
Developmental differences in S. litura
Through cross-comparison of DETs at three developmental stages, a total of 2,703 downregulated DETs and 4,270 upregulated DETs were identified in the larval group (LL) when compared to the chrysalis group (LC). There were 3,145 downregulated DETs and 1,847 upregulated DETs in the LC group when compared to the adult group (LA). In addition, 2,979 downregulated DETs and 3,629 upregulated DETs were identified in the LL group when compared to LA (Fig. 2A). The Venn diagram illustrated the number of common and unique DETs between pairwise comparison of groups (Figure S2).
In order to understand the biological role of DETs, GO and KEGG function enrichment analysis was performed. The upregulated DETs in the LL group (i.e., LL vs. LC) were mainly enriched in biological process (BP) and molecular function (MF) GO terms, such as lipid catabolic process (GO:0,016,042), carbohydrate metabolic process (GO:0,005,975), glutathione metabolic process (GO:0,006,749) and trehalose transport (GO:0,015,771) in BP terms, which all related to basic metabolism processes. Many enzyme-related MF terms were enriched by these DETs, such as glutathione transferase activity (GO:0,004,364) and oxidoreductase activity (GO:0,016,491) (Supplementary Table S2). In KEGG enriched metabolic pathways, detoxification metabolism-related KEGG pathways were also enriched such as drug metabolism—cytochrome P450 (ko00982) and drug metabolism—other enzymes (ko00983) (Fig. 2B and Supplementary Table S2). Downregulated DETs in the LL group were mainly enriched in the development BP terms, such as locomotion (GO:0,040,011) and muscle contraction (GO:0,006,936). In MP terms, many cell activity processes were enriched, such as structural constituent of muscle (GO:0,008,307), actin filament binding (GO:0,051,015) and actin-dependent ATPase activity (GO:0,030,898) (Supplementary Table S2). Moreover, many signaling KEGG pathways, such as ErbB signaling pathway (ko04012), Rap1 signaling pathway (ko04015) and Hippo signaling pathway (ko04390) were enriched (Supplementary Table S2).
In the LC group (i.e., LC vs. LA), upregulated DETs were mainly enriched in development BP terms, including muscle organ development (GO:0,007,517), tissue development (GO:0,009,888) and chitin metabolic process (GO:0,006,030). In MP terms, many DETs were associated with chitin metabolism, such as chitin binding (GO:0,008,061), structural constituent of pupal chitin-based cuticle (GO:0,008,011) and chitinase activity (GO:0,004,568) (Supplementary Table S3). Some basic metabolisms in the KEGG pathways were enriched, such as citrate cycle (TCA cycle) (ko00020) and oxidative phosphorylation (ko00190) (Supplementary Table S3). The majority of downregulated DETs were enriched in energy-related and behavior BP terms, including trehalose transport (GO:0,015,771), disaccharide transport (GO:0,015,766) and chorion-containing eggshell formation (GO:0,007,304). MP terms that were associated with enzyme activity were abundantly enriched, such as serine-type endopeptidase activity (GO:0,004,252) and serine hydrolase activity (GO:0,017,171) (Supplementary Table S3). For KEGG pathways analysis, many detoxification metabolism pathways were enriched, including metabolism of xenobiotics by cytochrome P450 (ko00980), drug metabolism—cytochrome P450 (ko00982) and drug metabolism—other enzymes (ko00983) (Supplementary Table S3).
Upregulated DETs in the LL group (i.e., LL vs. LA) were mainly enriched in developmental BP terms, such as axoneme assembly (GO:0,035,082) and cytoplasmic translation (GO:0,002,181). In MF terms, many terms that were associated with enzyme activity were enriched, including glutathione transferase activity (GO:0,004,364) and serine-type endopeptidase activity (GO:0,004,252) (Supplementary Table S4). In the KEGG analysis, many pathways related to energy metabolism and detoxification metabolism were enriched, such as citrate cycle (TCA cycle) (ko00020) and metabolism of xenobiotics by cytochrome P450 (ko00980) (Supplementary Table S4). Downregulated DETs in the LL group were mainly enriched in BP terms that were associated with the nervous system, including axon guidance (GO:0,007,411), neuron differentiation (GO:0,030,182) and nervous system development (GO:0,007,399) (Supplementary Table S4). There was only one KEGG pathway, lysosome (ko04142), which was enriched by these downregulated DETs (Supplementary Table S4).
In general, we detected many BP terms that were associated with development and basic metabolism and MF terms related to enzyme activity in these comparisons. We also found many terms were involved in chitin metabolism. Furthermore, several KEGG pathways related to detoxification were found.
Gene expression differences between S. litura and S. frugiperda
In order to analyze the gene expression differences between S. litura and S. frugiperda at the same age stage, we performed DETs analysis on the shared single-copy orthologous genes. The orthologous genes of two species were selected by the OrthoFinder, and 6,735 pairs of single-copy orthologs were identified after stringent screening. Among them, 6,728 were co-expressed by S. litura and S. frugiperda. Gene expression between the two species can be quantified and compared after the length correction of corresponding gene CDS region. After that, we compared differential expressions of these single-copy orthologous genes within developmental stages between the two species. We identified 454 downregulated DETs and 440 upregulated DETs in the LL group when compared to the larval group of S. frugiperda (hereafter FL, consisting of 5th instar larvae and 6th instar larvae). A total of 493 downregulated DETs and 589 upregulated DETs were identified in the LC group when compared to the chrysalis group of S. frugiperda (i.e., FC). In addition, 658 downregulated DETs and 611 upregulated DETs were identified in the LA group when compared to the adult group of S. frugiperda (i.e., FA). A heatmap was generated using DETs of single-copy orthologous genes between S. litura and S. frugiperda across the three developmental stages (Fig. 2C).
We further performed GO and KEGG enrichment analysis of DETs of single-copy orthologous genes. The results showed that upregulated DETs in the LL group (i.e., LL vs. FL) were only enriched in one MF term, structural constituent of cuticle (GO:0,042,302). There was no significant KEGG pathway enriched in upregulated DETs in the LL group. The GO analysis showed that downregulated DETs in the LL group were enriched in one BP term, olfactory learning (GO:0,008,355) and one CC term, synapse (GO:0,045,202), and they were associated with sensory and nervous system development. There were four KEGG pathways enriched by these DETs, Rap1 signaling pathway (ko04015), axon regeneration (ko04361), adherens junction (ko04520) and tight junction (ko04530) (Supplementary Table S5).
The GO analysis showed that upregulated DETs in the LC group (i.e., LC vs. FC) were mainly enriched in, for example, extracellular space (GO:0,005,615), extracellular region (GO:0,005,576) and serine-type endopeptidase activity (GO:0,004,252). There was only one enriched KEGG pathway, fanconi anemia pathway (ko03460). Downregulated DETs in the LC group were enriched in nervous system developmental related CC terms, such as neuromuscular junction (GO:0,031,594). The nervous development and cell activity KEGG pathways that were associated with these downregulated DETs included the axon regeneration (ko04361) and tight junction (ko04530) (Supplementary Table S6).
Lastly, the GO annotation analysis of upregulated DETs in the LA group (i.e., LA vs. FA) were mainly enriched in energy related BP terms such as mitochondrial translation (GO:0,032,543), ATP synthesis coupled proton transport (GO:0,015,986) and aerobic respiration (GO:0,009,060), as well as binding-related MF terms, including chitin binding (GO:0,008,061) and heme binding (GO:0,020,037) (Fig. 2D). The metabolism-related KEGG pathways associated with these DETs included citrate cycle (TCA cycle) (ko00020), oxidative phosphorylation (ko00190) and 2-oxocarboxylic acid metabolism (ko01210) (Supplementary Table S7). The GO and KEGG enrichment analysis showed that there was no significant pathway enriched of downregulated DETs in the LA group.
In summary, the GO and KEGG enrichment analysis showed that the DETs of single-copy orthologous genes between S. litura and S. frugiperda were involved in basic metabolism and development. At the larval stage, downregulated DETs in S. litura were enriched in a few pathways related to nervous system development and cell activity. At the chrysalis stage, there were more pathways associated with development in S. frugiperda than in S. litura. At the adult stage, the upregulated DETs of single-copy orthologous genes in S. litura were enriched with more terms and pathways related to energy metabolism than S. frugiperda.
Next, we focused on the detoxification genes. In the past studies, it was found that some enzymes in S. litura, such as cytochrome P450, carboxylesterase and glutathione-s-transferase, which would contribute to the host's adaptability of detoxification [22]. In order to further explore the differences of detoxification ability between S. litura and S. frugiperda, we compared the expression levels of their homologous genes involved in detoxification pathways related to P450, GST, and carboxylesterase. The heatmap showed the expression of homologous genes involved in detoxification pathways shared by both species at different stages, including glutathione transferase activity (Fig. 3A), metabolism of xenobiotics by cytochrome P450 (Fig. 3B), drug metabolism—cytochrome P450 (Fig. 3C) and glutathione metabolism (Fig. 3D) (Supplementary Table S13). The results showed that the gene expression levels involved in detoxification pathways were different in the three stage groups of the same species. For example, the expression of genes involved in the glutathione transferase activity pathway was significantly higher in the larval stage than other stages in S. litura (Fig. 3D). In addition, the gene expression levels involved in detoxification pathways were also different between species within stages. For example, gene expression of xenobiotic metabolism by the cytochrome P450 pathway was highest at the larval stage in S. litura, while the expression of genes involved in the same pathway was high at all three stages in S. frugiperda. These results clearly indicated that the expression patterns of detoxification related genes during development may differ between the two lepidoptera (Fig. 3D).
Comparative microbiota composition of S. litura and S. frugiperda
Metagenome sequencing analysis of the larval (fifth instar larvae and sixth instar larvae) stage of S. litura and S. frugiperda obtained 474,122,459 and 479,748,778 valid reads from six samples each of S. litura and S. frugiperda, respectively (Supplementary Table S8). A total of 27 phyla, 49 classes, 114 orders, 152 families, 750 genera and 2,047 species were identified in S. litura, while 26 phyla, 52 classes,124 orders, 167 families, 881 genera and 2,631 species were found in S. frugiperda samples. Firmicutes, Proteobacteria and Actinobacteria were the dominant phyla in all samples. The proportion of Firmicutes was higher in S. litura compared to S. frugiperda (95.88% and 93.00%), whereas the relative percentage of Proteobacteria composition was higher in S. frugiperda (1.30% and 3.00%; Fig. 4A). At the genus level, Enterococcus, Alphabaculovirus and Corynebacterium were the dominant genera in all samples. Enterococcus, Pediococcus and Weissella were the dominant genera of S. litura, while Enterococcus, Corynebacterium and Bacillus were the most abundant genera of S. frugiperda (Fig. 4B). Linear discriminant analysis effect size was performed to identify differences in gut bacteria at genus level between S. litura and S. frugiperda at the larval stage. The relative abundance of Pediococcus was significantly higher in S. litura than in S. frugiperda, however, the relative abundance of several genera, such as Corynebacterium, Clostridium and Glutamicibacter were significantly higher in S. frugiperda compared to S. litura (LDA > 3, p < 0.05; Fig. 4C). A heatmap that was constructed based on the abundance of the top 50 genera of each sample (Fig. 4D) also showed significant differences at the genus level.
Metagenomic functional analysis
To analyze the functional relationships of the composition of S. litura and S. frugiperda metagenomes, genes were predicted from the CAZy, Humann3 and CARD databases. A total of 302 annotated genes from six families, consisting of glycoside hydrolases (GHs), glycosyltransferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydrate binding modules (CBMs), were detected in the gut microbiota based on the Carbohydrate-Active Enzymes (CAZy) database (http://www.cazy.org). In S. litura, we obtained 170 GH subfamilies, 74 GT subfamilies, 48 CBM subfamilies, 17 CE subfamilies, 24 PL subfamilies, and 16 AA subfamilies. However, S. frugiperda samples identified 134 GH subfamilies, 64 GT subfamilies, 40 CBM subfamilies, 17 CE subfamilies, 14 PL subfamilies and 11 AA subfamilies. The most abundant CAZyme families detected were represented by subfamilies GH38 and GH85 in GH families. Moreover, LEfse identified GT31, GH13_25 and GT10 as the top three significantly more abundant CAZymes in S. frugiperda than in S. litura. While GH13_14, GH67 and GH94 were the top three significantly more abundant CAZymes in S. litura than S. frugiperda (LDA > 2, P < 0.05; Figure S3 and Supplementary Table S9).
Metagenomic functional profiling was performed using Humann3 and the pathway abundances were compared using LEfSe with the LDA > 2. In the unstratified pathway analysis, 54 pathways differed between the S. litura and S. frugiperda groups. Among them, 48 pathways, such as mycothiol biosynthesis, were more abundant in the S. frugiperda group. In contrast, six pathways, including the pentose phosphate pathway, were more abundant in the S. litura group (Fig. 4E). The relative abundance of all pathways is summarized in Supplementary data, Table S10.
Next, a total of 407 antibiotic resistance genes (AGRs) were detected in the gut microbiomes. We detected that S. frugiperda had 330 resistance gene types, while S. litura had less, only 117 gene types (Fig. 4F). Of all the AGRs, the relative abundances of AbaF, Trimethoprim, and MexD were most abundant in S. frugiperda, and Trimethoprim, vanRC, and efrA genes were most abundant in S. litura. We only found 289 genes unique to S. frugiperda and 76 genes were unique to S. litura (Supplementary Table S11). We then matched each resistance gene type to its corresponding antibiotic, and found that the number of cephalosporin and tmacrolide antibiotics were the most in all samples.