Transcriptome of the Maize Leafhopper (Dalbulus maidis) and Its Transcriptional Response to Maize Rayado Fino Virus (MRFV), Which It Transmits in a Persistent, Propagative Manner.

The corn leafhopper (Dalbulus maidis) is an important vector of maize rayado fino virus (MRFV), a positive-strand RNA (+ssRNA) marafivirus which it transmits in a persistent propagative manner. The interaction of D. maidis with MRFV, including infection of the insect and subsequent transmission to new plants, is not well understood at the molecular level. To examine the leafhopper-virus interaction, a D. maidis transcriptome was assembled and differences in transcript abundance between virus-exposed and naive D. maidis were examined at two time points (4 h and 7 days) post exposure to MRFV. The D. maidis transcriptome contained 56,116 transcripts generated from 1,727,369,026 100-nt paired-end reads from whole adult insects. The transcriptome of D. maidis shared highest identity and most orthologs with the leafhopper Graminella nigrifrons (65% of transcripts had matches with E values of <10-5) versus planthoppers Sogatella furcifera (with 23% of transcript matches below the E value cutoff) and Peregrinus maidis (with 21% transcript matches below the E value cutoff), as expected based on taxonomy. D. maidis expressed genes in the Toll, Imd, and Jak/Stat insect immune signaling pathways, RNA interference (RNAi) pathway genes, prophenoloxidase-activating system pathways, and immune recognition protein-encoding genes such as peptidoglycan recognition proteins (PGRPs), antimicrobial peptides, and other effectors. Statistical analysis (performed by R package DESeq2) identified 72 transcripts at 4 h and 67 at 7 days that were significantly responsive to MRFV exposure. Genes expected to be favorable for virus propagation, such as protein synthesis-related genes and genes encoding superoxide dismutase, were significantly upregulated after MRFV exposure. IMPORTANCE The transcriptome of the corn leafhopper, D. maidis, revealed conserved biochemical pathways for immunity and discovered transcripts responsive to MRFV-infected plants at two time points, providing a basis for functional identification of genes that either limit or promote the virus-vector interaction. Compared to other hopper species and the propagative plant viruses they transmit, D. maidis shared 15 responsive transcripts with S. furcifera (to southern rice black-streaked dwarf virus [SRBSDV]), one with G. nigrifrons (to maize fine streak virus [MFSV]), and one with P. maidis (to maize mosaic virus [MMV]), but no virus-responsive transcripts identified were shared among all four hopper vector species.

INTRODUCTION

The corn leafhopper, Dalbulus maidis (Delong and Wolcott) (Hemiptera: Cicadellidae), is an important vector of maize rayado fino virus (MRFV). The virus species Maize rayado fino virus belongs to the genus Marafivirus within the family Tymoviridae (1). MRFV is transmitted by D. maidis in a persistent, propagative manner. After leafhopper acquisition of MRFV, a latent period of at least 1 week must elapse before leafhoppers become capable of transmitting the virus to new plants (2, 3). During the latent period, MRFV is thought to replicate within the insect and cross several molecular and/or physical barriers before it can be transmitted to plants, similar to other persistent, propagatively transmitted viruses (see review [4]). However, interactions of MRFV with its leafhopper vectors during the latent period, including the triggering or suppression of insect immune and antiviral responses, are not well delineated. To begin to understand these processes, we determined the D. maidis transcriptome and its transcriptional response to feeding on MRFV-infected maize.

Understanding the transcriptomic response of D. maidis to MRFV could provide key insights into how leafhoppers react to viral invasion and propagation. MRFV propagation in D. maidis has not been associated with negative phenotypes, suggesting that the virus is in what we term “benign equilibrium” with its host (see also reference 5). The equilibrium is likely to involve an interaction with the vector immune system that is sufficiently permissive to allow virus replication without disrupting essential functions required for robust vector capability.

In addition to understanding the specific interactions of D. maidis with MRFV, one of our goals was to identify immunity genes and pathways present in D. maidis. Insect immunity involves both cellular and humoral immune responses. The cellular immune response results in encapsulation of pathogens in hemocytes by phagocytosis (6, 7). The humoral immune response involves the production of antimicrobial peptides (AMPs) synthesized through several signaling pathways, such as the Toll, immune deficiency (IMD), and Janus tyrosine kinase-signal transducer and activator of transcription (Jak-Stat) pathways (8). The humoral immune response also includes melanization and encapsulation of pathogens through the prophenoloxidase-activating system (9). The genes associated with the Toll, IMD, JNK, and Jak/Stat signaling pathways have been well addressed in many insects, such as the fruit fly Drosophila melanogaster (10, 11), the red flour beetle Tribolium castaneum (12), and the mosquito Anopheles gambiae (13), and orthologs of genes in all of these pathways were identified in the transcriptome of another important leafhopper vector, the blackfaced leafhopper, Graminella nigrifrons (14). The Toll signaling pathway is activated in response to fungi and Gram-positive (G+) bacteria. In this pathway, serine proteases are triggered to cleave Spaetzle, which further activates Toll and leads to the synthesis of antimicrobial peptides (AMPs) and immune effectors (12, 15–18) (Fig. 1), while serpin is a negative regulator of the Toll signaling pathway (8). The IMD and JNK signaling pathways are triggered by Gram-negative (G−) bacteria, which will activate relish to lead the expression of other AMPs (19, 20) (Fig. 1). In response to virus infection, an insect may activate the Jak/Stat signaling pathway (reviewed in reference 21), which has been identified in D. melanogaster, Aedes aegypti, A. gambiae, T. castaneum, and other insects (reviewed in reference 22) (Fig. 1).

FIG 1

Immune response signaling pathways in Dalbulus maidis. The Toll, Imd, and Jak/Stat immune signaling pathways and prophenoloxidase-activating system pathways from other insects were compared with gene orthologs identified in the assembled D. maidis transcriptome. Gene names are based on the names used for similar genes in Drosophila (8, 18) with the number of putative genes from transcriptomic data of D. maidis given in parentheses. The yellow highlighted genes were not identified from the transcriptomic data of D. maidis.

RNA interference (RNAi) is also a component of insect response to pathogens and is an endogenous cellular process in insects triggered by double-stranded RNA (dsRNA) and leading to cleavage of mRNA complementary to small interfering RNAs, resulting in transcriptional gene silencing (see references 14, 23, and 24 and references therein). RNAi likely plays a vital role in protecting leafhoppers from virus infection (25). Many genes, proteins, and pathways associated with RNAi have been identified in D. melanogaster (23), T. castaneum (24), A. gambiae (26), and G. nigrifrons (14).

In this study, we assembled a transcriptome for D. maidis and compared it with those of other hopper species that transmit maize- and rice-infecting viruses: the leafhopper G. nigrifrons and the planthoppers Sogatella furcifera and Peregrinus maidis. From the transcriptome, we identified genes in D. maidis which are likely to be involved in insect immune responses, including components of Toll, Imd, Jak/Stat, other immune signaling pathways, the prophenoloxidase-activating system, the melanization pathway (9, 27), and RNAi-associated genes and pathways. Further, we examined the transcriptome response of D. maidis to feeding on MRFV-infected maize, which the leafhopper persistently and propagatively transmits, shortly after exposure (4 h) and in response to active virus replication at the end of the latent period (7 days) (see references 2, 3, and 28). Genes highly regulated in response to MRFV exposure were identified and compared with those identified in other hopper vector-virus systems.

RESULTS

Dalbulus maidis transcriptome.

After trimming adapter and poly(A) sequences, the D. maidis transcriptome was assembled from 1,727,369,026 ca. 100-bp paired-end transcript reads from 32 samples, each sample a pool of 25 whole adult insects. Reads were deposited in GenBank as a short sequence read archive (SRA) BioProject PRJNA579843. De novo assembly of reads yielded 45,281 (± 1,949) contigs for naive D. maidis (insects fed on healthy plants, 16 samples) with an average length of 784.4 nucleotides (nt; standard deviation [stdev] = 7.3), and N50 of 1,265.6 nt (± 22.6 stdev). For D. maidis fed on MRFV-infected plants (16 samples), 47,340 (stdev = 4,037) contigs with an average length of 758.3 nt (stdev = 13.9) and N50 of 1,184.3 nt (stdev = 41.3) were assembled. The combined D. maidis nonredundant transcriptome of 56,116 transcripts of ≥200 nt was deposited as a transcriptome shotgun assembly (TSA; SUB7857194, DDBJ/EMBL/GenBank accession GITV00000000). Benchmarking Universal Single-Copy Orthologs (BUSCO) assessment for assembly and annotation completeness showed that the combined D. maidis transcriptome contained 1,126 (82.4%) homologs of the expected 1,367 conserved insect genes (29), with homologs of another 39 (2.9%) of the expected genes present as fragments. Approximately 15% (202 transcripts) of the expected conserved gene set was missing from our transcriptome. Some of these missing genes for D. maidis may be expressed in earlier life stages, including the eggs and nymphs, which were not represented in our adult-only transcriptome. More completeness could be expected if multiple life stages and enriched tissues were represented (29, 30).

