Introduction

The rose-grain aphid, Metopolophium dirhodum (Walker) (Hemiptera: Aphididae) is one of the most common and economically important aphid pests of cereals, including wheat, barley, rye, and oat, worldwide (1–4). M. dirhodum is native in the Holarctic and then introduced into South America, South Africa, Australia, and New Zealand (5,6). Under the continental climate of central Europe, M. dirhodum is usually the most abundant aphid species on cereals (5,7,8). In China, M. dirhodum was first recorded in the 1980s and then gradually spread from the western to eastern side of the wheat growing regions, resulting in increased crop yield reduction (9). M. dirhodum damages cereals by sucking the juice from wheat leaves, stems, and young ears, further resulting in the deterioration of plant nutrition (10). This aphid defecates sticky honeydew that further obstructs photosynthesis and reduces wheat quality (11) and transmits a number of pathogenic plant virus, including the barley yellow dwarf virus (12). The nymphs and adults of this aphid may cause yield losses of 27%–30% during the latter part of flowering stages of wheat (10,13).

Wing polymorphism is commonly observed in insects of various orders, including Hemiptera, Coleoptera, Hymenoptera, Orthoptera, Diptera, Lepidoptera, Isoptera, Psocoptera, and Dermaptera (14). Similar to most aphids, M. dirhodum, can produce wing morphs when exposed to crowding, poor nutrition, and temperature or photoperiod changes (15–17). Wing dimorphism in insects is an adaptive switch to environmental changes. In particular, wingless morphs allocate additional resources to reproduction, enabling rapid colony growth. Meanwhile, winged morphs focus on dispersal, which enables them to look for new habitats and food resources. Compared with wingless ones, winged morphs are better at long distance migration and host alternation, thus causing more serious host damage and virus transmission (14,17). Great breakthrough has been recently achieved in the exploration on the molecular mechanism underlying the wing differentiation of planthopper. Transcription factors, Zfh1 and FoxO, were verified to regulate alternative wing morphs by faithfully relaying the insulin signaling activity in Nilaparvata lugens, providing a new dimension to the molecular explanations of this unique feature in insects (18–20). Furthermore, miR-34 could regulate wing dimorphism by mediating the cross-talk among insulin signaling pathway, 20-hydroxyecdysone (20E), and juvenile hormone in a positive feedback manner in N. lugens (21).

Aphids and planthoppers may develop two different strategies in wing dimorphism. Planthopper nymphs can respond to environmental cues to develop into long- or short-winged morphs (14). Meanwhile, wing polymorphism in aphids is usually transgenerational. The parthenogenetic mother could produce either winged or wingless offspring in response to different environmental factors (22). Hence, different regulation mechanisms may exist between planthoppers and aphids. Vellichirammal et al. (2017) observed that a relatively high number of wingless offspring were produced after treatment with 20E or its analogue. Conversely, many winged offspring were produced when ecdysone signaling was suppressed by RNA interference (RNAi) targeting the ecdysone receptor (EcR) (23). A hypothesized model suggested that the elevated maternal ecdysone signaling may suppress the embryonic insulin/insulin-like growth factor (IIS) signaling and result in the high expression of FoxO-targeted genes for producing wingless offspring; the opposite is true for producing winged offspring (24). These results indicated that 20E is an essential regulatory factor underlying transgenerational plasticity in wing-dimorphic aphid.

In RNA editing, the sequence of an RNA is posttranscriptionally altered (25, 26). Adenosine (A) to inosine (I) editing is the most prevalent type of RNA posttranscriptional modification catalyzed by adenosine deaminase (adenosine deaminase acting on RNA, ADAR) that uses double-stranded RNAs (dsRNAs) as substrate. Owing to the highly similar structures of I and guanosine (G), cellular machineries usually recognize I as G during translation (27, 28). High-throughput sequencing has greatly facilitated the genome-wide identification of RNA editing events, and A-to-I RNA editing sites have been systematically characterized in various organisms, such as humans (29), macaques (30), mice, (31), worms (32), and cephalopods (33). In this work, we generated a high-quality reference genome of M. dirhodum and performed transcriptome analyses to systematically identify the A-to-I RNA editing sites of all transcripts.

The A-to-I mRNA editing events of numerous genes in 20E biosynthesis and signaling were examined, and two synonymous Ato-I RNA editing sites on CYP18A1 were identified to be closely associated with transgenerational wing dimorphism induced by crowding. CYP18A1 plays key roles in 20E metabolism, and its abnormal expression could significantly change the 20E titer, resulting in developmental disorders in multiple insect species (34). Furthermore, we identified that one synonymous A-to-I mRNA editing site may change the binding ability of miR-3036-5p to CYP18A1, thus affecting the CYP18A1 expression, 20E metabolism, and wing dimorphism of offspring in M. dirhodum. These findings revealed a new crucial regulatory mechanism of phenotypic plasticity in aphids.

Results

Chromosomal-level de novo genome of M. dirhodum

A total of 41.17 Gb of high-quality paired-end reads were obtained by Illumina genomic sequencing (~92.22X coverage, Table S1). The genome size of M. dirhodum was estimated to be 457.2 M based on k-mer counting. k-mer distribution analysis revealed a peak at 79.8x of the sequencing depth that suggested a moderate level of heterozygosity (0.445%) and highly repetitive sequence content (59.20%) in the genome (Fig. 1A), similar to those in other Aphidinae insects with low or moderate level of heterozygosity (11, 35). To obtain a reference genome for M. dirhodum, we generated 161.53 Gb of PacBio long reads using a CCS model (Table S1) and subsequently corrected them into 10.34 Gb HiFi reads. The genome was initially assembled using hifiasm, resulting in 296 contigs with a contig N50 of 7.82 Mb and the longest contig of 23.64 Mb (Table 1). The 41.17 Gb of short reads generated by Illumina NovaSeq 6000 platform were then mapped against our assembly, resulting in a mapping rate of 92.18%. BUSCO analysis showed that 96.9% (singlecopied gene: 92.5%, duplicated gene: 4.4%) of the 1,367 single-copy genes in the insecta_odb10 database were identified as complete, 0.4% were fragmented, and 2.7% were missing in the assembled genome. The percentage of complete single-copy genes is higher than those in the genomes of some other insect species, such as Sitobion miscanthi (90.2%) (11), Rhopalosiphum maidis (94.5%) (35), Acyrthosiphon pisum (93.5%) (36), and Eriosoma lanigerum (96.8%) (37). Given the moderate level of heterozygosity and the high level repetitiveness of the genome, these results indicate the high-quality genome assembly of M. dirhodum.

