Alternative splicing fine-tunes prey shift of Coccinellini lady beetles to non-target insect

Abstract

Background

Coccinellini lady beetles have been applied as biological control agent of aphids, however, not all of these species are obligately aphidophagous. Thus, a comprehensive understanding of the molecular mechanisms behind predaceous specificity of Coccinellini lady beetles can provide important clues for evaluating their performance and ecological risk assessment in biological control. Post-transcriptional regulations act a key role in shaping organisms’ rapid adaptation to changing environment, yet, little is known about their role in the acclimation of Coccinellini lady beetles to non-target preys.

Results

In this study, we conducted a genome-wide investigation to alternative splicing (AS) dynamics in three Coccinellini species Propylea japonica, Coccinella septempunctata and Harmonia axyridis in response to feeding shift from natural prey bean aphids (Megoura japonica) to non-target insect citrus mealybugs (Planococcus citri). Compared to aphid-feeding, all three lady beetles were subject to substantial splicing changes when preying on mealybugs. Most of these differentially spliced genes (DSGs) were not differentially expressed, and regulated different pathways from differentially expressed genes, indicating the functionally nonredundant role of AS. The DSGs were primarily associated with energy derivation, organ formation and development, chemosensation and immune responses, which may promote tolerance of lady beetles to nutrient deprivation and pathogen challenges induced by prey shift. The lady beetles feeding on mealybugs moreover downregulated the generation of splicing products containing premature termination codons (PTCs) for the genes involved in energy derivation and stimulus responses, to fine-tune their protein expression and rationalize energy allocation.

Conclusion

These findings unraveled the functional significance of AS reprogramming in modulating acclimation of Coccinellini lady beetles to prey shift from aphids to non-target insects and provided new genetic clues for evaluating their ecological safety as biological control agents.

Introduction

The tribe Coccinellini is a relatively large lady beetle group, which consists of 90 genera and over 1,000 recognized species worldwide [1]. Taking their advantages in aphidophagy in principle, many Coccinellini species have been applied as biological control agent of aphids [2]. However, not all of these species are obligately aphidophagous. When lacking of aphids, Coccinellini lady beetles can feed on alternative preys, such as other hemipterous insects (e.g., heteroptera, psyllids and whiteflies) and larvae of beetles and moths, as complement to overcome food shortage [1, 3]. However, their broad prey spectrum may also cause negative environmental effects. For example, Harmonia axyridis was first introduced for controlling aphids, but they have been found to indiscriminately attack non-target native species in introduced areas and cause detriment of local biodiversity [4, 5]. Generally, the lady beetles feeding on alternative diets show worse growth status than those feeding on aphids, in terms of delayed development, lower body weight and decreased survival rates [6, 7]. While there are studies revealing that larvae of some Coccinellini species are able to develop into adults relying solely on non-aphid diets [8]. Thus, a comprehensive investigation to the mechanism behind the prey shift of Coccinellini lady beetles from aphids to non-target insects can help in evaluating their risk of nonspecific predation and provide useful guidance to their application in pest control.

Recent genomic and transcriptomic studies have largely broadened our understanding of the molecular mechanisms underlying the acclimation of lady beetles to alternative diets [8,9,10,11,12]. Chen et al. [8] showed that, the digestions of three commonly-used biological control agents of Coccinellini, Propylea japonica, Coccinella septempunctata and H. axyridis, were substantially inhibited by the prey shift from aphids to mealybugs, while their expression of chemosensory and detoxifying genes were reprogrammed to combat the adverse impacts of diet change. In contrast, in Cryptolaemus montrouzieri, the digestion process was enhanced upon the transition from feeding on the primary prey mealybugs to the alternative prey aphids [10]. These findings provided new insights into the genetic basis underlying the physiological responses to prey shift in lady beetles. However, all these studies mainly focused on the regulatory mechanisms at transcriptional level. Little attention has been paid to the post-transcriptional processes underlying prey shift in lady beetles, which leaves a gap between gene transcription and protein expression.