D. maidis transcript sequences were compared to the NCBI nonredundant (nr) protein and nucleotide databases using blastx- and blastn. Of the 56,116 transcripts, 29,208 D. maidis transcripts (52.0%) matched sequences in the nr databases with E values less than 10−2, 41% (23,001) of them matching insect sequences (Table 1). The D. maidis transcript sequences were also blastx searched against a custom database that included the curated RefSeq database (31) plus NCBI nucleotide nr sequences identified from insects. Of D. maidis transcripts, 27,277 (48.6%) matched sequences in this database with E values less than 10−5, 26,885 (47.6%) with top hits to insect sequences in this insect sequence-enriched database (Table 1). Of the 26,885 sequences matching insects, 24,639 (92.3%) matched sequences from insects in the order Hemiptera (Table S1). The majority of top hits were to sequences from leafhoppers in the family Cicadellidae, including Graphocephala atropunctata, Cuerna arida, and Homalodisca liturata (Table S2). In these database comparisons, 42 to 52% of transcripts matched database sequences, while the remaining 48 to 58% of transcripts in the D. maidis transcriptome had no database matches at the selected E value cutoffs, such that homologs or putative functions could be assigned to only about half of the transcriptome. This is comparable to findings for transcriptomes of the leafhopper G. nigrifrons, for which approximately 42% of the transcripts had a significant hit to the Swiss-Prot database (E value of <10−6) (28), and the planthopper P. maidis, for which 49% of the sequences had significant matches to known protein sequences in GenBank (E value of <10−3) (32).

TABLE 1

Number of the total 56,112 Dalbulus maidis transcripts matching sequences from two sequence databases

Top hit sequence classification	NCBI protein and nucleotide nr databases (E values of ≤10−2)a			Insect genome with RefSeq database (E values of ≤10−5)b
	No. transcripts	No. species	%c	No. transcripts	No. species	%
Eukaryotes	28,368	1,046	50.6	27,074	279	48.2
Insects	23,001	382	41.0	26,685	211	47.6
Animals (noninsect)	4,062	420	7.2	226	27	0.4
Plants	979	103	1.7	121	27	0.2
Fungi	116	86	0.2	42	14	0.1
Bacteria	767	317	1.4	120	34	0.2
Viruses	54	28	0.1	3	3	0.0
Archaea	19	10	0.0	0	0	0.0
Others	210	55	0.4	80	2	0.1
Total	29,208	1,401	52.0	27,277	318	48.6

Transcripts were blastx and blastn searched against NCBI protein and nucleotide nr databases with E values of ≤10−2. Transcripts were first assigned matches based on blastx results. Transcripts with no blastx matches were then assigned matches based on blastn. The database was downloaded on 10 October 2020.

Transcripts were blastx searched against insect genome with RefSeq database with E values of ≤10−5. The database was retrieved from NCBI GenBank on 18 September 2018 and contained all insect nucleotide nr plus the entire Reference Sequence (RefSeq) collection, a curated, nonredundant set of sequences from >1,400 diverse organisms from genomic DNA, transcript (RNA), and protein sequences (31).

Percentage (%): the number of D. maidis transcripts that match each organism’s sequences divided by total D. maidis transcripts (56,116).

Comparison of D. maidis transcriptome with that of other hopper species.

The transcriptome of D. maidis was compared to those available for the leafhopper, G. nigrifrons, and two planthoppers, P. maidis and S. furcifera. Transcriptome-wide analyses showed that the translated sequences of D. maidis had more orthologs in common with G. nigrifrons than with either of the two planthoppers. Of D. maidis transcript predicted translations, 65% (36,659) had a conserved ortholog in G. nigrifrons compared with 21% (11,460) and 23% (12,914) of transcripts with identified orthologs in P. maidis and S. furcifera, respectively (E value less than 10−5; Table 2). Of the identified leafhopper orthologs, 16.8% (9,439) of translated D. maidis sequences were nearly identical to their G. nigrifrons ortholog (E value of ∼0). In contrast, highly conserved orthologs were identified for only 1.8% (1,037) of S. furcifera orthologs, and none were identified among the P. maidis orthologs (Table 2). These results correspond well with the taxonomic classification of these insects; D. maidis and G. nigrifrons are classified in the same superfamily, Cicadoidea, while S. furcifera and P. maidis taxonomically belong to the superfamily Fulgoroidea.

TABLE 2

Comparison of Dalbulus maidis transcriptome with those of other hopper species

E value cutoff	Percent match (transcript no.)a with:
	Graminella nigrifrons	Peregrinus maidis	Sogatella furcifera
10−5	65.3 (36,659)	20.7 (11,460)	23.0 (12,914)
10−10	60.4 (33,893)	14.7 (8,122)	17.7 (9,932)
∼0	16.8 (9,439)	0.0 (0)	1.8 (1,037)

Percentage of transcripts with blastx hits against transcriptome database of each hopper species, Graminella nigrifrons (28), Peregrinus maidis (32), and Sogatella furcifera (45), divided by the total transcripts of D. maidis at each E value cutoff.

Potential endosymbiont and viral sequences within the D. maidis transcriptome.

Among transcriptome sequences were those with high sequence identity to noninsect sequences, with top hits to other animals, plants, fungi, bacteria, archaea, and viruses (Table 1). Our findings provide data for further studies of endosymbiont bacteria and viruses hosted by D. maidis (Table S3). Sequences with top hits to known bacterial endosymbionts were identified (E value < 10−5), including 56 transcripts with matches to 11 Wolbachia spp., 5 transcripts matching 4 Rickettsia spp., 7 transcripts matching a Sodalis-like symbiont of the spittlebug Philaenus spumarius, 1 transcript most similar to a Burkholderia endosymbiont of Escarpia laminate, 2 transcripts most similar to an endosymbiont of Donacia bicoloricornis, and 2 transcripts most similar to a putative symbiont of Hydra magnipapillata (Table S3). Leafhoppers have been reported to contain facultative endosymbionts, including Wolbachia, Rickettsia, Burkholderia, and Diplorickettsia (33). We identified sequences with top hits to known viruses with an E value of <10−5, including transcripts of MRFV and transcripts with top hits to invertebrate-infecting viruses, including iflaviruses Graminella nigrifrons virus 1 and Scaphoideus titanus iflavirus 1 (Table S3). A wide range of viruses in different families, including Reoviridae, Iridoviridae, Dicistroviridae, and Iflaviridae were previously reported in leafhoppers (34).

Immunity genes and immune signaling pathways in D. maidis.

To characterize immune response genes and signaling pathways in D. maidis, the D. maidis transcriptome was mined for orthologs of known insect immune pathway genes using a data set of 12,255 sequences created using GenBank insect sequences. Of the 56,116 D. maidis transcripts, 3,679 had orthologs in this immunity gene data set as determined using BLASTX (E value < 10−5; Table S4). Among the identified transcripts were those for Toll, Imd, Jak/Stat, prophenoloxidase-activating system, and melanization immune signaling pathways as well as several antimicrobial peptides (AMP) and other immune effectors (Fig. 1; Table S4). Many immune recognition-associated genes were also identified: peptidoglycan recognition proteins (PGRPs), C-type lectins (CTLMA2 protein), hemomucin, bacteria-responsive protein 1, acetylcholine receptor, and a V-type proton ATPase. Eighteen types of PGRP, including LB, LC, LF, SA, SB1, SB2, SE, SC, SCa1, SC2, and SD were identified in D. maidis (Fig. 1; Table S4).