Fig. 1:

Table 1:

For the chromosome-level assembly, 38.09 Gb of clean reads (150 bp paired-end) were obtained from the Hi-C library (coverage: 85.31X, Table S1). A total of 118,367,396 (86.83%) reads were mapped to the draft genome. Among these sequences, 96,331,684 (70.67%) were uniquely mapped and then analyzed with 3D-DNA software to assist genomic assembly. Sixty-eight scaffolds were assembled with an N50 length of 37.54 Mb (Table 1). Finally, 447.8 Mb genomic sequences (accounting for 98.50% of the whole assembled length) were located on nine chromosomes (Fig. 1B, Table 1, Table S2), which is identical to S. miscanthi (11). The contig N50 and scaffold N50 of M. dirhodum were also higher than those of previously reported aphid genome assemblies (Table 1). This is the first high-quality chromosome-level genome of M. dirhodum, which will be very helpful for the cloning, functional verification, and evolutionary analysis of genes in this important species or even other Hemiptera insects.

Genome annotation

RepeatMasker (38) and RepBase (39) were used to annotate repeat sequences. In total, 34.97% of the M. dirhodum genome was annotated as repeat sequences. Long terminal repeats (LTRs), long interspersed nuclear elements (LINEs), and DNA transposons accounted for 9.23%, 2.25%, and 10.33% of the whole genome, respectively, and 13.16% of repeat sequences were annotated as unclassified. A total of 286 tRNAs were predicted by trnascan-SE. Using infernal, we also identified 51 small nucleolar RNAs (snoRNAs), 586 ribosomal RNAs (rRNAs), 73 small nuclear RNAs (snRNAs), 59 microRNAs (miRNAs), 286 tRNAs, and 639 other types of ncRNAs.

After repeat sequences were masked, 18,003 protein-coding genes with a mean CDS length of 1,776 bp were identified from the M. dirhodum genome using de novo, homology-, and RNA sequencing-based methods. The number of genes in the M. dirhodum genome is comparable with that of several other Aphidinae species, such as S. miscanthi with 16,006 protein-coding genes (11) and Diuraphis noxia with 19,097 protein-coding genes (40), but far less than those of R. maidis (35), A. pisum (36), Myzus persicae (41), and E. lanigerum (37) with 26,286, 36,195, 23,910, and 28,186 proteincoding genes, respectively (Table 1). Functional annotation found that 16,548 (91.92%), 9,030 (50.16%), and 12,836 (71.30%) genes had significant hits with the proteins catalogued in NR, SwissProt, and eggNOG, respectively. A total of 9,260 (51.44%) and 6,254 (34.74%) genes were annotated to GO terms and KEGG pathway, respectively (Fig. S1).

Genome synteny and phylogeny analysis

Whole genome-based phylogenetic analysis was performed with eight other hemipteran insect species, namely, M. persicae (41), D. noxia (40), A. pisum (36), R. maidis (35), Melanaphis sacchari (GCA_002803265.2), N. lugens (42), Bemisia tabaci (43, 44), and Apolygus lucorum (45), to gain insights into an evolutionary perspective for M. dirhodum. Drosophila melanogaster (46) was used as the outgroup. A total of 209,881 genes were assigned to 22,945 orthogroups for the 10 species (Fig. 1D). A phylogenetic tree was constructed using the single-copy orthologous genes (Table S3). M. dirhodum and the five other Aphididae insects formed an Aphididae cluster, showing that M. dirhodum is close to A. pisum but separated from M. sacchari and R. maidis (35). Three other Hemiptera insects, namely, B. tabaci (43, 44), N. lugens (42), and A. lucorum (45), formed another cluster (Fig. 1D).

The syntenic relationships between M. dirhodum and A. pisum genome were compared. The results reveal high levels of genome rearrangement between the chromosomes of M. dirhodum and A. pisum and a number of fission and fusion events (36). Chr1 in M. dirhodum shares 81.9% of the syntenic blocks of chr X in A. pisum (Fig. 1C). Considering the conservation of X chromosome in Aphidini insects (47), we inferred that chr 1 might be the sex chromosome in M. dirhodum. In addition, chrA1 in A. pisum is mainly syntenic to chr2, chr4, chr5, and chr8 in M. dirhodum; ChrA2 in A. pisum is mainly syntenic to chr6, chr7, and chr9 in M. dirhodum; and ChrA3 in A. pisum is mainly syntenic to chr3 in M. dirhodum (36). However, many fusion events covering small regions occurred in all chromosomes between these two insect species (Fig. 1C).

Identification of A-to-I RNA-editing sites