Alternative splicing (AS) is one of the important post-transcriptional regulatory mechanisms in eukaryotic organisms [13,14,15]. AS refers to the alternative process of precursor mRNA splicing that the final spliced products are generated by different combinations of exons, and sometimes introns [16, 17]. Thus, AS largely boosts structural and functional complexity of the proteome of eukaryotes without increase in the number of genes. There are many studies having emphasized the important roles of AS in shaping organisms’ phenotypic evolution and rapid adaptation to environmental stresses [18,19,20]. In insects, the ubiquitous AS events have been demonstrated to make key contributions to individual growth and development [21], sex determination [21, 22], behaviors [23], and phenotypic plasticity [24,25,26]. For instance, in Plutella xylostella, distinct AS patterns were observed among different developmental stages and genders, indicative of the importance of splicing reprogramming in regulating insect development and sex determination [21]. The study in Drosophila melanogaster also revealed that the proper splicing of fruitless (fru) gene plays an essential role in male courtship behavior [23]. In another example, approximately 30% of alternatively spliced genes of small brown planthopper (Laodelphax striatellus) were found to participate in the response to infection of rice stripe virus, highlighting the contribution of AS to the interactions between insects and their pathogens [27]. Given these functional importances, we hypothesize that AS plays a key role in regulating the acclimation of lady beetles to alternative preys.

To examine our assumptions, in this study, we performed a genome-wide investigation of the AS dynamics in response to the shift from aphidophagy (the primary prey) to coccidophagy (the alternative prey) in three Coccinellini species, P. japonica, C. septempunctata and H. axyridis. With the findings, we hope to unravel the post-transcriptional mechanisms that fine-tunes their prey switch to non-target insects, and provide new insights into the genetic basis behind the predaceous specialization of Coccinellini lady beetles, which would help in evaluating their ecological safety as biological control agents.

Materials and methods

Data preprocessing

In this study, we conducted AS analysis for fourth instar larvae (24 h after molting and then 12 h after starvation) of P. japonica, C. septempunctata and H. axyridis feeding on bean aphids (Megoura japonica) and citrus mealybugs (Planococcus citri), respectively, using the RNA-seq data (BioProject ID of the NCBI Sequence Read Archive (SRA): PRJNA549114) published by Chen et al. [8]. For each species, two biological replicates were included for the AS analysis. For each dataset, the quality of raw sequencing reads was assessed using the FastQC program (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). For each database, low-quality bases (quality score < 20) were first trimmed from the end of reads using the Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), and low-quality reads containing > 15 ambiguous (N) bases, with error rate > 0.1 or shorter than 50 bp were discarded from the dataset. Clean reads of each species were then aligned to the corresponding reference genome using HISAT2 v.2.2.1 for downstream analysis [28].

Identification of alternative splicing events

For each species, genome-wide AS profiles were investigated under each feeding scenario using the rMATS v.4.0.2 [29]. The AS events are classified into five primary modes by rMATS: exon skipping (SE), intron retention (IR), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), and mutually exclusive exons (MXE) (Fig. 1A). Under each feeding condition, the inclusion level of alternative spliced fragment was computed for each AS event, and the events with an average inclusion ratio > 0 and < 1 were considered to be alternatively spliced. The nucleotide motifs of splicing sites were moreover identified and characterized for alternatively spliced genes. Taking the A3SS events in P. japonica as an example, the sequences of 5 bp upstream and 10 bp downstream of the 5′ splicing sites and 10 bp upstream and 5 bp downstream of both constitutive/proximal and alternative/distal 3′ splice sites were extracted from P. japonica genome using the FastaFromBed program of the Bedtools [30]. Motif patterns of each splicing site were then visualized across all the alternatively spliced genes by the R package ggseqlogo [31]. For each species, the inclusion level of alternative fragments of each AS event was compared between the individuals feeding on mealybugs and aphids. The events with False Discovery Rate (FDR) < 0.05 were considered to be significantly different upon different feeding scenarios, and genes hosting differential splicing events were denoted as “differentially spliced genes (DSGs)”.