Searching for orthologs of genes in the Toll signaling pathway, important in response to fungi and G+ bacteria and activated by a serine protease and Spaetzle but inhibited by serpin (Fig. 1) (10), we found 9 transcripts encoding putative serine proteases (with top hit matches to serine protease 14, immune-responsive serine protease-related protein [ISPR9], CLIPA14, and CLIPB1, a serine protease immune response integrator), 6 transcripts encoding putative serpins (with top hit matches to immune-responsive serpin-related protein I, serpin 9 inhibitory serine protease inhibitor, serpin 12, SRPN10 protein, SRPN2 protein, noninhibitory serine protease inhibitor, serine protease inhibitor [27A, 88Ea]), 4 transcripts encoding putative Spaetzle-like proteins (with top hit matches to Spaetzle, Spn43Ac, SPZ3, SRPN10, and SRPN2), and 9 transcripts encoding putative Toll proteins (with top hit matches to Toll 6, Toll 7, Toll 8, Toll 9, and Toll 10) (Fig. 1; Table S4). In addition, putative Toll signaling pathway genes encoding weckle and embryonic polarity protein dorsal were identified in the D. maidis transcriptome (Table S4). However, we did not find transcripts encoding orthologs of Tube, Pell, or Cactus in the Toll pathway.

We further searched the D. maidis transcriptome for orthologs of genes in the G− bacteria-responsive Imd and JNK signaling pathways, including PGRPs, IAP, TAK, NF-kappa-B, and relish for synthesizing AMP in the nucleus (Fig. 1) (12). Transcripts encoding putative orthologs of IAP2 protein, TAK1-binding protein, death-associated inhibitor of apoptosis 2, I-kappaB kinase beta, NF-kappa-B essential modulator, NF-κB, mitogen-activated protein kinase kinase kinase 7, and relish were all identified in D. maidis (Fig. 1; Table S4). Transcripts orthologous to genes encoding proteins associated with JNK pathways, including protein kinase JNK, DJNK, and putative mitogen-activated protein kinase ERK, were also identified in D. maidis.

Of particular interest to us was the Jak/Stat signaling pathway, typically associated with virus responses. Viruses are recognized by TEP/Mcr and then activate Jak and Stat for the synthesis of AMPs. Transcripts encoding proteins associated with Jak/Stat signaling pathway, such as immune deficiency signal transducer and transcription activator (Stat), ubiquitin-conjugating enzyme E2 K, TEP1, thioester containing protein (I, II, IV), and ubiquitin carboxyl-terminal hydrolase (Usp2, Usp36) were identified in D. maidis.

Transcripts with highest sequence identity to genes encoding prophenoloxidase-activating system proteins and melanization-associated proteins such as prophenol oxidase subunit 2, melanization protease 3, small heat shock protein, and heat shock protein 90 were identified in the D. maidis transcriptome. In addition, transcripts with identity to genes encoding antimicrobial peptides and effectors such as defensin, defense protein l(2)34Fc, eater, dicer 2, and NO synthases were identified (Table S4).

Orthologs of transcripts encoding proteins that may be associated with cellular immune response or other immune functions were also identified in D. maidis, including the cellular immune response proteins Down syndrome adhesion, Down syndrome, and integrin alpha-PS4 for cell adhesion (35), catenin alpha for cell-cell cohesion (36), and protein kinase C for cell proliferation and regulation (37). Orthologs of transcripts encoding proteins for suppressing the immune response, such as Argonaute 2, autophagy protein 5-like protein, dual specificity protein phosphatase MPK-4, and protein Diedel were identified in D. maidis. Argonaute 2 mediates specific antiviral immunity (38), autophagy protein 5-like protein can negatively regulate the innate antiviral immune response, dual specificity protein phosphatase MPK-4 suppresses bsk/JNK activation during the immune response, and protein Diedel suppresses the IMD pathway in Drosophila (39). Other genes encoding proteins such as draper splice (required in engulfing follicle cells and activating the JNK pathway), ubiquitin carboxyl-terminal hydrolase (playing a role in the innate immune defense against viruses by stabilizing the viral DNA sensor and inhibiting its autophagic degradation), UDP glycosyltransferase precursor Dorothy (synthesizing most of the identified glycolipids, the important antigens in early pathogen infections) (40), vacuolar protein sorting-associated protein 16B (essential for phagosome maturation and the innate immune response to bacteria), and tripeptidyl-peptidase 2 (linking intracellular amino acids homeostasis with the immune response through glycolytic pathway) (41) were also identified in D. maidis (Table S4).

RNAi pathway genes.

RNA interference (RNAi) in insects is triggered by double-stranded RNA (dsRNA) and can play an important role in insect antiviral responses, as well as mortality or interference with pathogen transmission (23, 24, 42). Forty transcripts were identified as RNAi pathway genes by blastx against an RNAi gene data set constructed from sequences of D. melanogaster, A. gambiae, T. castaneum, and G. nigrifrons (E value < 10−10; Table S5). Several transcripts with identity to short interfering RNA (siRNA) pathway, microRNA (miRNA) pathway, and Piwi-interacting RNAs (piRNA) pathway genes were identified in D. maidis (43). We found genes encoding R2D2, the cofactor of Dicer-2 and Argonaute-2 for the siRNA pathway, Argonaute-1, Argonaute-3, and R3D1 for miRNA pathway, and Aubergine, Piwi, and Maelstrom for the piRNA pathway. However, some genes for the miRNA pathway, such as genes encoding Drosha and Loquacious (43), were not identified in D. maidis. Orthologs of RNA interference-deficient (SID) transmembrane channel-mediated proteins, which were discovered in Caenorhabditis elegans and thought to be involved in dsRNA uptake (44), were identified in D. maidis (Table S5).

Gene regulation of D. maidis in response to MRFV.

Gene regulation was compared between MRFV-exposed and naive leafhoppers at two time points at 25°C. Statistical analysis (performed by R package DESeq2) identified 72 transcripts with significantly differential expression, with either relative upregulation or downregulation observed between virus-exposed and naive leafhopper treatments at 4 h postexposure, and 67 at 7 days (d) postexposure (adjusted P value (adj. P) of ≤0.1). These time points were selected to assess responses to immediate virus exposure/entry (4 h) and response to active virus infection and replication at the earliest end of the latent period (7 d). Among the differentially expressed (DE) transcripts, 51 at 4 h (Table 3) and 23 at 7 d (Table 4) shared high identity (E value < 10−5) to an insect sequence ortholog in the NCBI nr database to suggest sequence name and/or functional annotation, while the rest with no matches were assigned no sequence name or functional annotation. No virus-responsive genes were shared between the 4 h and 7 d time points (Tables 3 and 4). Differentially expressed genes with ortholog matches were classified with putative cellular functions in immunity, protein synthesis, transcription, cell rescue, defense, virulence, protein fate, metabolism, cellular communication signaling, and other processes (Fig. 2). Over time, more genes were highly upregulated in response to virus, such that more gene upregulation was observed at the 7 d than at the 4 h time point (Fig. 3; Tables 3 and 4).

TABLE 3

Significantly up- or downregulated transcripts in response to MRFV exposure in D. maidis at 4 h (adj. P ≤ 0.1) based on R package DESeq2 analysis

Contig	Functional category/protein namea	Top hit taxon	E value	LFC	P adj. (R-stat)b
	Immunity
14129	Gram-negative bacteria binding protein 1	Drosophila melanogaster	1.9E−46	0.947	6E−07
3972	Serine/threonine-protein phosphatase 4 regulatory subunit 1-like isoform X4	Thrips palmi	3.0E−143	0.668	8E−03
4144	Inter-alpha-trypsin inhibitor heavy chain H3-like	Sipha flava	4.2E−138	0.421	2E−02
5916	Secreted serpin protein	Pristhesancus plagipennis	3.8E−112	0.491	5E−02
11603	Monoacylglycerol lipase ABHD12 isoform X4	Halyomorpha halys	8.9E−101	0.182	2E−02
16040	Cathepsin L1	Halyomorpha halys	4.3E−99	0.519	5E−02
16335	Secreted serpin protein	Pristhesancus plagipennis	4.5E−120	0.505	2E−02
52673	Beta-1,3-glucan recognition protein 4a	Anasa tristis	5.8E−79	0.712	2E−02
41208	Peptidoglycan recognition protein LC, partial	Graminella nigrifrons	2.4E−56	0.69	0.09
	Transcription
4184	Protein male-specific lethal-3	Frankliniella occidentalis	9.3E−54	0.171	8E−02
11639	Transcription factor IIIA	Cryptotermes secundus	5.5E−43	0.465	2E−02
14378	G patch domain-containing protein 1 homolog	Frankliniella occidentalis	2.9E−14	0.415	5E−02
16806	Longitudinals lacking protein, isoforms N/O/W/X/Y isoform X20	Onthophagus taurus	2.9E−10	0.559	6E−03
	Cell rescue, defense, and virulence