RNA-seq was performed among winged and wingless M. dirhodum third- and fourth-instar nymphs and adults using our assembled genome as a reference to identify potential A-to-I RNA-editing sites in this aphid. After processing, a total of 11678 A-to-I RNA-editing sites were obtained in M. dirhodum (Table S4); of which, 1557 are located on chromosome-1, 1397 on chromosome-2, 1472 on chromosome-3, 1300 on chromosome-4, 1257 on chromosome-5, 898 on chromosome-6, 1300 on chromosome-7, 1197 on chromosome-8, 1299 on chromosome-9, and 1 on scaffold 30 (Fig. 2A, B). Among the A-to-I RNA-editing sites, 4356 (37.3%) are in intergenic regions, and 7323 (62.7%) in 3090 protein-coding genes, 342 (2.9%) in 5′ UTRs, 4515 (38.7%) in introns, and 466 (4.0%) in 3′ UTRs; 515 (4.4%) are nonsynonymous (in CDS regions and causes amino acid changes when edited) and 1485 (12.7%) are synonymous (in CDS regions but do not cause amino acid changes) (Fig. 2B). The average level of these editing sites on chromosome-1 to -9 ranges from 57.2% to 59.3%, showing only a marginal difference (Fig. 2C). Given the different lengths of the nine assembled chromosomes, the A-to-I RNA-editing density of chromosome-1 (the sex chromosome, chromosome-X) is lower than those of the other eight chromosomes (Fig. 2A).

Fig. 2:

A-to-I RNA-editing on CYP18A1 is linked to transgenerational wing dimorphism under crowding

20E is an essential control factor underlying transgenerational plasticity in wing-dimorphic aphids. An environmentally regulated maternal ecdysteroid hormone can mediate wing dimorphism in the next generation (23, 48). The A-to-I RNA editing of numerous genes in 20E biosynthesis and signaling pathway (Fig. 3A, Table S4) (49), including CYP306A1, CYP302A1, CYP315A1, CYP18A1, ECR, FTZ-F1, and E74, was identified. Further verification was performed to confirm whether the A-to-I RNA editing sites on these genes are involved in wing dimorphism (Fig. 3B).

Fig. 3:

Transgenerational wing dimorphism was observed in M. dirhodum in which crowding of the parent (100 mother aphids in a 10 cm³ tube) increased the winged offspring (a total of 255 offspring were used to calculate the proportion of winged and wingless individuals in the crowding group, n=255) by 36.4% (Fig. 3E) compared with that under normal conditions (10 mother aphids in a 10 cm³ tube) (n=277). Nevertheless, 20E treatment for the crowded parent significantly decreased the number of winged offspring (n=272) by 31.2% (Fig. 3E). The proportions of two synonymous Ato-I RNA editing sites (editing site 528 and 759, the 528th and 759th nucleotide in the CDS region of CYP18A1) (Fig. 3C, D) on CYP18A1 were significantly higher in the parents treated with crowding (93.3% and 36.4%) and crowding+20E (86.6% and 35.6%) than those in the normal group (21.3% and 20.0%) (RNA extraction was individually performed from 150 aphids and then used for reverse transcription and gene amplification in normal, crowding, and crowding+20E groups) (Fig. 3F). No evident difference was observed for the A-to-I RNA editing and expression levels of the other genes involved in 20E biosynthesis and signaling pathway after crowding treatment (Fig. S2, Table S5). Given that CYP18A1 is an essential enzyme in 20E metabolism, these results indicated that the A-to-I RNA editing of CYP18A1 might be important in crowding induced wing dimorphism in M. dirhodum.

The transcriptional level of CYP18A1 significantly increased by 1.75 fold in the crowded parents compared with that in the normal group, and 20E treatment can further increase its transcription (5.97 fold) (Fig. 3G). Western blot assay results showed that CYP18A1 expression in crowding (2.90 fold) and crowding+20E (2.80 fold) treatment groups was also higher than that in the normal group (Fig. 3H). Moreover, the titer of 20E was significantly decreased in the crowded parent (Fig. 3I). Thus, we inferred that crowding could affect CYP18A1 expression and alter its A-to-I RNA editing level in M. dirhodum, thereby controlling the 20E titer in mother aphid to regulate the wing dimorphism of offspring.

RNAi-mediated knockdown of ADAR2 could regulate the A-to-I RNA-editing level and expression of CYP18A1

Insects have completely lost ADAR1, and ADAR2 homologue from mammals may be extremely critical in A-to-I RNA editing (Fig. S3) (50). One ADAR2 located on chromosome-2 was identified from our genome annotation result. The transcriptional level of ADAR2 was 2.19 fold higher in the crowding+20E treatment parent than that in the normal group, but no significant difference was identified between the crowding and normal groups (Fig. 3J). These results indicated that ADAR2 expression may not be affected by crowding but can be significantly induced by 20E.

The RNAi-mediated knockdown of CYP18A1 and ADAR2 was performed in crowded parent aphids to further explore whether Ato-I RNA-editing could regulate CYP18A1 expression and to determine the role of CYP18A1 in transgenerational wing dimorphism. RT-qPCR results showed that compared with that of the control (dsEGFP), the relative CYP18A1 expression decreased by 51.4% and 42.2% (Fig. 4A) at 48 h post feeding of dsCYP18A1 and dsADAR2, respectively. Similar results were also obtained by Western blot assay (Fig. 4B). The relative ADAR2 expression decreased by 46.7% (Fig. 4C) at 48 h post feeding of dsADAR2. These findings indicated that the RNAi-mediated knockdown of ADAR2 can affect the abundance of CYP18A1. Moreover, the percentages of Ato-I RNA editing individuals for the two synonymous editing sites on CYP18A1 significantly declined by 61.3% and 22.4% (Fig. 4D) (RNA extraction was individually performed from 150 aphids and then used for reverse transcription and gene amplification in dsEGFP, dsCYP18A1, and dsADAR2 groups), respectively, after ADAR2 silencing.

Fig. 4:

The RNAi-mediated knockdown of CYP18A1 and ADAR2 can significantly increase the titer of 20E (Fig. 4E) and reduce the number of winged offspring by 29.6% and 24.4% (Fig. 4F) (273, 248, and 250 offspring were used to calculate the proportion of winged and wingless individuals among dsEGFP, dsCYP18A1, and dsADAR2 groups, respectively), respectively. All these results showed that ADAR2 can regulate the A-to-I RNA-editing level of CYP18A1 to affect its expression and subsequently control the maternal 20E titer, finally influencing the proportion of winged offspring.