Functional characterization of genes with exon skipping events

In animals, SE events are generally of the highest abundance [32]. Thus, in the present study, we mainly focused on the functional importance of SE events. According to their changes in the inclusion ratio of the alternative exons, the DSGs were assigned into two categories: the DSGs with higher inclusion ratio of the alternative exons when feeding on mealybugs over aphids were denoted as “more exon-inclusive DSGs” for convenience, while those with lower inclusion level were denoted as “less exon-inclusive DSGs”. For each of the three species, Gene Ontology (GO) enrichment analysis was performed on more exon-inclusive DSGs by hypergeometric test using the enricher program of the ClusterProfiler package [33]. The GO terms with p-value < 0.05 were considered to be significantly overrepresented. Some representative GO terms were then clustered and visualized using the treeplot program of the enrichplot package [34].

We further took a deep look into the sequence features of the proteins produced by each spliced isoforms of the ATP-binding cassette (ABC) transporter gene, which plays a key role in toxicant transport out of cells to lower their toxicity to insects [35]. The amino acid sequence coded by the longest isoform was extracted as the representative protein sequence for each lady beetle species and orthologous genes among the three species were identified using the OrthoFinder [36] with default parameter settings. Multiple Expectation-Maximization for Motif Elicitation (MEME) software [37] was used to visualize conserved motifs and the their distribution patterns. The amino acid sequences with no gap were then used for protein structure prediction in the PHYRE2 [38].

Premature termination codons introduced by intron retention

Despite of the low frequency, IR events may substantially affect protein expression by introducing premature termination codons (PTCs) to mRNAs [39]. Here, we investigated the impacts of IR-induced PTCs on the completeness of the resulting proteins in each of the three lady beetle species. For each tested isoform, a stop codon was denoted as a PTC if it occurred earlier than the constitutive stop codon, which was expected to appear at the end of the isoform. Impact of PTCs on proteins depend on their locations in mRNAs [40], that is, a PTC occurring father away from the 3′ end of transcripts is supposed to cause a larger impact on protein structure. GO enrichment analysis was then conducted for the genes with significantly lower inclusion of the PTC-contained isoforms upon predation on mealybugs than aphids, which might be biologically relevant to the responses to prey shift in lady beetles.

Comparison between alternative splicing and differential gene expression

We further compared the DSGs to the differentially expressed genes (DEGs) identified by Chen et al. [8], to dissect the different roles of AS and gene transcription in responding to prey shift. The genes regulated specifically by differential splicing or differential expression were denoted as “DS-specific genes” and “DE-specific genes”, respectively, while those under the regulations of both differential splicing and differential expression were denoted as “DE-DS genes”. GO enrichment analysis was then performed for the genes of each category with p-value cutoff of 0.05.

Results

Genome-wide features of alternative splicing in three Coccinellini lady beetles

To investigate the role of post-transcriptional mechanisms in modulating the responses to prey shift in Coccinellini lady beetles, we performed genome-wide analysis of AS in four instar larvae of three species under different predations. When feeding on the primary prey aphids, a total of 13,408, 14,283 and 14,458 AS events of 5126, 4485 and 4324 genes were detected in P. japonica, C. sptempunctata and H. axyridis, respectively (Fig. 1B-D). Of them, A3SS was the most abundant AS type in all the three lady beetle species, followed by A5SS and SE. Comparatively, IR and MXE events occurred at relatively low frequencies, which accounted for 9.94 to 15.37% and 4.47 to 7.16% of the total AS events, respectively. It was noteworthy that the occurrence frequency of AS events varied apparently when switching to the alternative prey mealybugs. Particularly, compared to the predation on aphids, the total number of AS events (10,902 events of 4402 genes) in P. japonica was reduced by 18.69% when feeding on mealybugs, where the three major types, SE, A3SS and A5SS, decreased in number by up to 23.05% (Fig. 1B). By contrast, the numbers of AS events were 5.18% and 15.20% higher in C. septempunctata and H. axyridis, respectively, under the alternative feeding scenario than natural prey (Fig. 1C and D).

