The studies were carried out under the glasshouse of ICAR-Indian Institute of Wheat and Barley Research (IIWBR), Karnal, Haryana, India during 2018. The seeds of chosen wheat T. durum genotype, A-9-30-1, and tolerant wheat T. aestivum genotype, HD2967 were obtained from the Germplasm Resources Unit (GRU) of ICAR-IIWBR, Karnal. For the study, each genotype was planted in a plastic pot containing a mixture of farmyard manure, vermicompost, and cocopeat (2:2:1). Corn leaf aphids (R. maidis) colonies were established by collecting aphids from wheat fields and rearing them on wheat susceptible genotype A-9-30-1 in a glass house. The aphids were maintained for several generations before the initiation of the studies to get a pure aphid population. Host resistance response was recorded by measuring nymphiposition, nymphal duration, and nymphal mortality of aphids on the susceptible genotype A-9-30-1 and the tolerant genotype HD2967. The host plant resistance response studies (nymphiposition, nymphal duration, and nymphal mortality) of R. maidis on both the genotypes HD2967 and A-9-30-1 were recorded during January-March, 2018.

To record nymphiposition, five adult aphids of R. maidis were released on each test genotype inside a glass chimney for 24 h to get pre-conditioned nymphs, and then the adult aphids were removed after the laying of young nymphs. These young ones were kept on each test genotype for 10 days until they underwent their last molt and attained a reasonable size. The aphids were transferred to a seven-day-old test plant of each genotype. The developing aphids were regularly observed and the numbers of nymphs deposited/plant/day were counted and removed daily. Nymph production was recorded until the death of the mother aphid. The total number of nymphs laid throughout the life of an aphid was counted, and plants were randomized by genotype with 10 replications. The number of nymphs deposited by each aphid till death was summed up. Mean fecundity and the number of nymphs/aphid/day were calculated.

For measuring the nymphal duration and their survival, five adult aphids were kept on seven-day-old plant of each tested genotype. After 12 h, the adult aphids were removed after nymphiposition on the tested plants, and only 10 neonate nymphs were allowed to grow, and the rest of the nymphs were removed. After 10 days of total feeding period, the number of nymphs survived and dead was recorded. The nymphal duration period was calculated from birth till the last molt or laying of first young one by an aphid adult.

At the two-leaf stage (12-day old plants), 20 apterous adult aphids of R. maidis were confined to the first leaf of wheat seeding of each test genotype (A-9-30-1 & HD2967) inside a glass chimney and were allowed to feed on test genotypes for 48 h. New-born nymphs produced by aphid adults were carefully removed every 12 h using a brush. After 48 h of feeding, all the aphids were removed, and leaf tissues of approximately 2.5 × 2.5 cm² from the aphid feeding sites of each plant were harvested and flash frozen with liquid nitrogen for further processing of RNA extraction. Three leaf sections covering the aphid feeding sites were collected from three independent plants of each chosen genotype and pooled to form one biological replicate. Three biological replicates were performed for each treatment. Similarly, sample collection was also done from uninfested plants of each genotype that served as a control.

Total RNA was isolated from young infested and uninfested leaves of susceptible genotype, A-9-30-1 and tolerant genotype, HD2967 infested through the ZR (control) plant RNA miniprep (ZYMO Research) kit. The quality and quantities of the isolated RNA were checked on a 1% denaturing RNA agarose gel and Nanodrop (NanoDrop technology, USA), respectively. RNA samples having RIN values of more than 8.0 and an OD 260/280 ratio of more than 1.7 were used for RNA-sequencing. The isolated mRNAs were purified using oligo (dT) magnetic beads (Illumina, Inc) from total RNA.RNA-seq pair-end libraries were prepared using the Illumina TruSeq stranded mRNA sample prep kit. These short fragments of mRNA were used as templates for first-strand cDNA conversion using SuperScript II and Act-D mix to facilitate RNA-dependent synthesis for each sample (Hyun et al., 2014; Gutierrez et al., 2017). The first-strand cDNA was then used as a template strand for second-strand synthesis by using a second-strand mix. These cDNA fragments were then purified by using Ampure XP beads followed by A-tailing and adaptor ligation. The products of the ligation reaction were purified on agarose gel electrophoresis that was enriched by PCR for final cDNA library generation (Sharma et al., 2019). These PCR-enriched libraries were further analyzed in the 4200 tape station system of Agilent technologies by using D1000 screen tape (high sensitivity) for quality and quantity checks. The mean of the library distribution was 547 bp for A-9-30-1 (infested leaf), 555 bp for A-9-30-1 (un-infested leaf), 520 bp for HD2967 (infested leaf), and 512 bp for HD2967 (un-infested-leaf), respectively. These cDNA libraries were loaded onto the NextSeq 500 for cluster generation and sequencing. Pair-end sequencing allows the template fragments to be sequenced in both the forward and reverse directions on the NextSeq 500 sequencing platform (Illumina, Inc.) to obtain 2×100 bp paired-end reads.