30485	Heat shock factor-binding protein 1 isoform X1	Cimex lectularius	1.6E−27	0.519	5E−02
4722	Predicted: superoxide dismutase [Cu-Zn]	Polistes canadensis	6.8E−79	0.340	8E−02
7302	Predicted: peroxisomal multifunctional enzyme type 2 isoform X2	Bemisia tabaci	7.1E−105	0.317	4E−02
14274	CDGSH iron-sulfur domain-containing protein 3, mitochondrial	Diachasma alloeum	2.9E−42	0.461	8E−02
48263	Predicted: zinc finger protein 664-like (LOC114325064), transcript variant X3, mRNA	Diabrotica virgifera virgifera	4.6E−09	−0.082	2E−03
	Protein fate
1691	Mitochondrial-processing peptidase subunit beta	Halyomorpha halys	0.0E+00	0.459	9E−02
1818	Glutamate dehydrogenase, mitochondrial	Thrips palmi	0.0E+00	0.306	2E−02
4174	Venom serine carboxypeptidase-like	Halyomorpha halys	6.7E−153	0.439	2E−02
4697	Secreted aspartyl glucosaminidase-like protein	Pristhesancus plagipennis	8.9E−86	0.869	2E−02
6979	Probable serine hydrolase isoform X2	Halyomorpha halys	1.6E−73	0.210	6E−02
	Metabolism
4262	NADP-dependent malic enzyme-like isoform X2	Nilaparvata lugens	1.9E−149	0.401	8E−02
13683	Lipase 3	Zootermopsis nevadensis	5.8E−110	0.563	2E−02
1872	V-type proton ATPase 21 kDa proteolipid subunit	Culicoides sonorensis	1.0E−72	0.869	8E−02
	Protein with binding function or cofactor requirement
1491	NcSP70 mRNA for transferrin	Nephotettix cincticeps	0.0E+00	1.340	8E−02
40821	Aminopeptidase N-like protein	Nilaparvata lugens	6.3E−70	0.758	5E−02
	Cellular communication signal
2029	Growth factor receptor-bound protein 10-like (LOC111865116)	Cryptotermes secundus	2.3E−71	1.229	8E−02
5430	Fatty acid-binding protein, muscle isoform X2	Cryptotermes secundus	1.8E−48	0.735	9E−02
6798	CD63 antigen-like	Nilaparvata lugens	6.1E−80	0.210	2E−02
6907	Predicted: histidine triad nucleotide-binding protein 1	Papilio polytes	3.0E−64	0.381	2E−02
9888	Predicted: inactive dipeptidyl peptidase 10 isoform X2	Bemisia tabaci	0.0E+00	0.526	2E−02
14429	Phosphoesterase 24, partial	Laodelphax striatellus	1.3E−14	0.321	5E−02
	Regulation of metabolism and protein function
6496	Protein-tyrosine sulfotransferase isoform X2	Nilaparvata lugens	2.0E−109	0.456	2E−02
42206	Juvenile hormone acid O-methyltransferase-like	Belonocnema treatae	8.9E−11	0.322	7E−02
	Cellular transport
7433	ATP-binding cassette transporter subfamily C member 4 isoform 1	Laodelphax striatellus	6.5E−111	0.576	4E−02
	Unknown
743	Hypothetical protein LSTR_LSTR005741	Laodelphax striatellus	0.0E+00	0.439	1E−01
1800	Hypothetical protein LSTR_LSTR002520	Laodelphax striatellus	0.0E+00	0.255	9E−02
2450	Hypothetical protein LSTR_LSTR001311	Laodelphax striatellus	1.1E−155	0.434	1E−01
6309	Hypothetical protein LSTR_LSTR008156	Laodelphax striatellus	0.0E+00	0.316	5E−03
6842	Hypothetical protein B7P43_G09429	Cryptotermes secundus	3.1E−141	0.124	2E−02
7351	Hypothetical protein LSTR_LSTR004310	Laodelphax striatellus	1.9E−132	0.407	5E−02
7894	Hypothetical protein B566_EDAN002302	Ephemera danica	2.8E−92	0.579	2E−02
8967	Hypothetical protein LSTR_LSTR005033	Laodelphax striatellus	3.4E−27	0.059	8E−03
10466	Uncharacterized protein LOC110834514	Zootermopsis nevadensis	1.5E−96	0.146	4E−02
11144	Hypothetical protein LSTR_LSTR008826	Laodelphax striatellus	1.7E−45	0.279	8E−02
12616	Uncharacterized protein LOC110828818	Zootermopsis nevadensis	8.7E−10	0.407	8E−02
34166	Hypothetical protein YQE_06069, partial	Dendroctonus ponderosae	5.7E−14	0.271	8E−02
38705	Uncharacterized protein LOC108632082	Ceratina calcarata	1.4E−14	0.649	4E−02
30451	Uncharacterized LOC108365256 (LOC108365256), ncRNA	Rhagoletis zephyria	2.2E−12	0.931	5E−03
	Hits to other organisms
5474	Cathepsin L	Penaeus japonicus	4.9E−131	0.35	0.02
22137	30S ribosomal subunit protein S4 (S9e)	Candidatus Sulcia muelleri	9.3E−100	−0.42	0.02
25497	30S ribosomal subunit protein S4 (S9e)	Candidatus Sulcia muelleri	2.4E−126	−0.32	0.02
38707	Polyprotein	Graminella nigrifrons virus 1	0.0E+00	0.75	0.08
38980	Cathepsin B-like	Daphnia magna	1.2E−16	0.22	0.08
52215	Replicase polyprotein	Maize rayado fino virus	0.0E+00	1.89	0.00
	No hits
12823				−0.30	0.09
18136				0.75	0.08
23542				0.11	0.00
28549				0.30	0.00
31340				0.74	0.05
32068				0.44	0.07
39256				1.24	0.01
39750				0.61	0.08
39990				0.66	0.02
40074				1.03	0.00
41473				−0.22	0.08
42729				0.65	0.00
45919				0.79	0.08
48016				−0.29	0.08
52215				1.89	0.00
55352				0.27	0.00

Assignment of putative functional category was based on Gene Ontology and UniProtKB.

P adj.: adjustment of P value statistic performed by R package and DESeq2 software, obtained by the Wald test and corrected for multiple testing by the Benjamini and Hochberg method.

TABLE 4

Significantly up- or downregulated transcripts in response to MRFV exposure in D. maidis at 7 d (adj. P ≤ 0.1) based on R package, DESeq2 analysis

Contig	Functional category/protein namea	Top hit taxon	E value	LFC	P adj. (R-stat)b
	Immunity
1730	26S proteasome regulatory subunit	Drosophila melanogaster	7.7E−57	2.480	4E−06
1795	1-Phosphatidylinositol 4,5-bisphosphate phosphodiesterase epsilon-1-like (LOC111050143)	Nilaparvata lugens	0.0E+00	3.250	4E−06
4646	Signal peptidase complex catalytic subunit SEC11A (LOC111512019)	Leptinotarsa decemlineata	8.5E−40	1.889	4E−05
15291	Serine protease 14	Anopheles quadriannulatus	1.1E−07	2.653	8E−09
	Protein synthesis
6560	60S ribosomal protein L22 (LOC111047508)	Nilaparvata lugens	2.7E−43	1.575	3E−08
5558	40S ribosomal protein S26 (LOC111056161)	Nilaparvata lugens	8.3E−62	3.238	4E−06
45832	60S ribosomal protein L29-1 (LOC113213475)	Frankliniella occidentalis	1.3E−32	2.276	2E−08
54392	60S ribosomal protein L38 (LOC112051453)	Bicyclus anynana	5.4E−43	2.483	2E−08
	Transcription
1705	Pre-mRNA-splicing factor Slu7 (LOC110833152)	Zootermopsis nevadensis	0.0E+00	1.291	7E−05
17037	Small nuclear ribonucleoprotein E	Agrilus planipennis	3.0E−52	2.950	1E−11
	Cell rescue, defense, and virulence
10958	Cysteine dioxygenase type 1 (LOC106648916)	Trichogramma pretiosum	4.8E−42	1.980	6E−06
	Protein fate
4873	Legumain	Colius striatus	2.5E−42	1.824	4E−06
7611	Gamma-glutamyltranspeptidase 1-like (LOC111053586)	Nilaparvata lugens	1.6E−109	1.758	1E−08
	Metabolism