In all the three species, similar distributions of the inclusion level of the alternative fragments were observed between the two predation scenarios (Fig. S1A and B). In particular, for most SE events, the inclusion ratio of the alternative exon was larger than 0.8 (Fig. S1A and B), indicating that most genes prefer not to skip the alternative exons during splicing, with the advantages of producing more complete proteins than the exon-skipped form. In contrast, the isoforms with retained intron mainly remained at low abundance, which allows to avoid the potential adverse impacts caused by the extra sequence, for example, introducing PTC to the mRNA (Fig. S1A and B). Furthermore, the pattern of splicing site motifs was also highly consistent between different feeding scenarios. For instance, a conserved GT motif was detected at the 5′ ends of the spliced introns across all the A3SS events, either feeding on aphids or mealybugs (Fig. S2A and B). In addition, an AG motif with a polypyrimidine tract starting with C/T were recognized at both proximal and distal 3′ splicing sites (Fig. S2A and B), which is consistent with the common observations in animals [41].

Fig. 1

Alternative splicing modes and number of differential splicing events in three Coccinellini lady beetles Propylea japonica, Coccinella septempunctata, and Harmonia axyridis upon different feeding scenarios. (A) Schematic diagram showing five primary types of alternative splicing events we analyzed in this study. (B-D) Bar charts showing the numbers of alternative splicing events of each type when feeding on aphids (orange) or alternative prey mealybugs (blue) in P. japonica (B), C. septempunctata (C), and H. axyridis (D)

Differential splicing regulates responses to prey shift in Coccinellini lady beetles

In P. japonica, C. sptempunctata and H. axyridis, 863, 1087 and 1111 splicing events from 518, 613 and 612 genes, respectively, were significantly changed between the two feeding scenarios. Of them, SE events were accounted for the highest abundance (40–41%) among the five AS types in all the three species (Fig. 2A). GO enrichment analysis showed that, the more exon-inclusive DSGs in P. japonica were mainly enriched in growth process, neuro and muscle development, functions of circulation system, flight and responses to reactive oxygen species (ROS) and pH alteration (Fig. 2B, left panel). In C. sptempunctata, genes involved in muscle development and organization, morphogenesis of wing disc and salivary gland, sensory perception, and epithelial cell differentiation were differentially spliced upon prey shift (Fig. 2B, middle panel). Muscle development, organization and differentiation, neuro morphogenesis and specification, compound eye development, and response to decreased oxygen levels were significantly overrepresented in H. axyridis (Fig. 2B, right panel). These intuitively reflected the impacts of diet alteration on the basic growth, development, and physiological processes of the three lady beetles, and highlighted the important role of AS reprogramming in responding to prey shift.

Fig. 2

Differential splicing events in three Coccinellini species and functional enrichment for more exon-inclusive genes. (A) Number of differentially spliced genes (DSGs) of five alternative splicing modes in Propylea japonica (left panel), Coccinella septempunctata (middle panel), and Harmonia axyridis (right panel), respectively. (B) Tree diagram showing the functional group of representative Gene Ontology (GO) terms enriched for more exon-inclusive DSGs in three Coccinellini lady beetles when feeding on mealybugs (from left to right: P. japonica, C. septempunctata and H axyridis). Dot size corresponds to the number of genes in the GO term, and the color corresponds to the p-value

Alternative splicing shows ample specificity among different Coccinellini species