For obtaining high-quality reads, reads with adapter sequences, ambiguous reads (reads with unknown nucleotides “N” larger than 5%), and low-quality reads (reads with more than 100% quality threshold (QV)<20 phred score) were removed by the Trimmomatic v0.35 software. FastQC was used for quality analysis of cleaned reads from all samples. After this, a reference genome index was set by using Bowtie2 software, whereas the remaining clean reads from paired-end sequencing were mapped to the reference genome (http://www.ensembl.org/index.html) through TopHat2. HtSeq software was used to count the read numbers mapped to every gene. Differential expression of genes was analyzed through the DESwq R package (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html). The significant differences in gene expression were calculated by the absolute value of log₂ fold change>1 and P-value<0.05. After an aphid attack, functional enrichment analyses such as Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) were used to identify the DEGs. The Blast2GO program was used for gene annotation against the Gene Ontology database (http://geneontology.org/) (Ivamoto et al., 2017). KEGG pathways against the KEGG database (http://www.genome.jp/kegg) with a P-value < 0.05 were analyzed, which were significantly deemed in the DEG analysis.

For the validation of sequencing results, total RNA was isolated from young and healthy leaves of A-9-30-1 infested, A-9-30-1 uninfested, HD2967 infested, and HD2967 uninfested by using the Real genomics Hiyield™ total RNA mini kit (Real Biotech Corporation Ltd.). Then, the Thermo Scientific Revert Aid first-strand cDNA synthesis kit was used for the reverse transcription. Pathogenesis-related gene-specific primers were developed by using Primer3 software (Rozen and Skaletsky, 2000). The reactions of real-time PCR were carried out in triplicate using the ABI PRISM 7700 Real-time PCR system with the ABI 7700 Fast sequence detector (Applied Biosystems, USA). The expression of variant selected genes was normalized by the housekeeping gene ‘Actin’ of wheat. The designed primers ( Supplementary information: Table 1 ) were tested for single-band amplification through conventional endpoint PCR. Melting curve analysis was performed using qPCR in a reaction containing 1X SYBR Green Master mix (Applied Biosystems, USA), 1 µl cDNA template, 10 pmol of each primer, and a final volume of 20 µl was prepared by adjusting with nuclease-free water (Sharma et al., 2019). The conditions that were followed for the real-time PCR experiment included initial denaturation at 95˚C for 30 s, followed by 35 cycles of denaturation at 95˚C for 10 s, and extension for 10 s based on primer Tm followed by a thermal dissociation curve. The relative expression level was analyzed using the 2^-ρρct method (Schmittgen and Livak, 2008).

The data for aphid resistance identification among genotypes were analyzed by one-way analysis of variance (ANOVA, P < 0.05). The results of biological parameters i.e. nymphiposition, and nymphal duration and survival are reported as mean values with standard errors (SE). Fisher’s least significant difference (LSD) test was used to separate means (P < 0.01). The software SPSS v17.0 (IBM Inc., Chicago) was used for data analyses.

There was a significant difference in the nymphiposition of R. maidis recorded on wheat susceptible genotype A-9-30-1 and tolerant genotype HD2967. The nymphiposition (fecundity) of R. maidis was higher (40.72 ± 0.54 nymphs) on A-9-30-1 as compared to HD2967 (10.81 ± 0.75 nymphs). There were also significant differences in nymphal duration and survival between the two genotypes studied. The nymphal duration of 12.75 days was recorded on A-9-30-1 as compared to 20.89 days on genotype HD2967. The nymphal mortality rate on genotype HD2967 was higher (45.79%) than on A-9-30-1 (29.36%) ( Figure 1 ).

By using Illumina sequencing technology, a total of 59,691,508 (8.99 Gb), 56,052,968 (8.43 Gb), 47,363,964 (7.13 Gb), and 56,403,641 (8.49 Gb) numbers of reads were generated for A-9-30-1 infested, A-9-30-1 uninfested, HD2967 infested, and HD2967 uninfested, respectively. These high-quality reads were mapped onto the reference genome through TopHat v2.1.1 with default parameters. The mapping statistics for A-9-30-1 infested leaf and A-9-30-1 uninfested leaf were 47% and 87.8%, respectively. However, it was 66.8% and 89.3% for HD2967 infested leaf and HD2967 uninfested leaf, respectively ( Table 1 ).

Transcriptome analysis of two selected genotypes, A-9-30-1 and HD2967 revealed a significant difference in gene expression after aphid attack. In combination 1 (tolerant genotype HD2967 infested vs uninfested), 1183 and 959 genes were significantly upregulated and downregulated, respectively. Similarly in combination 2 (susceptible A-9-30-1 infested vs uninfested), a reduction in the number of upregulated (1055) and downregulated genes (765) was observed. On the other side, comparative analysis of susceptible and tolerant genotypes without aphid infestation (combination 3), it was noticed that 353 and 1175 genes were upregulated and downregulated, respectively. Further in combination 4 when both genotypes were infested with aphids, only 212 genes were significantly upregulated and 1009 genes were significantly downregulated. ( Table 2 ). The proprietary R script was used to depict the graphical representation and distribution of differentially expressed genes in control as well as treated samples. The ‘volcano plot’ arranges expressed genes along dimensions of biological as well as statistical significance ( Figure 2 ).

The DEGs obtained were subjected to GO enrichment analysis for the identification of the main biological functions of DEGs in wheat under biotic stress conditions (Young et al., 2010). All DEGs were classified into three categories, i.e. cellular components (CC), molecular function (MF), and biological processes (BP) (Ashburner et al., 2000; Sharma et al., 2019) Figures 3A–C ). Annotation of DEGs of HD2967 (infested vs uninfested) revealed that a total of 4,364 GO terms corresponding to 2062 molecular functions, 1323 biological processes, and 979 cellular components were assigned. “Transferase activity,” “catalytic activity (acting on a protein),” “organic cyclic compound binding,” and “heterocyclic compound binding activity” were the most dominant GO terms in the MF category, while “nitrogen compound metabolic process,” “regulation of the metabolic process,” “oxidation-reduction process,” “cell communication,” and “cellular component organization” were the most dominant GO terms in the BP category. However, in the case of CC, most of the DEGs were classified into “cell part”, “intracellular part”, “cell”, and “protein-containing complex”, followed by “membrane-bounded organelle”. However, the annotation of DEGs of A-9-30-1 (infested vs uninfested) showed a total of 3,255 GO terms, which contained 1531 molecular functions, 1012 biological processes, and 712 cellular components. The dominant terms in the MF category were “catalytic activity, acting on RNA”, “lyase activity”, “carbohydrate binding” and transporter activity. In the case of BP, “cellular metabolic processes”, “catabolic processes”, “biosynthetic process” and “small molecule metabolic process”, whereas, “cell part”, “intracellular part”, “endomembrane system”, and “protein-containing complex” were the most dominant GO terms ( Figures 3A–C ).