miR-3036-5p targets on CYP18A1 in M. dirhodum

The RNA editing of miRNA binding sites within target transcripts could alter miRNA targeting (51), that is, a miRNA that preferentially targets the edited or unedited version of the transcript and then participates in editing-specific miRNA-mediated transcript degradation. Here, two miRNA-target prediction software programs, miRanda and RNAhybrid, were used to identify the miRNAs that potentially act on CYP18A1. The results showed that miR-3036-5p could bind to the sequence containing edited position (editing site 528) of CYP18A1 in M. dirhodum (Fig. 5A).

Fig. 5:

MiR-3036-5p was initially identified in A. pisum. Meanwhile, the precursor and mature sequences of miR-3036-5p, which was considered a conserved miRNA in aphids, can be excellent matches to the genome of multiple aphids (52, 53). MiR-3036-5p expression significantly decreased by 48.8% (Fig. 5B) after crowding, indicating the potential role of miR-3036-5p in crowding-induced transgenerational wing dimorphism.

RNA immunoprecipitation (RIP) assay using a monoclonal antibody against the Ago1 protein was initially performed to validate the binding between miR-3036-5p and CYP18A1. CYP18A1 was significantly enriched in the Ago1-immunoprecipitated RNAs treated by miR-3036-5p agomir compared with those treated by agomir-NC (Fig. 5C). Moreover, the colocalization signals of CYP18A1 and miR-3036-5p were detected by fluorescence in situ hybridization (FISH) assay. CYP18A1 and miR-3036-5p colocalized in multiple tissues of the whole aphid body, and no signal was observed in the negative control (Fig. 5E).

We performed reporter assays using luciferase constructs fused to the binding region of CYP18A1 to further confirm the interaction between CYP18A1 and miR-3036-5p in vitro. When miR-3036-5p agomir was cotransfected with the pmir-GLO-CYP18A1-wt (inserted sequence containing the unedited site) into HEK293T cells, the luciferase activity significantly declined by 52.7% (Fig. 5D) compared with that of the mutated control (pmir-GLO & GSTu1-long-3′UTR-mt, the binding sites complementary to the “seed” sequences of miR-3036-5p were complete mutated). When miR-3036-5p agomir was cotransfected with the pmir-GLO-CYP18A1-edit (inserted sequence containing the edited site) into HEK293T cells, the luciferase activity was not significantly changed relative to that of the mutated control but significantly declined compared with that of the empty vector control (cotransfected with the pmir-GLO) (Fig. 5D). The luciferase activity in miR-3036-5p agomir & pmir-GLO-CYP18A1-wt group significantly declined by 39.2% (Fig. 5D) compared with that in the miR-3036-5p agomir & pmir-GLO-CYP18A1-edit group, indicating the potential influence of A-to-I RNA editing to the binding between miR-3036-5p and CYP18A1.

miR-3036-5p regulates transgenerational wing dimorphism by targeting CYP18A1 in M. dirhodum

To determine the role of miR-3036-5p and its target on transgenerational wing dimorphism in M. dirhodum, we combined different treatments by feeding mother aphids with miR-3036-5p agomir and antagomir.

When miR-3036-5p was inhibited at 48 h post feeding of miR-3036-5p antagomir (Fig. 6B) under normal condition (Fig. 6A), the transcriptional and translational levels of CYP18A1 were significantly upregulated by 2.89 and 2.24 fold (Fig. 6C, D), respectively. Meanwhile, the 20E titer was declined by 31.0% (Fig. 6E), and the proportion of winged offspring increased (n=249 in antagomir-NC group and n=267 in miR-3036-5p-antagomir group) by 19.2% (Fig. 6F). Although CYP18A1 expression was significantly reduced (Fig. 6C, D) and the 20E titer was increased (Fig. 6E), the proportion of winged offspring did not change (Fig. 6F) when treated with miR-3036-5p agomir (Fig. 6B). These results showed that miR-3036-5p inhibition could induce winged offspring reproduction by regulating CYP18A1 expression, but its overexpression could not play an effective role possibly due to the extremely low proportion of winged offspring (2.80%) under normal growth condition.

Fig. 6:

Under crowding condition (Fig. 6G), miR-3036-5p agomir and antagomir treatments (Fig. 6H) had minimal effects on CYP18A1 expression (Fig. 6I, J) and could not bring about further changes in the 20E titer (Fig. 6K) and proportion of winged offspring (Fig. 6L). Given the high A-to-I RNA editing level on CYP18A1 after crowding treatment, the binding between CYP18A1 and miR-3036-5p might be destroyed. As a consequence, the inhibition or overexpression of miR-3036-5p could not effectively regulate the CYP18A1 expression and wing dimorphism of offspring.

Discussion

Aphids belong to the superfamily Aphidoidea, which is part of the insect order Hemiptera. More than 5000 species of aphid have been described, and approximately 100 of them are important agricultural pests (54). Aphids have gradually become important models to study symbiosis, insect–plant interactions, and developmental polyphenism (55). Owing to the remarkable advances in sequencing technologies, the genomes of over 20 aphid species have been assembled (56). In this research, a high-quality chromosome-level genome of M. dirhodum was first produced using PacBio long HiFi reads and Hi-C technology. A total of 447.8 Mb genomic sequences were located on nine chromosomes, and 18,003 protein-coding genes were annotated. These genomic resources developed for M. dirhodum are valuable for understanding its genetics, development, and evolution and provide important references for the study of other insect genomes.

A-to-I RNA editing is the most prevalent form of posttranscriptional modification in animals, plants, and other organisms (57). It acts through multiple mechanisms, including the alteration of protein-coding capacity, generation of diverse protein isoforms, influence on the miRNA binding ability on RNA targets, and alteration of mRNA recognition by RNA binding proteins (58). A-to-I RNA editing also plays an important role in normal physiological processes, such as body development (59), reproduction (60), and responses to environmental changes (61). However, research on A-to-I RNA editing in insect, especially in agricultural pest, is limited. In the present study, 11678 A-to-I RNA-editing sites were systematically identified in M. dirhodum by using the present assembly genome. This work is also the first systematic identification of RNA editing events in aphids.