Differential splicing profiles induced by prey shift displayed substantial specificity among the three Coccinellini lady beetles. For more exon-inclusive genes, only 48 out of the 573 enriched GO terms were overlapped among the three species, while the majority of GO terms were specifically overrepresented in one species (Fig. 3). In C. sptempunctata, the 160 specific terms were related to muscle organ development, locomotory behavior, regulation of reproductive process and cellular component assembly involved in morphogenesis. Melanization defense response, head development and flight were exclusively enriched in P. japonica, and sensory perception, wing disc and salivary gland development, and instar larval or pupal morphogenesis were found to be specific to H. axyridis (Fig. 3). In a particular instance of the gene of ABC family, exclusion of the alternative exon disrupted the open reading frame and led to loss of three motifs at the C-terminus of the product protein in both H. axyridis and P. japonica (Fig. 4), which may affect outward transfer of toxicants and increase their accumulation in cells [35]. But no AS occurred in this gene in C. sptempunctata. These results reflected that different regulatory mechanisms were employed by the three lady beetle species to cope with the impacts of prey switch on their growth and development.

Fig. 3

Similarity among the differentially spliced genes in three Coccinellini species when feeding on alternative prey (mealybugs). Upset intersection diagram illustrated the number of the significantly enriched Gene Ontology (GO) terms shared between species or specific to a certain species. Representative GO terms in each category are listed around the upset diagram with the corresponding color

Fig. 4

Alternative splicing impacts on the structure of ATP-binding cassette (ABC) transporter proteins. Left panel: The unrooted tree was constructed based on amino acid sequences. PJAPO, Propylea japonica; CSEPT, Coccinella septempunctata; HAXYR, Harmonia axyridis. The significantly differential splicing event is highlighted by *. Middle panel: Motif organization of each ABC protein. Different conserved motifs are represented in different colored boxes. Right panel: Predicted 3D structure for the protein products of the two genes with alternative splicing (ABC.1 gene of P. japonica, and ABC.2 gene of H. axyridis)

Alternative splicing affects protein structures and functions via introduction of premature termination codons

Besides exon skipping, IR events may also pose substantial effects on protein production by introducing PTCs to mRNAs [42, 43]. In the three lady beetle species, more than 93.31% of the retained introns were found to give rise to PTCs. More important, 60–67% of these PTCs occurred in the first 80% length of transcripts, which led to a protein truncation with > 20% loss in length (Fig. 5A). In addition, the inclusion level of PTCs across genes were largely different among the three species. More PTCs in P. japonica (~ 25%) were located within the first 10% of the length of the alternative transcripts than the other two species (~ 7% and 12% in C. sptempunctata and H. axyridis, respectively) (Fig. 5A), indicative of a larger impacts of IR events on P. japonica. To alleviate the load caused by PTCs, the inclusion ratio of the retained introns in some genes were significantly reduced along with the transition from feeding on mealybugs to aphids. In P. japonica, the less intron-inclusive DSGs were mainly enriched in the biological processes related to hormone regulation, ion homeostasis and transmembrane transport, and organ development and formation, for example, wing, heart and ventral furrow (Fig. 5B). In C. sptempunctata, the GO terms of intracellular signal transduction, external stimulus responding, energy derivation and carbohydrate metabolism were overrepresented (Fig. 5C), while those associated with adult locomotory behavior, carbohydrate homeostasis and oxidative stress resistance were substantially enriched in H. axyridis (Fig. 5D). The downregulation of PTC-contained isoforms in these genes can reduce the production of aberrant proteins and guarantee the proper activities of these physiological processes upon diet shift.

Fig. 5

Premature termination codons (PTCs) introduced by intron retention (IR) in three Coccinellini species. (A) Distribution of PTCs along the transcripts in Propylea japonica (red), Coccinella septempunctata (green), and Harmonia axyridis (blue), respectively. In each panel, each bar represents a location interval of 10% of the length of transcripts, and the farther to the right, the closer to the end of transcripts. (B-D) Featured Gene Ontology (GO) terms enriched for less intron-inclusive genes in P. japonica (B), C. septempunctata (C), and Harmonia axyridis (D) when feeding on alternative diet (mealybugs)

Alternative splicing is functionally nonredundant for acclimation to prey shift