On a comparison of infested HD2967 and infested A-9-30-1, a total of 621 GO terms was identified which corresponds to 257 MF, 193 BP, and 171 CC. The GO assignment for MF involves “binding”, “organic cyclic compound binding”, “heterocyclic compound binding” and “small molecule binding” whereas, BP involved “biosynthetic process”, “cellular metabolic process”, “nitrogen compound metabolic process” and oxidation-reduction process. In the case of CC, the dominant GO terms include “cell part”, “intracellular” and “transferase activity” ( Figures 3A–C ).

The KEGG pathway analysis was performed for the understanding of active biological pathways during aphid attack in wheat ( Supplementary information: Excel-Sheets showing all pathways). The analysis revealed that out of 127 KEGG pathways, 25 were significantly active enriched at a p-value <0.05, covering five main KEGG categories such as metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems. The KEGG pathways associated with DEGs under aphid attack conditions were widely associated with metabolism, biosynthesis of secondary metabolites, and plant hormone signal transduction. The genes involved in these 25 pathways were upregulated in HD2967 as well as in A-9-30-1 after aphid attack ( Table 3 ).

R. maidis attack on HD2967 and A-9-30-1 wheat showed its negative effect on the photosynthesis process in which a large number of genes associated with light-harvesting and photosystem, i.e. chlorophyll a-b binding proteins, photosystem I and II, and ferrochelatase were significantly downregulated ( Table 4 ). The expression levels of transcripts of carbonic anhydrase and ribulose bisphosphate carboxylase oxygenase (RuBisCO) genes were also significantly downregulated, as these genes are mainly involved in the Calvin cycle. Transcription profiles for sucrose and starch metabolism were also analyzed for a few genes, which showed that one trehalose-6-phosphate synthase gene, one sucrose synthase-3 gene, and six beta-glucosidase genes were also significantly downregulated during the attack of the R. maidis. However, the sucrose-phosphate gene is upregulated during the attack of R. maidis on HD2967 wheat leaves. The level of nitrate reductase transcripts was strongly downregulated by R. maidis in both, i.e., in HD2967 and A-9-30-1, whereas glutamate dehydrogenase transcript levels were significantly upregulated in HD2967 wheat leaves infested with R. maidis ( Table 4 ).