Functional wing polymorphism is commonly observed in insects as a form of adaptation and a source of variation for natural selection (14). For many aphids, polymorphism is transgenerational (23), that is, the mother senses the environment, and her offspring responds with or without wings. High population density is a key environmental factor that induces winged morphs in aphid species (23, 62). Here, we also found that a large number of winged offspring of M. dirhodum would be generated under crowding conditions. CYP18A1, a key enzyme of 20E inactivation, is involved in crowding-mediated wing dimorphism. CYP18A1 is a cytochrome P450 enzyme with 26-hydroxylase activity, a prominent step for ecdysteroid catabolism in D. melanogaster. When CYP18A1 was transfected in Drosophila S2 cells, 20E was extensively converted into 20-hydroxyecdysonoic acid (34). In Bombyx mori, CYP18A1 is a 20E-inducible gene, and its ectopic overexpression in transgenic individuals could induce 20E titer reduction, growth arrestment, and larval lethality (63). Considering the essential role of 20E in the wing dimorphism of aphids, we speculated that CYP18A1 might be a key factor mediating this important process. Two synonymous A-to-I RNA editing sites on CYP18A1 could be affected by crowding and thus modify the expression of this gene. A-to-I RNA editing levels are usually affected by changing environments (61, 64). For example, one RNA editing site in K⁺ channels from polar octopuses greatly accelerates gating kinetics by destabilizing the open state so this species can adapt to the extremely cold environment (64).

The A-to-I RNA editing events within mRNA transcripts have the potential to create or destroy the binding site of miRNAs (51). For example, 54,707 A-to-I RNA editing sites were identified in Tianzhu white yak, and 202 A-to-I editing sites altered 23 target genes of 140 miRNAs as determined by miRNA–mRNA interaction analysis (65). Nakano et al. (66) showed that ADAR1 could positively regulate the expression of dihydrofolate reductase (DHFR), which plays a key role in folate metabolism by editing the miR-25-3p and miR-125a-3p binding sites in the 3′-UTR of DHFR, resulting in enhanced cellular proliferation and resistance to methotrexate. Here, we also found that one conserved miRNA named miR-3036-5p in this aphid can target on CYP18A1, and one synonymous A-to-I RNA editing site (editing site 528) on CYP18A1 could destroy this binding. Under normal condition, this site tends to maintain a low RNA editing level, resulting in the accurate regulation of CYP18A1 expression and wing dimorphism of offspring by miR-3036-5p. However, under crowding condition, this site has an extremely high editing level, so miR-3036-5p could not affect CYP18A1 expression. Therefore, the RNA editing level of CYP18A1 changes with the population density of mother aphids to affect the accurate regulation mediated by miR-3036-5p and the wing dimorphism of offspring in M. dirhodum. Epigenetic modifications to the genome such as DNA methylation and RNA editing are proposed to be a mechanism to regulate phenotypic plasticity by affecting the changes in transcription and subsequent phenotypes (67). RNA editing, a widespread epigenetic process, is hypothesized to be an adaptive strategy to generate phenotypic plasticity in cephalopods. Squid could rapidly employ RNA editing in response to changes in ocean temperature, and RNA editing variants of kinesin generated in cold seawater display enhanced motile properties in single-molecule experiments conducted in the cold (68). Our findings supported that insects could also use RNA editing as a mechanism to adapt to environmental pressure and generate phenotypic plasticity.

In conclusion, a high-quality chromosome-level genome of M. dirhodum was first assembled, and A-to-I RNA editing events were identified to be extensively available in this important agricultural pest. Furthermore, we found that two synonymous A-to-I RNA editing sites on CYP18A1 could be induced by population density, which is important in crowding-mediated wing dimorphism. One synonymous A-to-I RNA editing at site 528 could inhibit the binding of miR-3036-5p to CYP18A1, thus increasing CYP18A1 expression, decreasing 20E titer, and leading to the wing dimorphism of offspring in M. dirhodum (Fig. 7). Meanwhile, crowding can also inhibit miR-3036-5p expression and further increase CYP18A1 abundance, resulting in winged offspring (Fig. 7). This work is the first to report that A-to-I RNA editing is involved in insect developmental polyphenism, and the results shine a new light on the functions of RNA editing in phenotypic plasticity in agricultural pests.

Fig. 7:

Materials and methods

Insects

The Metopolophium dirhodum used in the present study was originally collected from Langfang in Hebei Province, China in 2018, and then reared on wheat seedlings in our laboratory at 22 ± 2 °C and 60% relative humidity with a 16 h light:8 h dark cycle for more than 4 years. All of the progeny were produced by asexual reproduction through parthenogenesis.

Sample preparation, library construction and sequencing

Isogenic colonies were started from a single parthenogenetic female of M. dirhodum, and maintained alone on wheat seedlings prior to the collection of insects for sequencing, respectively. 200 mg of fresh mixed M. dirhodum, (including first- to fourth-instar nymphs and winged and wingless adults) was collected for DNA extraction and genome sequencing. Total genomic DNA was extracted using a Blood & Cell Culture DNA Mini Kit according to the manufacturer’s protocol (Qiagen, Hilden, Germany). For shortread sequencing, a paired-end library (2×150 bp) with short insert sizes of approximately 500 bp was constructed using the VAHTSTM Universal DNA Library Prep Kit for Illumina V2 (Vazyme, Nanning, China) and then sequenced on an Illumina NovaSeq 6000 platform (San Diego, CA, USA). For long-read genomic sequencing, the PacBio SMRTbell 15 kb library was constructed using a SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, CA, USA) and then sequenced on the PacBio Sequel II SMRT Cell 8M platform for circular consensus sequencing (CCS) (Pacific Biosciences, CA, USA).