We then compared the AS-mediated responses to alternative prey to those under the regulation of transcriptional changes. In all the three species, only 33–46 DSGs were also differentially expressed upon diet shift, which accounted for 2.63–7.44% of the DEGs and 5.69–8.13% of the DSGs (Fig. 6A; Figs. S3A and S4A). Most genes were only under the regulation of gene transcription or mRNA splicing. Functional enrichment analysis showed that, in P. japonica, the DE-specific genes were enriched in behavioral response to starvation and innate immune response to bacteria and fungi, while the DS-specific genes were overrepresented for developmental growth and muscle cell development and differentiation (Fig. 6). Comparatively, in C. sptempunctata, detection of abiotic stimulus, hormone biosynthesis, eye pigmentation, and turning behavior were mainly under the regulation of differential expression, while responses to nutrient levels and starvation, immune system process, and development of salivary gland, wing disc, compound eye and muscles were mainly regulated by AS (Fig. S3B and C). In H. axyridis, the genes involved in energy derivation processes (carbohydrate, ketone and fatty-acyl-CoA metabolisms), reproductive behavior, and xenobiotic metabolism were differentially expressed, whereas those associated with sensory organ development, locomotory behavior, carbohydrate metabolism, response to toxic substance, larval midgut histolysis, digestion, and response to amino acid starvation were specifically subject to different levels of alternative splicing (Fig. S4B and C). Together, these results indicated a nonredundant role of AS in modulating the acclimation of Coccinellini lady beetles to alternative prey.

Fig. 6

Comparison between differential splicing and differential gene expression in Propylea japonica when feeding on alternative prey (mealybugs). (A) Venn diagram showing the overlap between upregulated differentially expressed genes (DEGs) and more exon-inclusive differentially spliced genes (DSGs). Circle sizes correspond to the total number of DEGs and DSGs. (B) Representative Gene Ontology (GO) terms enriched for the genes that were only differentially expressed (DE-specific genes, left panel), both differentially expressed and differentially spliced (DE-DS genes, middle panel), and only differentially spliced (DS-specific genes, right panel)

Discussion

For predaceous insects, a broad feeding spectrum is a double-edged sword for their application in biological control: a high predation plasticity would be profitable to control multiple types of pests via a single introduction, but may also allow the agents to attack non-target organisms and threaten the biodiversity of local ecosystem [44]. A comprehensive understanding of the relationships between predator and prey allows us to evaluate the predaceous efficiency and specificity of control agents, and prevent the potential side effect [45]. In this study, we presented a genome-wide characterization of AS dynamics in three Coccinellini lady beetles in response to the shift to non-target insects. All the species were subject to high levels of differential splicing upon the transition to alternative prey (Figs. 1B and 2A). Moreover, the majority of the genes of differential splicing were not differentially expressed, and these DSGs were involved in different biological pathways from DEGs (Fig. 6; Figs. S3 and S4), signifying the complementary role of splicing plasticity in regulating the acclimation of Coccinellini lady beetles to diet shift. Such a low degree of overlap between the two regulatory processes was also observed in other animals, such as butterfly [46], salmonid fish [47] and wild house mouse [48]. Compared to gene transcription, AS evolves independently and is likely to be under weaker selective constraint, thus making an important contribution in adaptive phenotypic diversification [47]. One recent study showed that, compared to non-predaceous bugs (Apolygus lucorum and Riptortus pedestris), the diet-related genes specific to assassin bug Rhynocoris fuscipes possessed more exons and were enriched for AS, but non-diet genes displayed no significant difference among the three species, suggesting that AS-mediated transcriptome diversification plays an important role in the formation of predaceous-related phenotype of R. fuscipes [49]. Thus, our study not only provides new insights into the genetic basis behind the prey switch of Coccinellini lady beetles from aphids to other non-target insects, but also opens a new avenue for investigating evolutionary history and origin of aphid-eating coccinellids. When treated with mealybugs, all the three species displayed substantial AS changes related to organ formation and development (e.g., compound eyes, midgut, salivary glands, wings, muscles and circulation system), which intuitively reflected the influences of diet shift on their growth and development. It is mainly attributed to the different nutrient compositions between natural and alternative preys, and simply relying on mealybugs is not able to provide enough nutrients required by these aphidophagous lady beetles [50]. To combat the negative impacts, the genes involved in energy derivation and starvation responses were differentially spliced in all the three Coccinellini lady beetles upon the diet shift (Figs. 2B and 5; Figs. S3B and S4B). Of them, a kinase-encoding gene Ck1alpha (casein kinase I isoform alpha) was found to be a conservative regulator of glucose metabolism in fruit fly and mammal, and the loss of function of Ck1alpha lowered the glucose tolerance of mice and disrupted their homeostasis of blood sugar [51]. The genes of sensory perception, chemosensation and locomotory behavior also showed an increased inclusion level of alternative exons when feeding on mealybugs (Fig. S4B), and these changes may facilitate the food locating and hunting of lady beetles, and help them overcome the nutrient depletion caused by undesired diet change [52]. Additionally, different insects have distinct microorganism compositions on their surface and body, thus the diet shift from aphids to mealybugs may cause pathogen challenges to these lady beetles [53]. Accordingly, immune-related genes in all the three Coccinellini species were found to be differentially spliced accompanying with the prey transition (Fig. 3). Together, these findings supported our hypothesis that genome-wide AS changes play a key role in modulating physiological and biochemical reprogramming that promotes the tolerance of Coccinellini lady beetles to prey shift.