Characterization of genes involved in jasmonic acid (JA), salicylic acid (SA), and ethylene (ET) defense pathways were used to estimate plant defense response to R. maidis attack ( Table 5 ). It was found that in SA biosynthesis pathway, 15 phenylalanine ammonia-lyase (PAL) genes and 9 pathogenesis-related genes (PR) were significantly upregulated in response to R. maidis feeding in both, i.e., in T. durum (A-9-30-1) and T. aestivum (HD2967). When the expression of these PAL genes and PR genes was compared, it was observed that it was significantly lower in T. durum (A-9-30-1) than their expression in T. aestivum (HD2967). Furthermore, when R. maidis attacked T. durum (A-9-30-1) and T. aestivum (HD2967), the expression level of two AOS genes varied from 1.08 to 2.07-fold in A-9-30-1 and 4.14 to 4.95-fold in HD2967. The expression of three LOX expressions varied from 3.13 to 5.83-fold in A-9-30-1 and 4.10 to 6.02-fold in HD2967, indicating its significant up-regulation in JA biosynthesis. The expression level of two PI genes (JA-defense responsive genes) was found to be higher in the tolerant genotype, T. aestivum (HD2967), as compared to the susceptible genotype, T. durum (A-9-30-1). It varied from 3.12 to 4.69-fold in A-9-30-1 and 3.14 to 5.71-fold in HD2967. Furthermore, it was found that in the ET signaling pathway, one 1-aminocyclopropane-1-carboxylate synthase (ACS) gene, four 1-aminocyclopropane-1-carboxylate oxidase 1 (ACO), three genes belonging to ET-responsive transcription factors and ET-insensitive protein genes were all significantly upregulated in response to R. maidis feeding in A-9-30-1 and HD2967 ( Table 5 ).

A large number of pathogenesis/defense-related genes were found to be differentially expressed under aphid-attack conditions. In total, 21 important pathogenesis/defense-related genes were selected for expression analysis in control as well as infested samples of wheat, i.e., T. durum and T. aestivum against R. maidis ( Figure 4 ). It was found that a serine/threonine kinase DEG showed upregulated expression in infested HD2967 whereas downregulated expression was in uninfested A-9-30-1. Similarly, a gene encoding signaling-associated proteins and other stress tolerance were observed, namely “Q5ZD81 Calmodulin-like protein”. It showed increased expression in infested HD2967 and decreased expression in uninfested A-9-30-1. The DEGs related to plant cell protection from oxidative damage by ROS scavenging, such as “Q8S702 Glutathione S-transferase” and “Q43212 Peroxidase precursor” showed their increased expression in infested HD2967, whereas downregulated expression was seen in uninfested A-9-30-1. Overall, Q9XEN7 β-1, 3-glucanase had the highest expression in the infested T. aestivum and T. durum whereas P12940 Bowman-Birk trypsin inhibitor had the lowest expression in the infested and uninfested T. aestivum and T. durum ( Figure 4 ).