To assist chromosome-level assembly, the Hi-C technique was applied to capture genome-wide chromatin interactions. Approximately 200 mg of fresh mixed M. dirhodum (including first-to fourth-instar nymphs and winged and wingless adults) was ground in 2% formaldehyde to allow cross-linking of cellular protein, and approximately 100 μg of DNA was extracted. Subsequently, chromosome integrity and cross-linked protein residues were assessed. Chromatin digestion was performed with the restriction enzyme Mbo I. Biotinylated residues were added during repair of the sticky ends, and the resulting blunt-end fragments were ligated under dilute conditions (69, 70). The DNA was extracted and randomly sheared to fragments of 300–500 bp. The biotin-labeled fragments were isolated with magnetic beads. The next four steps, including end repair, dA tailing, adapter ligation and DNA purification, were accomplished by adding the corresponding reaction components sequentially. The library quantity was estimated by Qubit 2.0, an Agilent 2100 instrument (Agilent Technologies, Santa Clara, CA, USA), and quantitative PCR. The Hi-C library was then sequenced using the Illumina NovaSeq 6000 platform with paired-end 150 bp reads.

For PacBio full-length transcriptome sequencing, total RNA was isolated from fresh mixed M. dirhodum (including first-to fourthinstar nymphs and winged and wingless adults of equal quality) using an EASYspin Plus Cell/Tissue RNA Isolation Kit (Aidlab Biotechnologies, Beijing, China) and quantified using a NanoDrop ND-2000 spectrophotometer (NanoDrop products, Wilmington, DE, USA). Ten micrograms of total RNA were reverse transcribed into cDNA using a SMARTer PCR cDNA Synthesis Kit (Takara, Dalian, China) following the manufacturer’s protocols. The SMRT library was constructed using the SMRTbell template prep kit (Takara) following the manufacturer’s protocols. The library was sequenced on the PacBio Sequel II SMRT Cell 8M platform, and SMRTlink was used to obtain full-length consensus isoform sequences.

For Illumina transcriptome sequencing, total RNA was isolated from winged or wingless M. dirhodum of third-and fourth-instar nymphs and adults of equal quality using an EASYspin Plus Cell/Tissue RNA Isolation Kit (Aidlab Biotechnologies, Beijing, China) and then quantified using a NanoDrop ND-2000 spectrophotometer. cDNA libraries were constructed using a VAHTSTM mRNA-seq V3 Library Prep Kit (Vazyme, Nanjing, China). A total of 18 libraries were constructed with three biological replicates per sample. Sequencing was performed on an Illumina NovaSeq instrument (Illumina, San Diego, CA, USA), and 150 bp paired-end reads were generated.

Genome survey and assembly

The K-mer distribution was analyzed to estimate the genome size, heterozygosity, and repeat content using Illumina paired-end reads. The K-mer distribution was analyzed using the Jellyfish and GenomeScope tools based on a k value of 17 (71).