In all the three Coccinellini species, AS, especially IR events, were found to introduce PTCs to product mRNAs (Figs. 4 and 5A), which finally affected protein structures and functions. For example, in P. japonica and H. axyridis, a PTC was introduced to the ABC protein-coding gene by a skipped exon, which caused a deletion of motifs at the C-terminus of the product protein (Fig. 4). Notably, for most genes, the transcripts with included alternative intron were maintained at low abundance (Fig. S1). The genes associated with energy derivation (e.g., organic compound oxidation, carbohydrate homeostasis and lipid metabolism), organ development, and stimulus responses displayed a significant downregulation of IR events upon the feeding transition from aphids to mealybugs (Fig. 5). These changes might be attributed mainly to the surveillance of the nonsense-mediated decay (NMD) pathway, which can identify and degrade the aberrant PTC-containing transcripts [40, 54, 55]. Thus, a coupling of AS and NMD can reduce the energy waste on translation of non-functional proteins and allocate more resources for producing other functional proteins, especially under the scenario of limited nutrient availability [39].

There are two noticeable limitations of the current work. First, we investigated AS in Coccinellini lady beetles using short reads from RNA-seq, however, it might lead to a biased estimation of IR events. The reads produced by the next-generation sequencing (e.g., 150 bp pair-end as we used in this study) are not long enough to span the entire length of many introns, and the low coverage of intronic regions may lower the accuracy of quantification of IR events [56]. Moreover, reads from non-coding RNAs, such as small nucleolar RNAs, microRNAs or circular RNAs, may also give rise to signals that mimic IR events and thus pose difficulties for distinguishing genuine IR events from false positive calls [57]. One solution to mitigate the spurious detection of IR events is using long reads generated by the third-generation sequencing platform, such as PacBio and Oxford Nanopore, which can improve the ability to capture the persistent introns in transcripts [57]. Second, this study investigated AS changes primarily based on bioinformatic analyses using published RNA-seq data. A future experimental validation would help us confirm the functional impacts of splicing isoforms and open new avenues for understanding the role of AS in diet adaptation of Coccinellini lady beetles.

In summary, our study provides a new insight into the post-transcriptional blueprint for the acclimation of aphid-feeding coccinellids to non-target insects. We showed that AS responded dynamically to the negative influence caused by the feeding of undesired preys. The little overlap between differential splicing and differential expression suggested the irreplaceable role of AS in regulating the protein repertoire of lady beetles. The coupling of AS and NMD further eliminated the mRNAs consisting of PTCs to reduce production of truncated proteins and energy waste during diet shift. These novel findings broaden our understanding of the molecular mechanisms associated with the predaceous specificity of Coccinellini lady beetles, and also provide theoretical guidance for evaluating their ecological safety as biological control agents.