1693	Protein arginine N-methyltransferase 1-like (LOC111063517)	Nilaparvata lugens	0.0E+00	1.763	2E−11
10321	Long-chain-fatty-acid–CoA ligase 3 (LOC6592464)	Drosophila persimilis	7.8E−95	2.585	2E−06
14111	Clone CFSNI1556 alpha-amylase	Coptotermes formosanus	2.4E−21	2.364	4E−06
	Cellular communication
52876	Predicted: CD63 antigen-like, partial	Atta colombica	1.6E−14	1.597	4E−06
	Regulation of metabolism and protein function
14919	Troponin I protein	Lethocerus indicus	4.8E−06	2.884	2E−05
	Cellular transport
1701	Importin-5 (LOC111064460)	Nilaparvata lugens	0.0E+00	2.085	6E−09
	Unknown
15127	Hypothetical protein LSTR_LSTR004460	Laodelphax striatellus	4.1E−86	1.714	5E−08
55795	Hypothetical protein B566_EDAN012453	Ephemera danica	2.5E−28	2.474	1E−08
12675	TR19232_c0_g1_i2 transcribed RNA sequence	Graphocephala atropunctata	3.1E−23	0.788	3E−03
13940	Uncharacterized LOC111875770	Cryptotermes secundus	1.0E−13	2.169	9E−05
	Hits to other organisms
21724	Hypothetical protein CAPTEDRAFT_190440	Capitella teleta	2.4E−27	0.11	0.05
50276	DDE-type integrase/transposase/recombinase, partial	Gammaproteobacteria bacterium	1.9E−27	1.17	0.00
52215	Replicase polyprotein	Maize rayado fino virus	0.0E+00	3.11	0.00
	E value cut more than 1E−5/no hit
13799				0.64	0.00
17553				2.08	0.00
21090				0.13	0.08
27462				3.22	0.00
29228				0.31	0.01
30671				0.60	0.07
31685				1.69	0.00
33672				1.38	0.00
33790				1.49	0.00
33834				2.01	0.00
34274				0.33	0.07
34415				0.20	0.06
36156				1.40	0.00
36277				1.53	0.00
36707				1.23	0.00
36742				2.81	0.00
37053				0.79	0.00
39273				0.81	0.00
39372				1.18	0.00
39493				1.97	0.00
40542				1.39	0.02
41250				0.41	0.00
41307				−0.15	0.04
41373				2.43	0.00
42400				0.68	0.00
45046				1.45	0.00
45095				1.50	0.00
45390				0.88	0.00
45402				0.22	0.00
45868				0.50	0.03
47633				2.27	0.00
47635				1.97	0.00
47644				3.69	0.00
49354				1.67	0.00
53755				2.09	0.00
53826				1.09	0.00
54113				2.07	0.00
54179				2.52	0.00
54185				NAc	0.00
54402				1.30	0.00
54799				−0.14	0.09

Assignment of putative functional category was based on Gene Ontology and UniProtKB.

P adj.: adjustment of P value statistic performed by R package and DESeq2 software, obtained by the Wald test and corrected for multiple testing by the Benjamini and Hochberg method.

NA, not applicable; naïve expression level is zero, and therefore no calculation of LFC can be made.

FIG 2

The number of the significantly up - or downregulated transcripts based on R package DESeq2 analysis with different functional categories in response to MRFV in Dalbulus maidis at 4 h and 7 d at 25°C. A total of 139 out of 56,116 MRFV-responsive transcripts in D. maidis were assigned to functional categories based on Gene Ontology and UniProtKB.

FIG 3

Heatmaps of the significantly up- or downregulated transcripts in response to MRFV exposure in Dalbulus maidis at 4 h or 7 d based on R package DESeq2 analysis. Differential expression (DE) at 4 h: the left four columns of the 72 genes that are differentially expressed at 4 h; DE at 7 d: the right four columns of 67 genes that are differentially expressed at 7 d. The colors stand for different log10(FC + 1), (i.e., LFC) with deepest blue at −2 and deepest red at 2. MRFV/naive at 4 h and 7 d is the LFC between MRFV-exposed and naive leafhopper at 4 h and 7 d; 7 d/4 h for naive and MRFV is the LFC between 7 d and 4 h for MRFV-exposed and naive leafhopper. Heatmaps were created using R (version 3.6.2, 2019-12-12) statistical software with ggplot2 (56).

Immunity and RNAi pathway genes responsive to MRFV.

Nine immunity-related transcripts [log10(fold change + 1) [LFC]: 0.2 to 0.9] at 4 h and four immunity-related transcripts (LFC: 1.9 to 3.3) at 7 d were differentially expressed in naive and MRFV-exposed D. maidis (R package in DESeq2, adj. P ≤ 0.1; Tables 3 and 4). Those with increased expression had high identity with insect transcripts encoding recognition proteins (Gram-negative bacteria binding protein 1 [GNBP]; beta-1,3-glucan recognition protein 4a) and serpin proteins (inter-alpha-trypsin inhibitor heavy chain H3-like; secreted serpin protein) were significantly upregulated (Tables 3 and 4). No RNAi pathway-related genes were significantly up- or downregulated at either time point (adj. P ≤ 0.1; Tables 3 and 4).

Transcription, translation, and other differentially expressed transcripts in D. maidis in response to MRFV exposure.

Compared to naive D. maidis, MRFV-exposed leafhoppers showed significant upregulation of genes associated with metabolism, energy, transcription, protein synthesis, protein fate, cell cycle, DNA processing, cellular communication/signal transduction mechanism, cell rescue, defense, virulence, cellular transport, transport facilitation, and transport route putative functions assigned based on similarity to genes with previously annotated GO terms via UniProtKB. Transcripts in these categories were even more highly expressed at 7 d than at 4 h at 25°C (see LFC values in Tables 3 and 4).

MRFV genomic RNA can serve as mRNA for viral replication proteins required for viral RNA replication and is sufficient to launch infection in inoculated plants (2, 3). Translation of viral proteins requires host ribosomes (reviewed in reference 4). D. maidis ribosomal proteins, including L22, L29, L38, and S26, were significantly upregulated by 7 d after leafhoppers were exposed to MRFV (P ≤ 0.1; Table 4). RNA virus replication and gene expression utilize a virus-encoded RNA-dependent RNA polymerase. However, genes associated with leafhopper transcription were also significantly upregulated in MRFV-exposed leafhoppers (Table 3 and 4). For example, transcripts encoding the splicing factor Slu7 and small nuclear ribonucleoprotein E were significantly upregulated at 7 d (P ≤ 0.1; Table 4), with high LFC values as well (LFC is 0.81 at 4 h and 1.29 at 7 d for the transcript encoding splicing factor and 1.36 at 4 h up to 2.95 at 7 d for the transcript encoding nuclear ribonucleoprotein).

The gene encoding cysteine dioxygenase, which may be involved in oxidative response (LFC: 1.98), was significantly upregulated at 7 d (P ≤ 0.1) (Table 4). Also, MRFV significantly triggered upregulation of genes likely involved in the oxidative response, protein folding, and protein fate, including genes encoding heat shock factor-binding protein 1 isoform X1 (LFC: 0.52), zinc finger protein 664-like (LFC: 0.46), and superoxide dismutase [Cu-Zn] (LFC: 0.34; Tables 3 and 4). Heat shock factor-binding protein 1 isoform X1 can negatively regulate the heat shock response, which may have a role in the suppression of the activation of the stress response. Zinc finger protein triggers the expression of H2O2 and superoxide, important reactive oxygen species. However, the gene for superoxide dismutase [Cu-Zn] was also significantly upregulated after leafhopper exposure to MRFV, which may degrade reactive oxygen species and will be beneficial for the virus propagation.

Genes associated with metabolism proteins such as amylase and ligase were significantly upregulated (adj. P ≤ 0.1), especially at 7 d compared with 4 h (Tables 3 and 4). Genes associated with proteins functioning in cellular transport, transport facilitation, transport route, cellular communication, and signal transduction mechanism were also significantly upregulated at 4 h and 7 d at 25°C (Tables 3 and 4). In addition, transcripts with unknown function (14 at 4 h and 4 at 7 d) or blast hits to other organisms (7 at 4 h and 3 at 7 d) or without hits (16 at 4 h and 41 at 7 d) were also significantly up- or downregulated (adj. P ≤ 0.1; Tables 3 and 4).

Comparison of MRFV-responsive transcripts in D. maidis with those identified in other virus-vector systems.