PacBio subreads were obtained from the raw polymerase reads after removal of short and low-quality reads and the adaptor sequences, which were then filtered and corrected using the pbccs pipeline with default parameters (https://github.com/PacificBiosciences/ccs). The resulting HiFi reads (high-fidelity reads) were subjected to hifiasm for de novo assembly (https://github.com/chhylp123/hifiasm). BWA v0.7.15 (https://sourceforge.net/projects/bio-bwa/files/) (72) and SAMtools v1.4 (https://sourceforge.net/projects/samtools/files/samtools/) (73) were used for read alignment and SAM/BAM format conversion. Genome assembly and completeness were assessed using the conserved genes in BUSCO v3.0.2 (https://busco.ezlab.org/) (74).

Chromosome assembly using Hi-C

The Hi-C sequence data were aligned against the draft genome using JUICER v1.6.2 (https://github.com/aidenlab/juicer) (75). The uniquely mapped sequences were analyzed using 3D-DNA software (https://github.com/theaidenlab/3d-dna) to assist genomic assembly (76). The algorithms “misjoin” and “scaffolding” were used to remove the misjoins and obtain scaffolds at the chromosomal level. The algorithm “seal” was employed to find the scaffolds that had been incorrectly removed by the “misjoin”. The heatmap of chromosome interactions was constructed to visualize the contact intensity among chromosomes using JUICER v1.6.2.

Genome Annotation

Tandem repeats and interspersed repeats were identified using Tandem Repeats Finder (TRF) v4.09 (http://tandem.bu.edu/trf/trf407b.linux64.download.html) (77) and RepeatModeler v2.0 (http://www.repeatmasker.RepeatModeler/) (78), respectively. RepeatMasker v4.1.0 (http://www.repeatmasker.org/RMDownload) was used to mask the predicted and known repeated sequences (38). tRNAscan-SE v1.4alpha (79) was used to predict tRNAs, and Infernal v1.1.3 (http://eddylab.org/) was used to search the Rfam database v11.0 with an E-value cutoff of 10^-5 to identify other types of noncoding RNAs (ncRNAs) (80). Protein-coding genes were predicted through the combination of homology-based, RNA sequencing-based, and ab initio predictions. For the homolog-based approach, the protein sequences of several related species, including A. pisum (36), R. maidis (35), Diuraphis noxia (40), Aphis gossypii (81), Aphis glycines (82) and Myzus persicae (83), were downloaded from NCBI and aligned against the assembled genome using Gene Model Mapper (GeMoMa) v1.6.1.jar (http://www.jstacs.de/index.php/GeMoMa) to refine the blast hits to define exact intron/exon positions. For the RNA sequencing-based method, the PacBio full-length transcriptome, which was obtained from the pooled sample of M. dirhodum, was used to predict the open reading frames (ORFs) with PASA (https://sourceforge.net/projects/pasa/files/stats/timeline) using default settings. For the ab initio method, two de novo programs, Augustus v3.2.2 (http://augustus.gobics.de/binaries/) and SNAP (http://snap.stanford.edu/snap/download.html), were employed with default parameters to predict genes in the repeat-masked genome sequences. All predicted genes from the three approaches were integrated with EVidenceModeler (EVM) (https://sourceforge.net/projects/evidencemodeler/) to generate high-confidence gene sets, and the untranslated regions and alternative splicing were predicted using PASA.

The gene set was annotated by aligning protein sequences to functional databases, including NR (nonredundant sequence database), Swiss-Prot, eggNOG (evolutionary genealogy of genes: Nonsupervised Orthologous Groups), GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes), using BLAST with a threshold e-value ≤ 1e-5.

Phylogeny and comparative genomics

Orthologous groups were identified using the OrthoFinder pipeline (https://github.com/davidemms/OrthoFinder) with default parameters for M. dirhodum and nine other species, including A. pisum (36), R. maidis (35), D. noxia (40), M. persicae (83), Melanaphis sacchari (GCA_002803265.2), N. lugens (42), Bemisia tabaci (43) and Apolygus lucorum (45). Drosophila melanogaster (46) was used as an outgroup. MAFFT (https://mafft.cbrc.jp/alignment/software/) was used to align each orthologous gene sequence with default parameters. RAxML was used to infer the maximum-likelihood tree with the best-fit substitution model and 1000 bootstrap replicates. Mummer (https://github.com/mummer4/mummer) was applied for the detailed collinearity analysis between M. dirhodum and A. pisum genomes.

Identification of A-to-I RNA-editing sites

Raw data of FASTQ format (18 libraries constructed from third-and fourth-instar nymphs and adults) were first processed through primary quality control. In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N, low-quality reads (lower than 5), and contaminants from the raw data. Paired-end clean reads were aligned to the assembled genome of M. dirhodum using TopHat with default parameters, and all the mapped reads were used for downstream analyses. The basic principle for identifying an A-to-I RNA-editing site is that the site must be homozygous for gDNA and a mismatch must occur between RNA and DNA. We required that (i) a candidate A-to-I RNA-editing site must be supported by at least three RNA reads that were mapped to overlapping but not identical positions in the reference genome and the site had an editing level > 5% (The editing level of a candidate editing site was calculated as the number of reads supporting editing divided by the total number of reads covering that site); (ii) RNA reads with map quality score < 30 for the candidate editing positions were discarded; and (iii) candidate sites with multiple editing types were discarded (27). Given that the aphid species in this study is parthenogenesis and produces haploid progenies, the removal of false positives resulting from heterozygous genomic SNPs is not necessary for subsequent analysis.

Effect of crowding on winged offspring production

Wingless adults newly emerged within 2 hours were collected and separately assigned to two groups: 10 aphids in a 10 ml tube (normal) and 100 aphids in a 10 ml tube (crowding) with wheat seedlings. After treatment for 24 h, the mother aphids were individually moved onto wheat seedlings in a 60 mm diameter Petri dish and allowed to produce nymphs for 12 h. Subsequently, the adult aphids were removed, and the wing phenotypes of offspring were observed after 6-7 days of development.

Double-stranded RNA (dsRNA) and miRNA treatment

The dsRNA of CYP18A1 and ADAR2 was prepared in vitro using the MEGAscript RNAi kit (Ambion, Foster City, CA, USA). Genespecific primers containing a T7 polymerase promoter sequence were designed on the E-RNAi website (a tool for the design and evaluation of RNAi reagents for a variety of species, http://www.dkfz.de/signaling/e-rnai3/). The dsRNA of enhanced green fluorescent protein (EGFP) was used as a control. All of the synthesized dsRNAs were dissolved in nuclease-free water and then quantified using a NanoDrop 2000 (Thermo Scientific, Wilmington, DE, USA). MiR-3036-5p mimics (agomir), inhibitors (antagomir), and the respective negative controls (NC agomir, NC antagomir) were synthesized by Shanghai Genepharm Co., Ltd. (Shanghai, China). The subsequent treatment was performed by feeding dsRNA or miRNA agomir/antagomir through wheat seedlings following the method of Shang et al. (2020) with some modifications (62). Wheat seedlings were first dried for 12 hours and then completely immersed in water solution containing dsRNA (2,000 ng/μL) or miRNA mimic/inhibitor (2 μM) for 1 min. The wheat seedlings were allowed to air dry for 30 min, and their roots were subsequently inserted into a 250 μL PCR tube containing 200 μL of dsRNA (2,000 ng/μL) or miRNA mimic/inhibitor (2 μM). The tube was finally transferred to a 100 mm diameter Petri dish. Fourth-instar nymphs were inoculated to these treated wheat seedlings until they reach the adult stage, and the emerged wingless adults were then transferred under normal or crowding condition for continuous dsRNA or miRNA treatment.

20-Hydroxyecdysone (20E) treatment and measurement

In brief, 20E (Cayman Chemicals, USA) was dissolved in 95% ethanol as the stock solution, then diluted to 12 mg/mL with ultrapure water, and used as the working solution (84). The wheat seedlings treated with 20E were used to feed the fourth-instar nymphs, and the treatment method was the same as above.

For quantification, 20E was extracted from adults under different treatment conditions and then detected via competitive EIA (Cayman Chemical, Ann Arbor, MI, USA) using anti-20E rabbit antiserum (Cayman Chemicals). Samples were prepared following the method of Koyama et al. (2021) (85). (I) The collected adults were washed twice in double-distilled water to remove any contaminants. (H) After the adults were dried briefly on a small piece of paper towel for about 1 min, a group of samples for one biological replicate was weighed on an ultramicro balance (20-30 mg is usually sufficient for accurate ecdysone quantification). (ffl) All weighed samples were placed in a 1.5 ml microcentrifuge tube, added with threefold volume of absolute methanol (30 μl methanol for 10 mg samples), and frozen immediately on dry ice. (V) The frozen samples were homogenized carefully using disposable pestles and a cordless hand-pestle motor and then centrifuged in tubes at 4 °C at maximum speed for 5 min. (V) The supernatant was carefully transferred into new 1.5 ml microcentrifuge tubes, and centrifugation was repeated 1-2 times until no precipitate is visible. (M) Methanol was evaporated completely using a centrifugal vacuum concentrator at room temperature for 1-2 h. Quantification assay was performed in accordance with the manufacturer’s instructions.

MiRNA target prediction and dual-luciferase reporter assay

Two miRNA target prediction software programs, miRanda (http://www.microrna.org/microrna/getDownloads.do) and RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/welcome.html), were used with default parameters to predict miRNAs that potentially target CYP18A1. The miRNAs identified from A. pisum were downloaded and used for target prediction analysis. For luciferase activity assay, the luciferase reporter plasmid was constructed by inserting the wildtype, edited, or mutated target sequences of CYP18A1 between the firefly luciferase open reading frame (ORF) and SV40 poly (A) into the pmirGLO vector (Promega, Leiden, Netherlands). The constructed vectors, miRNA agomir or NC (negative control) agomir, were transferred into HEK293T cells using the calcium phosphate cell transfection kit (Beyotime) following the manufacturer’s instructions. The activities of firefly and Renilla luciferases were measured at 48 h post transfection with the Dual Glo Luciferase Assay System (Promega). For each transfection, the luciferase activity was averaged from the results of eight replicates.

Quantitative Real-Time PCR (qRT-PCR)

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s guidelines. The first strand of complementary DNA was synthesized from 1 pg of total RNA using a PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara Biotechnology, Dalian, China). The reaction was performed on an ABI 7500 Real Time PCR system (Applied Biosystems, Foster City, CA, USA). The expression levels for each gene were normalized to β-actin and GAPDH and calculated using the 2^-ΔΔCt method (86).

RNA immunoprecipitation (RIP)

A Magna RIP Kit (Millipore, Germany) was used to perform RIP assay following previous studies. Adults were first fed with miR-3036-5p agomir and then subjected to RIP analysis 24 h later. Approximately 100 aphids were collected and homogenized in icecold RIP lysis buffer. The lysates were centrifuged at 13,600 g for 10 min at 4 °C, and the supernatant (100 μl) was incubated with 5 μg of RIP Ab+ Ago-1 antibody (Millipore) or normal mouse IgG (Millipore, negative control) beads for 12 h at 4 °C. The beads were then washed with RIP wash buffer 2-3 times. Finally, the transcript enrichment ratio for the purified RNAs was determined by qRT-PCR.

Fluorescence in situ hybridization (Fish)

Antisense nucleic acid detection probes (5‣-CAACGAACTAATCACGTTGGTGATGGCGAGACACAGCGAACCGGC-3‣) labeled with Cy3 and (5‣-TGCTGAGATCTAGAAACCAAT-3‣) labeled with FAM (GefanBio, China) were designed to detect CYP18A1 and miR-3036-5p. The random shuffled probe (5‣-UUGUACUACACAAAAGUACUG-3‣) was used as negative control. For FISH, the adults were fixed in 4% paraformaldehyde for 2 h and then treated with 0.2 M hydrochloric acid and proteinase K. The treated samples were incubated with the mRNA or miRNA probes at 65 °C for 48 h in the dark and washed in PBS five times at room temperature. Fluorescence signals were finally analyzed, and images were recorded using a Nikon Eclipse Ci microscope (Tokio, Japan).

Western blot

The antibody against CYP18A1 was synthesized by Beijing Protein Innovation Co., Ltd. (Beijing, China), and the antibody against β-actin (TransGen Biotech, Beijing, China) was used as an internal control. Total proteins of each sample were extracted using a Tissue Protein Extraction Kit (Cwbio, Beijing, China), and the concentration was determined using a BCA Protein Assay Kit (Cwbio) following the manufacturer’s protocol.

The extracted total proteins were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and then transferred onto polyvinylidene fluoride (PVDF) membranes (Millipore, Germany). The membrane was blocked with 5% skim milk (Biotopped, China) for 2 h and subsequently incubated with specific antibody. Immunoreactivity was imaged with the multifunctional molecular imaging system (Azure). Quantitative analysis of the Western blot results was performed using the program Image J.

Statistical analysis

Statistical analysis was performed using GraphPad Prism version 8.0. One-way analysis of variation (ANOVA) followed by Tukey’s multiple comparisons was used for multiple comparisons (P < 0.05), and the Student’s t-test test (*, P < 0.05; **, P < 0.01; ns, no significance) was used for pairwise comparison. All primers used in this study are listed in S6 Table.

Data Availability

Raw genome sequencing reads and RNA-seq reads were deposited in the National Center for Biotechnology Information using BioProject Accession no. PRJNA751716 and PRJNA751719. The whole genome shotgun sequencing projects have been deposited at NCBI GenBank under accession no. JAIOUA000000000.

Acknowledgements

This work was supported by the China Agriculture Research System (Grant Number: CARS-05-03A-16). We thank Dr. Yuange Duan (China Agricultural University) for his critical comments on RNA editing analysis, Dr. Yaoguo Qin (China Agricultural University) for her assistance on aphid rearing, Miss Jiajia Song (China Agricultural University) for her drawing, and Berry Genomics Corporation for technical support in Illumina, PacBio and Hi-C.

Author contributions

P. Liang and B. Zhu planned and coordinated the project. R. Wei and B. Zhu prepared the samples for PacBio and Hi-C. W.J. Hua and B. Zhu prepared the samples for Illumina sequencing. R. Wei, W.J. Hua and L. Li, performed research; W.L. Zhang, B. Zhu and P. Liang analyzed the data. B. Zhu wrote the manuscript. P. Liang revised the manuscript, and all authors approved the final manuscript.

Additional Declarations:

The authors declare no competing interests.

Supplementary Information

Fig. S1:

Fig. S2:

Fig. S3:

Fig. S4: