Identification of Candidate Genes for English Grain Aphid Resistance from QTLs Using a RIL Population in Wheat

: The English grain aphid (EGA) ( Sitobion avenae F.) is one of the most destructive species of aphids in wheat-( Triticum aestivum L.) planting areas worldwide. Large quantities of insecticides are usually used to control aphid damage. The identification of new EGA-resistant genes is necessary for sustainable wheat production. The objective of this study was to identify candidate genes for EGA resistance from stable quantitative trait loci (QTLs). We previously constructed a genetic map of unigenes (UG-Map) with 31,445 polymorphic sub-unigenes via the RNA sequencing of ‘TN18 × LM6’ recombinant inbred lines (TL-RILs). The relative aphid index (RAI) for the TL-RILs was investigated for two growing seasons, with three measured times (MTs) in each season. Using the UG-Map, 43 candidate genes were identified from 22 stable QTLs, with an average of 1.95 candidate genes per QTL. Among the 34 candidate genes annotated in the reference genome Chinese Spring (CS) RefSeq v1.1, the homologous genes of seven candidate genes, TraesCS1A02G-319900


Introduction
Aphids are among the most devastating pests to crop production because of their short life cycles and enormous reproductive potential [1,2].Common wheat (Triticum aestivum L.) supplies approximately one-fifth of the calories and protein for human diets [3].Aphids are the most serious insect pests for wheat.The aphids that severely affect wheat production around the world are the English grain aphid (EGA) (Sitobion avenae F.), Russian wheat aphid (RWA) (Diuraphis noxia M.), greenbug (Schizaphis graminam R.), bird cherry-oat aphid (Rhopalosiphum padi L.), etc. [4][5][6].Among them, the EGA is one of the most destructive of the aphid species in wheat-planting areas worldwide [7,8].The EGA mainly feeds on wheat at the jointing and booting stages and causes considerable damage to wheat yield and quality by phloem sap sucking, excreting honeydew and transmitting plant viruses [9][10][11].
This results in annual wheat yield losses ranging from 15 to 30% and as high as 60% in severe cases [6,12].
Over the long course of their evolution, plants have evolved a variety of adaptations to reduce the damage caused by phytophagous insects.Their mechanisms of resistance against aphids have traditionally been divided into three classic categories: antibiosis, tolerance and antixenosis [13,14].A resistant cultivar may exhibit a combination of these categories of resistance [15].Antibiosis is the main resistance mechanism of aphid resistance in plants.Plants produce chemicals that are harmful to aphids or lack the nutrients that are necessary to aphids, resulting in increased aphid mortality and reduced reproductive rate [15,16].Antibiosis is associated with the nutritional composition and secondary chemicals of wheat such as soluble sugars, amino acids, DIMBOA, alkaloids, and phenolics [4,17].Tolerance is a characteristic that allows the host plant to grow and develop normally even after aphid infestation, showing a certain regeneration or compensation ability at the individual or population level and not greatly reducing production.This is a passive resistance mechanism that acts through repair or de novo synthesis of photosystem proteins [18].Antixenosis, or what is referred to as non-preference, involves morphological resistance in which the aphid is either repelled by or not attracted to the host plant, for feeding as well as for egg-laying.For example, the wax powder, glume color and ear density of wheat affect the feeding and reproduction of aphids [19,20].
At present, large quantities of insecticides are usually used to control aphid damage.However, these substances cause serious environmental pollution, threaten human health and destroy non-targeted beneficial insects [21,22].The use of resistant wheat cultivars can provide an efficient, economic and environmentally friendly approach to reducing the yield losses caused by aphid injury and can produce healthy, pesticide-free products [4,23,24].Therefore, improving the resistance of host-plants has become a high-priority target in wheat breeding programs.To date, a great number of wheat germplasms have been evaluated for EGA resistance using various methods, and many wheat cultivars have been identified as being resistant to EGA [4,[25][26][27][28].Identifying EGA-resistant genes and linked molecular markers is crucial for gene cloning and enhances the efficiency of molecular breeding in wheat.This approach enables precise selection and incorporation of EGA resistance into new wheat varieties.Furthermore, EGA resistance in wheat exhibited a typical quantitative trait inheritance pattern [6].Quantitative trait locus (QTL) analysis has proved an effective approach to dissect a complicated quantitative trait into component loci and is an effective strategy for crop improvement [29,30].
The wheat EGA resistant genes and their closely linked molecular markers have been reported.Fu et al. [31] assessed the natural infection of EGA in an F 2 population and identified the dnY gene for EGA resistance on chromosome 7D using SSR markers.Hu [32] artificially inoculated a RIL population in the field and found the Ra gene for EGA resistance on chromosome 7D.Wang et al. [24,33] reported the EGA-resistant Sa1 and Sa2 genes on chromosome 7D by counting the number of EGA in tillers of F 1 , F 2 , BC 1 , F 2:3 and F 3:4 populations under natural infection.Liu et al. [4] found EGA resistance gene RA-1 on chromosome 6AL of Triticum durum through artificial inoculation in field and greenhouse settings using F 1 , F 2 and F 2:3 populations.Using 70 bread wheat accessions, Li and Peng [34] identified four SSR loci significantly associated with EGA resistance, four with EGA tolerance and four SSR loci were detected continuously at different growing stages of wheat.Luo et al. [26] found five potential expressed sequence tags likely associated with EGA resistance and one was located on 7DL using two backcross progenies.Some QTLs and/or candidate genes associated with EGA resistance have also been identified.Zeng et al. [35] constructed a genetic linkage map with a total of 14,694 polymorphic SNPs of a RIL population.One major QTL, QSa3.haust-3D, was identified for EGA resistance on 3DL between 550.25 and 552.98 Mb in Chinese Spring (CS) RefSeq v1.1 [36], and eight putative candidate genes were identified from this region [36].Zhao et al. [6] identified a new locus of qSa-3A for EGA resistance in a narrow region of 34.52-37.50Mb on chromosome 3AS.They identified nine important loci and nine QTLs using genome-wide association study (GWAS) and QTL analysis, respectively, providing the physical intervals in the reference genome CS RefSeq v1.1 [36].A number of candidate genes were identified by expression analysis and gene-silencing technology.Fu et al. [10] reported that the expression of harpin protein Hpa1 inhibits EGAs from feeding from the phloem.This defensive mechanism was shown as enhanced expression of wheat genes encoding phloem lectin proteins (PP2-A1 and PP2-A2) and β-1,3-glucan synthase-like enzymes (GSL2, GSL10 and GSL12).Luo et al. [37] found that the expression levels of the candidate genes SaEST1 and SaEST2, which are the putative photosystem I assembly protein Ycf3 and vegetative cell wall gp1-like protein, respectively, up-regulated the EGA resistance in flag leaves.Wang et al. [38] identified lipoxygenase linoleate (LOX) genes, TaLOX5, 7, 10, 24, 29, 33, specifically expressed post-EGA infestation using RNA-seq data and qRT-PCR analysis.Zhai et al. [39] elucidated that the TaMYB19, TaMYB29 and TaMYB44 are co-regulators of wheat plant phloem-based defence (PBD) using qRT-PCR and the gene silencing method.Furthermore, a 13-lipoxygenase gene, designated as TtLOX, was cloned and functionally characterized.TtLOX is involved in wheat resistance against EGA and the regulation of jasmonate biosynthesis [17].
We previously constructed a genetic map of unigenes (UG-Map) via the RNA sequencing of 'TN18 × LM6' recombinant inbred lines (RILs) [40].The objective of this study was to perform QTL analysis for EGA resistance using this RIL population, and then to identify candidate genes from stable QTLs.

Plant Materials
A set of 184 RILs derived from the cross of 'TN18 × LM6' using a single-seed descent approach were used to perform QTL analysis [41,42].For convenience, we named this RIL population as TL-RIL.The female parent TN18 is a variety developed by our group.TN18 has become a core parent, with more than 60 authorized cultivars developed as of 2023.The male parent LM6 is an elite wheat line developed by the Linyi Academy of Agricultural Sciences.The two parents show differences in aphid resistance, with LM6 having higher resistance than TN18 overall.

Experimental Design and Trait Measurement
The trials were carried out at the Experiment Station of the Shandong Agricultural University (Tai'an, China).The test field was covered at the top and open around the sides.The rainfall during wheat growth was approximately 160 mm, and the soil structure was loamy soil.The compound fertilizer (N-15%, P-15% and K-15%) was applied at 750 kg per ha, watered two times at the jointing stage and flowering stage (600 m 3 per ha for one time) under full soil moisture sowing conditions and no pesticides were used.Seeds were sown on 10-11 October, and plants were harvested on 9-10 June the next year.

Data Analysis, QTL Detection and Candidate Gene Identification
Analysis of variance (ANOVA) was performed using SPSS 17.0 software (SPSS Inc., Chicago, IL, USA).Genotypes and environments were considered as two factors using the data from 184 lines under six MTs.All factors involved were considered sources of random effects.
The UG-Map for the TL-RILs was used for QTL analysis.Each line of the TL-RILs was sequenced using RNA-Seq technology.The clean reads were mapped to Chinese Spring (CS) RefSeq v1.1 [36].The polymorphic unigenes were used to construct the UG-Map based on their physical positions.The final UG-Map included 31,445 polymorphic sub-unigenes, which represented 27,983 unigenes and 118,120 SNPs/InDels [40].The mRNA expression levels for a SNP/InDel of candidate genes were determined by the number of reads mapped on the reference genome for all lines of TN18 or LM6 genotypes in the TL-RILs.When a gene included more than one SNP/InDel, we regarded the expression levels of the gene as different and at the least one SNP/InDel as different.
Windows QTL Cartographer 2.5 software (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm, accessed on 5 January 2019) was used to perform QTL mapping, and compositeinterval mapping (CIM) was selected to search for QTLs.The parameter setup was as follows: "Model 6 standard analysis" with a walk speed of 0.5 cM, "forward and backward" regression for the selection of the markers to control for the genetic background, up to five control markers and a blocked window with a size of 10 cM to exclude closely linked control markers at the testing site [42].The threshold for declaring the presence of a significant QTL in each MT was LOD ≥ 3.0, and the QTL interval was determined by LOD ≥ 2.0.A QTL was defined as a stable QTL when it was detected over (AV+2MTs).
The unigenes covered by or adjacent to the interval of QTLs were regarded as candidate genes for the corresponding QTLs.

Variant Types and Excellent Mutations for the Candidate Genes
For the 34 candidate genes annotated in CS RefSeq v1.1 [36], 20, 1 and 2 were nonsynonymous, UTR and synonymous mutation(s) in exons, respectively, and 11 were mutations in introns (Table S2).Among all the 43 candidate genes, the reads of 29 genes were significantly different between the TN18 and LM6 genotypes in the TL-RILs, indicating that their mRNA expression levels were changed.The mRNA expression levels of the other 14 candidate genes between the TN18 and LM6 genotypes were not different (Table S2).According to the DNA sequences of the parents and a part of RILs, the promoter region (−2000 bp from the start site of the 5 ′ UTR) of 16 candidate genes in CS RefSeq v1.1 [36] and three candidate genes in TL-RILs were mutated (Table S3).For the two genes with synonymous mutation (TraesCS4A02G287900 and TraesCS6A02G000600), mRNA expression levels were all changed, and the promoter region of TraesCS4A02G287900 was mutated.For the 11 genes with mutations in introns, six (TraesCS1B02G565900LC, TraesCS2D02G565600LC, TraesCS5B02G378700, TraesCS5D02G458800, TraesCS6B02G-644000LC, TraesCS7A02G123300) showed significantly different mRNA expression levels between the TN18 and LM6 genotypes in the TL-RILs, and one (TraesCS7A02G123300) was mutated in the promoter region (Tables S2 and 3).
For all 43 candidate genes, the excellent mutations of 33 candidate genes with 82 SNPs/InDels were TN18 type, for which the additive effects of corresponding QTLs were negative, with the parent TN18 decreasing RAI values.Whereas, the excellent mutations of ten candidate genes with 51 SNPs/InDels were LM6 type, for which the QTL effects were positive, with LM6 decreasing RAI values (Table S3).

Candidate Genes in This Study
In this study, we identified 43 candidate genes with resistance to EGAs on 14 chromosomes.Four resistance genes for EGA, dnY [31], Ra [32], Sa1 and Sa2 [24,33], were previously located on the chromosome 7D.Our investigation did not reveal any candidate genes on this chromosome.The work of Liu et al. [4], which identified RA-1 closely linked to the SSR markers Xwmc179, Xwmc553 and Xwmc201 on 6AL, falls outside our findings of four candidate genes on chromosome 6A (Tables 3 and S1).Similarly, the five SSR markers (Xgwm11, Xgwm294, Xgwm526b, Xgwm192b and Xgwm613b) highlighted by Li and Peng [34] were on the same chromosome as this study, but the physical positions are different.
In yeast, both synonymous and non-synonymous mutations disturbed the levels of mRNA expressions and affected the fitness effects, which negated the presumption that synonymous mutations are neutral or nearly neutral [43].In this study, we also found that the RNA expression levels were changed for two candidate genes with synonymous mutations and six genes with mutations in introns.Contrastingly, the RNA expression levels were not changed for five genes with intronic mutations, indicating that they might be less reliable candidates.

Seven Candidate Genes Have Been Reported to Play Roles in Aphid Resistance
We used the UG-Map of the TL-RILs that included 31,445 sub-unigenes to identify candidate genes for EGA resistance from QTLs.The number of candidate genes per QTL was 1.95, which can directly identify the candidate genes for EGA resistance from QTLs.Among the 34 candidate genes annotated in CS RefSeq v1.1 [36], the homologous genes of seven candidate genes (7/34 = 20.6%)have been reported to play roles in aphid resistance (Table S4).We infer that these seven genes are strongly associated with EGA resistance in wheat.This result demonstrated that the candidate genes found in this study were highly reliable.As far as we know, only one resistant EGA gene, TtLOX, has been isolated in wheat (Triticum turgidum) [17]; these highly reliable candidate genes for EGA resistance provided options for further confirmation.The candidate genes in this study need to be validated using transgenic methods such as CRISPR/Cas9 system.
The seven candidate genes, which are strongly associated with EGA resistance, were TraesCS1A02G319900, TraesCS1B02G397300, TraesCS2D02G460800, TraesCS4A02G-015600LC, TraesCS5B02G329200, TraesCS6A02G000600 and TraesCS6A02G418600LC (Table S4).The details of their resistance to aphids and the variation in TL-RILs of the seven genes were as follows.
TraesCS1A02G319900 for QRai-1A-8731 and TraesCS5B02G329200 for QRai-5B-10270 are annotated as defensins.Plant defensins have a negative effect on phloem-feeding aphids [44].Overexpression of biologically safe Rorippa indica defensin has been found to enhance aphid tolerance in mustard, with the mean number of surviving mustard aphids on the transgenic and control plants were found to be statistically significant [45].TraesCS1A02G319900 had three non-synonymous substitutions (Ile37Val, Gly62Asp and Trp69Arg), one synonymous variant in the exon and 16 SNPs and two InDels in the promoter (Tables S2 and S3, Figures S2A, S3A and S4A).TraesCS5B02G329200 had one non-synonymous substitution (Gly58Ser) and one synonymous variant in the exon, one SNP in the 3 ′ UTR, two variants in the introns and 10 SNPs in the promoter (Tables S2 and S3, Figures S2B, S3B and S4B).For both of the two genes, the expression levels were significantly different between the TN18 and LM6 genotypes of the TL-RILs (Table S2).
TraesCS1B02G397300 for QRai-1B-12814 is annotated as a cyclopropane-fatty-acylphospholipid synthase (CFAS).In cotton, CFAS was significantly upregulated in the Aphis gossypii-damaged plants when compared with undamaged plants [46].TraesCS1B02G-397300 had one non-synonymous substitution (His80Arg) in the exon, and eight SNPs in the promoter.The expression levels were significantly different between the TN18 and LM6 genotypes of the TL-RILs (Tables S2 and S3, Figures S2C, S3C and S4C).
TraesCS2D02G460800 for QRai-2D-6068 is annotated as a hydroxyproline-rich glycoprotein family protein, putative (HRGP).HRGPs are important structural components of plant cell walls and accumulate in response to infection as an apparent defense mechanism [47,48].Aphid feeding on the sugarcane-aphid-resistant sorghum line resulted in the increased expression of genes related to cell wall formation, including the HRGP gene [49].TraesCS2D02G460800 had one non-frameshift deletion (nine amino acid residue deletions) in the exon and no SNP/InDel in the promoter.The expression levels were significantly different between the TN18 and LM6 genotypes of the TL-RILs (Table S2, Figures S2D and S3D).
TraesCS4A02G015600LC for QRai-4A-371 is annotated as an auxin-responsive protein (ARP).In sorghum, the ARP gene involved in processes that halt accumulation of free IAA showed 6.6-and 3.9-fold increases in expression in response to aphid herbivory in resistant and susceptible genotypes, respectively.This suggests that up-regulation of this gene may have triggered a change in the allocation of resources from growth to defense and conferred greater resistance to the aphids [50].TraesCS4A02G015600LC had three SNPs in the intron and no SNP/InDel in the promoter.The expression levels were not significantly different between the TN18 and LM6 genotypes of the TL-RILs (Table S2 and Figure S2E).
TraesCS6A02G000600 for QRai-6A-80 is annotated as a protein disulfide-isomerase (PDI).PDIs play a role in the redox regulation of target enzymes and transcription factors [51].The expression of PDI was significantly different between aphid-resistant and aphid-susceptible wheat plants under aphid attack.Additionally, PDIs may be involved in wheat tolerance to RWA by participating in plant cell protection against oxidative stress [52].TraesCS6A02G000600 had two synonymous variants in the exon and no SNP/InDel in the promoter.The expression levels were significantly different between the TN18 and LM6 genotypes of the TL-RILs (Table S2 and Figure S2F).
TraesCS6A02G418600LC for QRai-6A-7415 is annotated as an F-box family protein.In cucumber, the expression level of the gene encoding the F-box family protein was up-regulated after aphid infestation.This indicated that aphid infestation enhanced the expression of the gene, suggesting that they could play an important role in aphid resistance in cucumber [53].F-box domain-containing proteins play a crucial role during aphid herbivory by fine tuning the ethylene pathway in melon [54].TraesCS6A02G418600LC had two non-synonymous substitutions (Trp24Arg and Ser38Ile) in the exon, and six SNPs in the promoter.The expression levels were not significantly different between the TN18 and LM6 genotypes (Tables S2 and S3, Figures S2G, S3E and S4D).

Conclusions
We used the UG-Map of the TL-RILs to identify candidate genes for EGA resistance from QTLs at six MTs over two growing seasons.A total of 43 candidate genes for 22 stable QTLs were identified.Most of the candidate genes were newly identified.Interestingly, among the thirty-four candidate genes annotated in the CS RefSeq v1.1 [36], the homologous genes of seven candidate genes have been reported to play roles in aphid resistance.These genes were strongly associated with EGA resistance in wheat, which demonstrated that the candidate genes found in this study were highly reliable.The candidate genes in this study should facilitate the cloning of EGA-resistant genes and improve resistance to English grain aphids in wheat breeding.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy14030637/s1, Figure S1: Distribution of relative aphid index (RAI) under different measured times (MTs) and average value (AV).Figure S2: Difference of DNA sequences between TN18 and LM6 for candidate genes have been reported to play roles in aphid resistance.Figure S3: Difference of amino acid sequences between TN18 and LM6 for candidate genes that have been reported to play roles in aphid resistance.Figure S4: Difference of promoter region (−2000 bp from the start site of the 5 ′ UTR) between TN18 and LM6 for candidate genes that have been reported to play roles in aphid resistance.Table S1: Stable QTLs and their candidate genes for RAI in the TL-RILs.Table S2: Mutant types and expression levels of reads for the candidate genes.Table S3: SNPs/InDels in the promoter regions of the candidate genes.Table S4: Homologues of candidate genes previously reported to play roles in aphid resistance.

Figure 1 .
Figure 1.Locations of the 22 stable QTLs and their 43 candidate genes for relative aphid index (RAI) using the 'TN18 × LM6′ RILs.The red letters are QTLs and the black letters are corresponding candidate genes.CS and TL are abbreviations of gene names of "TraesCSxx02G" and "TraesTLxx02G".

Figure 1 .
Figure 1.Locations of the 22 stable QTLs and their 43 candidate genes for relative aphid index (RAI) using the 'TN18 × LM6' RILs.The red letters are QTLs and the black letters are corresponding candidate genes.CS and TL are abbreviations of gene names of "TraesCSxx02G" and "TraesTLxx02G".

Table 1 .
Phenotypic parameters for relative aphid index (RAI) of the TL-RILs and their parents.

Table 2 .
Analysis of variance for RAI.