Differentially expressed transcripts of MRFV-exposed D. maidis were compared with those identified by various methodologies for other hopper species in response to the propagative plant viruses they transmit. D. maidis shared 15 MRFV-responsive transcripts with orthologs in S. furcifera differentially expressed in response to southern rice black-streaked dwarf virus (SRBSDV; of 4,611 determined by an absolute value of the log2 ratio of >1) (45). These 15 transcripts were associated with putative functions in protein synthesis (40S ribosomal protein S26, 60S ribosomal protein L38, 60S ribosomal protein L29), transcription (small nuclear ribonucleoprotein E), oxidative response (superoxide dismutase [Cn-Zn], peroxisomal multifunctional enzyme type), immunity (26S proteasome regulatory subunit), cellular communication signaling (phosphoesterase 24, inactive dipeptidyl peptidase 10 isoform X2, histidine triad nucleotide-binding protein 1, growth factor receptor-bound protein 10-like), protein fate (venom serine carboxypeptidase-like), metabolism (V-type proton ATPase 21 kDa proteolipid subunit), and unknown functions (two transcripts). One D. maidis MRFV-responsive transcript shared an ortholog among the 891 G. nigrifrons maize fine streak virus-responsive transcripts to MFSV, determined by Bayes-moderated t tests with a P value of <0.05 (28), which did not match any sequence in the NCBI nr database. One differentially expressed transcript was shared between MRFV-responsive D. maidis transcripts and maize mosaic virus (MMV)-responsive P. maidis transcripts (out of 144 determined as statistically significant using the Benjamini-Hochberg procedure for multiple tests with an adjusted P value cutoff of 0.1) (32) (Fig. 4; Table S6), which also had no NCBI nr matches. The larger number of matches of differentially expressed transcripts between virus-responsive D. maidis and S. furcifera than that of matches between D. maidis and G. nigrifrons and between D. maidis and P. maidis is likely due to the larger pool of differentially expressed transcripts reported for S. furcifera (Fig. 4; Table S6). No virus-responsive transcript orthologs were shared among all four virus-responsive vectors, nor were any virus-responsive transcripts shared among any set of three vector species that included D. maidis (Fig. 4; Table S6).

FIG 4

Venn diagram of shared virus-responsive transcripts among four hopper species. Venn diagram indicating the orthologous matches of Dalbulus maidis MRFV-responsive transcripts (139 in total from both time points with adj. P of ≤0.1), compared with propagative virus-responsive transcripts from three other hopper vector systems: Graminella nigrifrons to maize fine streak virus (MFSV) (28), 873 transcripts, Sogatella furcifera to southern rice black-streaked dwarf virus (SRBSDV) (45), 4,611 transcripts, and Peregrinus maidis to maize mosaic virus (MMV) (32), 144 transcripts, as identified using methodologies reported in each publication.

Validation of RNA-Seq using RT-qPCR and correlation of gene expression among replicates by PCA.

Twenty MRFV-responsive D. maidis transcripts belonging to different functional categories were selected for further response validation by quantitative real-time PCR (RT-qPCR). Of these, 13 transcripts showed differential expression (adj. P ≤ 0.1) or high up- or downregulation (LFC of ≤−0.78 or ≥0.78, corresponding to a 5-fold down- or upregulation) in response to MRFV exposure, and 7 transcripts had relatively high abundance in both naive and MRFV-exposed leafhoppers (Fig. 4). We selected these transcripts for validation based on two criteria: (i) the transcripts had a specific target abundant enough for RT-qPCR amplification with expression levels of >0.01 reads per kilobase per million (RPKM) (Tables 3 and 4), and the primer efficiencies were close to 2.0 (approximately ±0.1; efficiency defined here as E = 10−1/slope from the slope of the dilution series C values; https://www.bio-rad.com/webroot/web/pdf/lsr/literature/Bulletin_5279.pdf; Table S8), and (ii) the transcripts were annotated with a putative function related to the interaction between plant viruses and vector insects, such as immunity, oxidative response, protein synthesis, transcription, or metabolism (including Toll, ribosomal proteins, small nuclear ribonucleoprotein, heat shock protein, succinate-CoA ligase). After normalizing the RT-qPCR expression level using three reference genes, the relative expression (RE) and fold change between MRFV-exposed and naive leafhopper were calculated. The results showed that for most genes, the fold differences between MRFV-exposed and naive leafhoppers derived from RT-qPCR comparison matched well with those observed by Illumina transcript analyses (Fig. 5). Spearman’s correlation analysis found a significantly positive correlation (R = 0.835, P ≤ 0.0001) between the fold change for the relative expression of MRFV-exposed and naive in RT-qPCR with fold change in RNA-Seq, indicating that transcriptional changes measured from RNA-Seq data reasonably describe leafhopper responses to MRFV exposure within this study. Principal component analysis (PCA) showed transcript expression level variation between the replicates of each treatment and time point but also demonstrated correlation of biological replicates by clustering of datapoints by treatment (MRFV-exposed and naive) and postexposure time point (4 h or 7 d; Fig. 6). Two principal components accounted for 69.4% of the total variation, with most variation (60.8%) explained by PC1 showing clear separation of transcript expression by time point and the second most variation (8.6%) explained by PC2 showing treatment differences (Fig. 6).

FIG 5

The comparative fold change of gene expression between sequence analysis and RT-qPCR. (A) Fold change of MRFV-exposed versus naive Dalbulus maidis in response to MRFV exposure at 7 d at 25°C by sequence analysis. (B) The validation of gene expression by RT-qPCR. REMRFV and REnaïve are the relative gene expression of MRFV-exposed versus naive D. maidis in response to MRFV exposure at 7 d at 25°C by RT-qPCR. Note that there is a significantly positive correlation (R = 0.835, P ≤ 0.0001) between the fold change for the relative expression of MRFV-exposed and naive in RT-qPCR with fold change in RNA-Seq analyzed by Spearman’s correlation.

FIG 6

Principal-component analysis (PCA) of correlation of Dalbulus maidis gene expression levels between treatments, time points, and biological replicates. Each colored circle represents a single treatment replicate, and the same-colored circles represent the same treatment. EL1, EL2, EL3, and EL4 are the expression levels (EL) of four different replicates for naive leafhopper samplings at 4 h, EL5, EL6, EL7, and EL8 for MRFV-exposed leafhopper samplings at 4 h, EL9, EL10, EL11, and EL12 for naive leafhoppers at 7 d, EL13, EL14, EL15, and EL16 for MRFV-exposed leafhoppers at 7 d. A total of 56,116 expression levels for each replicate were analyzed. Figure created using R (version 3.6.2, 2019-12-12) statistical software with ggplot2 (56).

DISCUSSION

We report here the first transcriptome for the corn leafhopper, D. maidis, and the transcriptional response to exposure to a plant-infecting virus that it both transmits and is infected by. We found 139 total MRFV-responsive transcripts at 4 h and 7 d postexposure, assigned to several functional categories. These MRFV-responsive genes may be either reactionary and virus-limiting or beneficial to virus invasion and propagation.

Our results showed that D. maidis contains expected insect immunity and pathogen-responsive pathways, including Toll, Imd, and Jak/Stat immune signaling pathways, the prophenoloxidase-activating system, and several RNAi pathways (10–14). It is unclear whether the few pathway genes not identified in our transcriptome assembly reflect actual absence of these orthologs or incompleteness of the transcriptome.

Immune recognition genes, such as those that encode Gram-negative bacteria binding proteins (GNBP) and peptidoglycan recognition proteins (PGRP), were significantly upregulated in response to MRFV at 4 h, but this response was not durable to 7 d post MRFV exposure. Immune recognition genes are the first responders to pathogen invasion to trigger insect immune responses (8). On the other hand, genes encoding proteins expected to inhibit the immune response, which might contribute to a benign equilibrium with the virus, were not identified as differentially regulated. However, mean LFC values of immune response inhibitors such as Spn43Ac, serine protease inhibitor 88Ea, SRPN10 protein (46), and protein D2-like (35) were much higher in MRFV-exposed than naive D. maidis (LFC ≥ 0.78, i.e., fold change [FC] ≥ 5.0) at 7 d at 25°C. Serpins negatively regulate innate immunity by inhibiting serine proteases that trigger insect immune responses. In addition, several studies have revealed that insect serpins could also possess direct anti-pathogen activity which, upon infection, can inhibit immune response (10, 11, 13). These genes and associated pathways may be responsible for dampening the immune response triggered in the vector by 7 d postexposure, during virus propagation. Surprisingly, RNAi pathway genes were not highly responsive to MRFV at either time point, suggesting that MRFV may evade this antiviral response.

MRFV is a positive-strand RNA (+ssRNA) persistent propagative virus. The translation of viral proteins requires host ribosomes for MRFV propagation inside the vector (4). Genes associated with protein synthesis, including those encoding ribosomal proteins L22, S26, L29, and L38, were significantly upregulated at 7 d at 25°C. Upregulation of these genes could benefit MRFV propagation. Although transcription and replication of MRFV are expected to utilize virus-encoded RNA-dependent RNA polymerase as for other RNA viruses (47), transcription-associated genes encoding proteins such as the splicing factor and nuclear ribonucleoprotein were found significantly upregulated in D. maidis. Increased transcription machinery presumably supports additional transcription of leafhopper sequences responsive to virus exposure. The silencing of nuclear ribonucleoprotein resulted in a 90% reduction of viral RNA replication in cells infected with the Japanese encephalitis virus (47). Genes encoding proteins associated with metabolism, such as arginine N-methyltransferase 1-like, long-chain-fatty-acid–CoA ligase 3, and alpha-amylase, were also significantly upregulated, suggestive of major alterations to cellular metabolites and energy expenditures in response to MRFV. MRFV propagation requires and utilizes energy and metabolic resources from its hosts; thus, these changes are likely virus-supportive globally if not individually. Indeed, another +RNA virus, tomato bushy stunt virus, effectively exploits and hijacks host energy-producing enzymes to provide a readily available source of ATP for the viral replication process (48).

Oxidative and stress responses were induced by leafhopper exposure to MRFV. Genes involved in the oxidative and stress responses, such as those encoding zinc finger protein, were significantly upregulated. Zinc finger protein was reported to produce H2O2 and can trigger production of superoxide in the plant Arabidopsis thaliana and the insect Drosophila melanogaster (49). The high expression of these genes could trigger the leafhopper oxidative and stress response and produce reactive oxygen species (ROS). Thus, these genes may be reactionary and virus limiting. However, this production could be degraded by superoxide dismutase [Cu-Zn]. Superoxide dismutases are key antioxidant enzymes converting superoxide into oxygen and hydrogen peroxide (50). We found that superoxide dismutase was significantly upregulated after leafhopper exposure to MRFV at 4 h, suggesting that ROS production may be disarmed by 4 h after virus entrance vector. A superoxide dismutase was also found to be upregulated in G. nigrifrons in response to MFSV (28). Thus, it may be a common response of leafhopper vectors to propagative plant viruses.

Unsurprisingly, the D. maidis transcriptome was most similar to that of the leafhopper G. nigrifrons and less similar to those of the planthoppers P. maidis and S. furcifera. In hopes of identifying possible shared response pathways of leafhopper and planthopper vectors, we compared differentially expressed transcripts of MRFV-exposed D. maidis with those of other hopper species in response to the plant viruses they transmit in a propagative manner. Though this meta-analysis was imperfect because different methods were applied in each system, results did indicate some potential shared responses. More virus-responsive orthologs were shared between D. maidis and S. furcifera in response to SRBSDV (45) than between other pairwise comparisons, probably because of the larger number of transcripts for comparison than that of G. nigrifrons in response to MFSV (28) and P. maidis in response to MMV (32), which each shared a single virus-responsive ortholog of unknown function. These shared responsive genes had putative roles in protein synthesis, transcription, oxidative response, immunity, cellular communication signaling, protein fate, and metabolism. Together, the D. maidis transcriptome assembly and identification of MRFV-responsive genes provide insights into the molecular interaction of this virus with its leafhopper host and vector and, more broadly, into the interactions between hopper vectors and the viruses they propagate and transmit.

MATERIALS AND METHODS

Leafhopper colony and virus maintenance.

D. maidis was collected by Charles Summers, UC Riverside from Kings County, California, and maintained on maize cv. Early Sunglow (Schlessman Seed Co.) seedlings in 30 cm by 38 cm by 50 cm cages in an insect rearing room with 15 light (L)/9 dark photoperiod and 27°C (light) to 16°C (dark) temperatures. MRFV was collected by Lowell R. Nault’s lab from Harlington, Texas, in the Rio Grande Valley (51) and was maintained in maize cv. Early Sunglow by serial inoculation with leafhoppers. Virus presence in experimental plants was verified using enzyme-linked immunosorbent assays as described previously with virion antisera prepared by Donald T. Gordon at The Ohio State University (see reference 52).

MRFV-exposed and naive leafhopper sample preparation.

As described for previous experiments with G. nigrifrons (28), batches of 800 newly emerged D. maidis virgin females (at or less than 24 h after molt from 5th instar nymph stage) were collected, starved for 3 h, and divided evenly (200 insects/cage) between four 0.1 cubic meter cages. Each cage contained 49 3.5- to 4-week-old maize cv. Early Sunglow plants, with two of the cages containing MRFV-infected plants (maize cv. Early Sunglow prepared by serial inoculation of MRFV with leafhoppers and showing typical MRFV symptoms) and two cages containing healthy plants. Two hundred leafhoppers were released into each cage, and leafhoppers were allowed to feed on plants in cages in a growth chamber with a photoperiod of 15 L/9 D at ∼1,757 to 2,033 lumens/ft2. One cage for each of the two treatments (infected versus healthy plant exposure) was placed at 25°C and another pair at 30°C in each of four experimental replicates. Samples at 30°C were used for transcriptome generation but not for identification of differentially expressed transcripts. Twenty-five leafhoppers were randomly collected from each cage at each of two time points, 4 h and 7 d after placement on plants, from each of two separate cages for each time point. Collected leafhoppers were snap-frozen in liquid nitrogen and stored at −80°C until RNA isolation in pools of 25 insects/treatment (counted before freezing). Each treatment was replicated four times, resulting in a total of 32 samples (2 treatments by 2 temperatures by 2 time points by 4 experimental replicates).

RNA isolation and cDNA library preparation.

Total RNA was extracted from pools of 25 leafhoppers using the RNeasy minikit (Qiagen, Germantown, MD). RNA quantity and quality were assessed using the Experion Automated Electrophoresis system (Bio-Rad, Hercules, CA). For each sample, 200 ng total RNA was input into adaptor-ligated cDNA libraries generated from poly(A)-selected RNAs (captured by oligo-dT beads) at the Ohio State University Comprehensive Cancer Center for RNA-Seq using the TruSeq stranded mRNA library prep kit (Illumina, San Diego, CA) following the manufacturer’s protocol. Quantification and quality inspection of double-stranded cDNA (ds-cDNA) was again assessed by Bio-Rad Electrophoresis System and Qubit 2.0 fluorometer (Life Technologies, Carlsbad, CA). Samples were diluted to 20 nM to generate the multiplexed cDNA libraries from 32 adaptor-tagged sample pools.

Sequence assembly and annotation.

Sequences were generated using 50 fmoles of cDNA library prepared above on one flow cell lane in the Illumina HiSeq 2000 platform at the Ohio State University Comprehensive Cancer Center, for 100-nt paired-end reads. The Illumina raw data were then trimmed of adapter indexes, primers, and poly(A) tails with the default settings (ambiguity limit of 2 and quality limit of 0.05) in CLC Genomics Workbench v.11 (Qiagen Bioinformatics, Redwood City, CA), and trimmed reads of less than 20 nt were discarded. The contigs were obtained by de novo assembly with default settings (mismatch cost 2, insertion cost 3, deletion cost 3, length fraction 0.5, stability fraction 0.8) using CLC Genomics Workbench v.11. To further assemble de novo-assembled contigs into transcripts, all the contigs from 32 samples were assembled in Sequencer 5.4.6 (Gene Codes Corporation, Ann Arbor, WI) with the default assembly parameters (85% minimum match percentage and 20% minimum overlap). This collection of assembled contigs (56,116 transcripts) was used as the D. maidis reference transcriptome for expression and sequence analyses. The reference transcriptome sequences were blastx and megablast searched against the standard NCBI nonredundant protein and nucleotide databases on 12 June 2020 and also blastn searched against a database (31.5 GB) combining all insect genomic sequences and all RefSeq sequences downloaded from NCBI GenBank on 18 September 2018. The combined insect nucleotide nr plus RefSeq database included more than 9 million insect nucleotide sequences and over one million sequences from over 2,400 organisms, including archaea, bacteria, eukaryotes, and viruses, as the reference sequences (RefSeq) database (30). Illumina sequence reads were deposited in the NCBI short sequence read (SRA) archive as BioProject PRJNA579843.

Estimation of BUSCO scores.

Benchmarking Universal Single-Copy Ortholog (BUSCO) from OrthoDB (https://orthodb.org), used to assess the completeness of genome assemblies, annotated gene sets, and transcriptomes in terms of expected gene content (29, 30, 53), was applied. To assess the completeness of D. maidis transcriptome, BUSCO assessment was completed by comparing the D. maidis transcriptome to insect single-copy orthologs (insecta_odb10.2019-11-20, which contains 367 conserved genes for insects) using the BUSCO v4.1.1_cv1, which was developed by Swiss Institute of Bioinformatics (29). Here, the recovered matches were classified as complete (C) if their lengths were within the expectation of the BUSCO profile match lengths (default settings for the container version busco_v4.1.1_cv1). Alternatively, the matches were classified as fragmented (F) if their lengths were only partially recovered; if no matches passed the tests of orthology, BUSCO groups were classified as missing (M) (29, 30).

Comparison of D. maidis transcriptome with other hopper species.

Transcriptomes for G. nigrifrons (28), S. furcifera (45), and P. maidis (32) were kindly provided by the respective authors and were used as transcriptomic databases for each species. D. maidis transcripts (56,116) were blastn searched against each database, and the matches with E values lower than 10−5 were selected for analysis of conserved orthologs among hopper species.

Identification of immune response and RNAi genes.

Genes with known or putative immune and RNAi functions for D. melanogaster, A. gambiae, T. castaneum, and G. nigrifrons were retrieved from NCBI proteins database on 1 May 2019 to create an immune response gene database and an RNAi gene database containing 12,255 and 1,411 genes, respectively. D. maidis transcripts matching with E values less than 10−5 were completed using CLC genomics and selected as probable immune response and RNAi genes. In addition, some RNAi target genes were assigned to other functional categories. Thus, we removed them from the list and kept only the RNAi pathway genes.

Differential gene expression analyses, heatmapping, and PCA.

Quality-controlled Illumina sequence reads (adapter indexes, poly(A) tails, and low-quality sequence trimmed) from each of the 32 samples were mapped to the D. maidis reference transcriptome assembly (56,116 transcripts) using CLC Genomics Workbench v.11 with default settings (match score = 1, mismatch cost = 2, insertion cost = 3, deletion cost = 3, length fraction = 0.5, and similarity fraction = 0.8) for differential transcript expression analysis. Transcript expression for each sample was determined from the mapped read and total read counts data for each of four replicates per treatment. Mapped transcript read counts (paired only) were normalized by calculating the number of unique reads (i.e., mapping to only one transcript in the D. maidis transcriptome) per kilobase per million mapped reads (RPKM) as described previously (54). Transcript expression levels were calculated using the following formula: expression level = total reads mapped to the transcript/total mapped reads in millions for the sample/transcript length in kilobases. The expression fold change (FC) of transcripts in leafhoppers exposed to MRFV compared to those of naive leafhoppers was calculated for each transcript as follows, averaging expression values for four samples per treatment: (average expression level for MRFV-exposed samples – average expression level for naive samples)/average expression level for naive samples. In order to normalize the expression fold change for both up- and downregulated genes, log10(FC + 1) (LFC, in short) was calculated. Transcripts were defined as highly upregulated if LFC was ≥0.78 (i.e., FC ≥ 5.0) or highly downregulated if LFC was ≤−0.78 (i.e., FC ≤ −0.83). The expression level and FC values between MRFV-exposed and naive leafhopper transcripts for 56,116 total transcripts are listed in Table S7. Statistical tests comparing MRFV-exposed and naive leafhopper transcript read counts were performed by R statistical software version 4.0.2 (R Core Team; https://R-project.org) and DESeq2 package version 1.28.1 (55). Adjusted P values (adj. P) were determined by the Benjamini and Hochberg method, and transcripts with an adj. P of ≤0.1 were determined to be significantly differentially expressed.

Heatmaps were created using LFC comparing MRFV-exposed and naive leafhopper transcripts at 4 h and 7 d or comparing 7 d and 4 h for MRFV-exposed and naive leafhoppers using R statistical software (version 3.6.2, 12 December 2019) with ggplot2 (56). Principal-component analysis (PCA) was carried out using D. maidis gene expression levels between 4 biological replicates of 4 treatments by R statistical software with ggplot2.

Comparison of D. maidis transcripts with those of other hopper species infected with propagative plant viruses.

Differential transcript expression data for virus exposure of G. nigrifrons to maize fine streak virus (MFSV) (of 873, determined by Bayes-moderated t tests with a P value of <0.05 [28]), S. furcifera to southern rice black-streaked dwarf virus (SRBSDV) (of 4,611 determined by a false discovery rate [FDR] of <0.001 [45]), and P. maidis to maize mosaic virus (MMV) (of 144 determined by Benjamini-Hochberg procedure with an adjusted P value cutoff of 0.1 [32]) were kindly provided by the respective authors. These sequences (873, 4,611, and 144, respectively) were used as a virus-responsive transcript database for each hopper species. D. maidis transcripts were blastn searched against each database, and matches with E values lower than 10−5 were selected.

RT-qPCR for validation of gene expression calculations.

To validate observed patterns of gene expression, we used quantitative real-time PCR (RT-qPCR). Twenty transcripts (4 differentially expressed [adj. P ≤ 0.1], 9 highly up- or downregulated [LFC ≥ 0.78 or LFC ≤ −0.78], and 7 relatively abundant in both naive and MRFV-exposed leafhoppers but without differential expression) were selected for RT-qPCR validation. Transcripts were selected from among those with putative functions in immune response, protein synthesis and fate, metabolism, energy, transcription, and oxidative response to MRFV exposure. RNA was extracted using the RNeasy minikit as described above. Residual DNA was removed using an on-column DNase treatment, and then DNase was removed in subsequent wash steps according to the manufacturer’s protocol. For each sample, 1 μg of total RNA was reverse transcribed using SuperScript III transcriptase (Thermo Fisher Scientific, Columbus, OH, USA) according to the manufacturer’s protocol. Briefly, 1 μl 50 μM oligo(dT)20, 1 μl 10 mM deoxynucleoside triphosphate (dNTP), and 1 μg of total RNA with water to 13 μl were incubated at 65°C for 5 min and then chilled on ice for 1 min., after which 4 μl 5× first-strand buffer, 1 μl 0.1 M dithiothreitol (DTT), 1 μl RNaseOUT, and 1 μl of SuperScript III RT (200 units/μl) were added and incubated at 50°C for 50 min, after which the reaction was inactivated by heating at 70°C for 15 min. Quantitative PCR was performed using 100 ng of cDNA, 1 μl each of forward and reverse primers (10 μM; Table S8), and 10 μl of SsoAdvanced Universal SYBR green Supermix (Bio-Rad, Hercules, CA, USA) in a CFX96 real-time system (Bio-Rad, Hercules, CA, USA) with the following PCR conditions: 95°C denaturation for 3 min followed by 39 cycles of 95°C for 10 s, 60°C for 30 s, and a final temperature increment of 0.5°C for 5 s from 65°C to 95°C for melting curve analysis. Four biological replicates (samples) and three technical replicates (three RT-qPCR runs) of each sample were analyzed for each treatment. Analysis of relative gene expression for leafhoppers was carried out with Bio-Rad CFX manager 3.1 software (Bio-Rad, Hercules, CA, USA). The ΔΔCt method was employed for the calculation of the relative gene expression level (57, 58). Primer efficiency tests were performed (Table S8) and dissociation curves for each amplicon were analyzed to verify the specificity of each amplification reaction, where efficiency was measured from the slope of the dilution series Ct values such that efficiency E value is equal to 10−1/slope and percent efficiency is equal to (E − 1) × 100%; thus, an E value of 2.00 is 100% efficiency or doubling each cycle (https://www.bio-rad.com/webroot/web/pdf/lsr/literature/Bulletin_5279.pdf). Transcript expression levels were normalized against three reference genes, 60S ribosomal protein L27, tubulin beta chain-like, and 60S ribosomal protein L35a (Table S9). Fold change between MRFV-exposed and naive leafhoppers was calculated as (REMRFV-exposed − REnaive)/REnaive. The correlation between fold change in RT-qPCR and the fold change calculated from RNA-Seq differential expression analyses was statistically analyzed by Spearman rank-order correlation using Prism (Prism Software Corporation, Irvine, CA, USA).

Data availability.

Data supporting the results of this article are included. Sequence data are stored in the NCBI short sequence read (SRA) archive as BioProject PRJNA579843 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA579843). The assembled 56,116 transcripts were deposited at DDBJ/EMBL/GenBank (accession GITV00000000).