QTL Mapping by SLAF-seq and Expression Analysis of Candidate Genes for Aphid Resistance in Cucumber

Cucumber, a very important vegetable crop worldwide, is easily damaged by pests. Aphid is one of the most serious cucumber pests and frequently cause severe damage to commercially produced crops. Understanding the genetic mechanisms underlying pest resistance is important for aphid-resistant cucumber varieties breeding. In this study, two parental cucumber lines, JY30 (aphid susceptible) and EP6392 (aphid resistant), and pools of resistant and susceptible (n = 50 each) plants from 1000 F2 individuals derived from crossing JY30 with EP6392, were used to detect genomic regions associated with aphid resistance in cucumbers. The analysis was performed using specific length amplified fragment sequencing (SLAF-seq), bulked segregant analysis (BSA), and single nucleotide polymorphism index (SNP-index) methods. A main effect QTL (quantitative trait locus) of 0.31 Mb on Chr5, including 43 genes, was identified by association analysis. Sixteen of the 43 genes were identified as potentially associated with aphid resistance through gene annotation analysis. The effect of aphid infestation on the expression of these candidate genes screened by SLAF-seq was investigated in EP6392 plants by qRT-PCR. The results indicated that seven genes including encoding transcription factor MYB59-like (Csa5M641610.1), auxin transport protein BIG-like (Csa5M642140.1), F-box/kelch-repeat protein At5g15710-like (Csa5M642160.1), transcription factor HBP-1a-like (Csa5M642710.1), beta-glucan-binding protein (Csa5M643380.1), endo-1,3(4)-beta-glucanase 1-like (Csa5M643880.1), and proline-rich receptor-like protein kinase PERK10-like (Csa5M643900.1), out of the 16 genes were down regulated after aphid infestation, whereas 5 genes including encoding probable leucine-rich repeat (LRR) receptor-like serine/threonine-protein kinase At5g15730-like (Csa5M642150.1), Stress-induced protein KIN2 (Csa5M643240.1 and Csa5M643260.1), F-box family protein (Csa5M643280.1), F-box/kelch-repeat protein (Csa5M643290.1), were up-regulated after aphid infestation. The gene Csa5M642150.1, encoding probable LRR receptor-like serine/threonine-protein kinase At5g15730-like, was most likely a key candidate gene in cucumber plants in response to infestation. This study provides a certain theoretical basis of molecular biology for genetic improvement of cucumber aphid resistance and aphid resistant variety breeding.


INTRODUCTION
Cucumber, Cucumis sativus L. (2n = 2x = 14), is an agriculturally and economically important vegetable crop worldwide. The aphid, Aphis gossypii Glover, is found in the majority of cucumber producing areas in China and is a serious pest, often causing severe yield loss and reduced quality in cucumber production. Understanding the genetic mechanisms underlying aphid resistance is important for breeding aphid resistance and improving cucumber quality.
Determination of the genetic basis of crop resistance to insect pests has been widely studied. Investigation of sorghum resistance to the greenbug using different resistant plant lines and various aphid biotypes, revealed 3 to 9 associated genomic regions (Agrama et al., 2002;Wu and Huang, 2008). Quantitative trait loci (QTLs) identified as associated with aphid resistance include: five in barley, conferring resistance to the Russian wheat aphid (Mittal et al., 2008); two in cowpea, controlling soybean aphid resistance (Zhang et al., 2009); four additive and two couples of epistatic QTLs in melon, responsive to the meloncotton aphid infestation (Boissot et al., 2010); and one major candidate QTL on chromosome 7 controlling foxglove aphid resistance in soybean (Lee et al., 2015). In addition, a putative QTL in apple for resistance to the rosy apple aphid and another for resistance to the green apple aphid have been localized (Stoeckli et al., 2008), and in peach, two QTLs for resistance to the green peach aphid were identified (Sauge et al., 2004). However, the genetic mechanisms underlying aphid resistance in cucumber remain unclear.
Quantitative trait loci mapping is the main approach for genetic dissection of quantitative traits. It is conventionally conducted by genotyping a large number of individuals in segregating populations, which is labor-intensive, timeconsuming and sometimes costly (Salvi and Tuberosa, 2005). Bulked-segregant analysis (BSA) provides a simple and effective alternative technology to identify molecular markers linked to target genes or QTLs affecting a trait of interest, by genotyping only one pair of pooled DNA samples from two groups of individuals with distinct or opposite extreme phenotypes (Michelmore et al., 1991). Specific length amplified fragment sequencing (SLAF-seq) combining with high-throughput and reduced representation library sequencing is considered an efficient and high-resolution strategy for large-scale genotyping (Sun et al., 2013). In maize, three candidate regions related to genetic and molecular control of meristems on Chr3, containing 51 candidate genes and consisting of 3.947 Mb, were obtained by SLAF-seq (Xia et al., 2015). In cucumber, a major QTL controlling fruit flesh thickness was identified by SLAF-seq and confirmed by simple sequence repeat marker-based classical QTL mapping in 138 F 2 individuals .
In this study, pools of resistant and susceptible bulk samples (n = 50 each) were constructed from 1000 F 2 plants derived from a cross of JY30 (susceptible female parent, P 1 ) and EP6392 (resistant male parent, P 2 ) plants, and were used to detect regions in the cucumber genome harboring major aphid resistance QTLs by SLAF-seq, BSA and single nucleotide polymorphism index (SNP-index) methods. The effect of infestations with different numbers of aphids on the expression of genes identified by SLAF-seq were studied by qRT-PCR.

Aphid Culture
One aphid (Aphis gossypii Glover) was collected from the experimental fields of cucumber at Yangzhou University in the autumn of 2012, and reared on the susceptible cucumber line 'XiaFeng' at 25 • C (18 h)/18 • C (6 h) day/night in relative humidity of 50-60%. The offspring of this aphid were used to infest cucumber plants.

Plant Materials and Aphid Infestation
All plants were cultivated in potting substrate (nutrient availability: total N, P, and K nutrients, 40-60 g/kg; total humus content, ≥350 g/kg; pH 6.5-7.5) in environmental growth chambers maintained at 25 • C (18 h)/18 • C (6 h) day/night, at a light intensity of 12,000 lux (18 h)/0 lux (6 h), and relative humidity between 50 and 60%. An F 2 population of 1000 individuals derived from a cross of JY30 (susceptible female parent, P 1 ) and EP6392 (resistant male parent, P 2 ) was used for genetic analysis and molecular mapping of aphid resistance QTLs. Both EP6392 and JY30 are European cucumber. The commercial fruit length of EP6392 is about 25−30 cm, while 30−35 cm of JY30. The fruit peel of EP6392 is dark green, less thorn, however, the fruit peel of JY30 is green, spiny. The leave of EP6392is dark green, while the leaves of JY30 is green. Seeds of P 1 , P 2 , and F 2 generations were sown on 28 May, 2014. Ten days after sowing, five apterous adult aphids were transferred to the back of the first true leaf per seedling. The number of aphids on individual plants were counted on day 8 after infestation and recorded as aphid scores ranging from 1 to 5, where 1 ≤ 100 aphids, 2 = 101 to 200 aphids, 3 = 201 to 300 aphids, 4 = 301 to 400 aphids, and 5 ≥ 401 aphids per plant (scoring system was modified based on Jun et al., 2012).
The expression of genes identified by SLAF-seq after infestation of EP6392 host plants with 5 (T1) or 20 (T2) aphids per plant, or controls without aphids (CK), was determined on day 6 after infestation. Thirty seeds were sown on 6 June, 2015 and apterous adult aphids were transferred to the back of the first true leaf per seedling 10 days after sowing. Three replicate experiments, using three separate plants each replication, were conducted for each treatment (CK, T1, and T2).

DNA Isolation and SLAF Library Construction for High-Throughput Sequencing
Leaves from the two parent plants (JY30 and EP6392) and F 2 individuals were collected on day 8 after aphid infestation, and an equal amount (0.1 g per plant) of leaves from 50 resistant and 50 susceptible individuals from F 2 plants were pooled, frozen in liquid nitrogen, and used for DNA extraction. Total genomic DNA was prepared using the CTAB method, with a modified CTAB buffer (8.18 g NaCl and 2 g CTAB in a total volume of 100 mL of 20 mM EDTA, 100 mM Tris, pH 8.0). DNA concentration and quality were estimated using a BioPhotometer Plus spectrophotometer (Eppendorf, Hamburg, Germany) and by electrophoresis through 1% agarose gels. Two DNA pools, resistant (R-pool) and susceptible (S-pool), were constructed using DNA isolated from pooled leaves from 50 individuals each of the 1000 F 2 plants. DNA samples isolated from JY30 and EP6392 plant leaves, and the two DNA pools, were used for SLAF library construction and sequencing.
A pilot experiment was performed to establish the conditions required to achieve optimal SLAF yield, avoid repetitive SLAFs, and obtain an even distribution of SLAFs for maximum SLAFseq efficiency. We constructed the SLAF library based on the result of the pilot experiment. Genomic DNA from each sample was incubated with HaeIII and RsaI, T4 DNA ligase, ATP, and an RsaI adapter (all from New England Biolabs (NEB), Ipswich, MA, USA) at 37 • C. Then, the restrictionligation reaction solutions were diluted and mixed with dNTPs, Taq DNA polymerase (NEB) and primers containing barcode 1 for polymerase chain reaction (PCR). An EZNA R Cycle-Pure Kit (Omega, London, UK) was used to purify the PCR products, which were then pooled and incubated at 37 • C with MseI (NEB), T4 DNA ligase, ATP, and a Solexa adapter. After incubation, the reaction products were purified using a Quick Spin column (Qiagen, Hilden, Germany) and electrophoresed through a 2% agarose gel. After gel purification, DNA fragments (SLAFs, including adapter sequence indices and adaptors) of 264-364 bp were excised and diluted for paired-end sequencing . SLAF-seq was performed using the Illumina HighSeq 2500 platform, with samples prepared according to the Illumina sample preparation guide (Illumina, Inc, San Diego, CA, USA) at the Biomarker Technologies Corporation (Beijing, China).

Analysis of SLAF-seq Data
SNP_index is a marker association analysis method to determine differences in genotype frequencies between pooled samples (Abe et al., 2012), by calculation of the (SNP_index). The stronger the correlation between a marker and a trait, the closer the (SNP_index) value is to 1. In this study, 'M, ' 'P, ' 'aa, ' and 'ab' denote the female parent (susceptible, JY30), the male parent (resistant, EP6392), the R-pool and the S-pool, respectively. (SNP_index) values were calculated as follows: SNP_index(ab) = Mab/(Pab+Mab); SNP_index(aa) = Maa/(Paa+Maa); (SNP_index) = SNP_index(aa) -SNP_index(ab); where Maa indicates the depth of the aa population derived from M, Paa indicates the depth of the aa population derived from P, Mab indicates the depth of the ab population derived from M, and Pab indicates the depth of the ab population derived from P. The allelic frequency was calculated by Euclidean distance followed by Loess regression analysis which identifies regions in which QTL lies and generates a list of putative regions in the linked genomic segment.

RNA Isolation, Reverse Transcription, and qRT-PCR Analysis
Leaves from cucumber plants in CK, T1 and T2 treatments were collected on day 6 after aphid infestation, frozen in liquid nitrogen, and used for RNA extraction. Total RNA was isolated from leaves using TakaRa RNAiso Reagent kit (Dalian, China). RNA concentration and quality were estimated using a BioPhotometer Plus spectrophotometer (Eppendorf, Hamburg, Germany). Reverse transcription was performed using TakaRa PrimeScript RT Reagent Kit (Dalian, China).
Primers for amplification of 16 genes were designed using Primer Premier 5.0 and are listed in Supplementary Table S1. qRT-PCR was performed using a TaKaRa SYBR PrimeScript RT-PCR Kit (Dalian, China), according to the manufacturer's instructions. The cucumber actin gene was used as an internal standard and amplified with the following primers: forward: 5 -TCGTGCTGGATTCTGGTG-3 , and reverse: 5 -GGCAGT GGTGGTGAACAT-3 . Reactions were performed in 96-well plates and the PCR program consisted of 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s and 50-60 • C for 30 s. qRT-PCR analysis was performed on an iQ 5 multicolor real-time PCR detection system (Bio-Rad, USA). The relative expression levels of the genes were determined using the 2 − Ct method.

Analysis of Numbers of Aphids per Plant
As expected, EP6392 and JY30 plants differed greatly in their resistance to aphids; EP6392 showed significant resistance to aphid infestation, whereas JY30 was susceptible. Phenotypic manifestations of this included the fact that the leaves of infested EP6392 plants curled upward slightly, whereas those of infested JY30 plants were more obviously curled upward. Moreover, EP6392 plants grew and developed normally, whereas JY30 plants were severely stunted on day 8 after aphid infestation. The F 2 population consisted of both susceptible and resistant plants ( Table 1); consistent with the results of our previous study (Liang et al., 2015).
When EP6392 plants were infested with different numbers of aphids, the average number of aphids per plant on day 6 after infestation correlated with the number of aphids used to infest the plants. The most significant difference in the number of aphids per plant was observed between T1 and T2 treatment. The average number of aphids per seedling on day 6 after aphid infestation were 80.20 and 167.20 for T1 and T2 treatment, respectively ( Table 2).

Analysis of SLAF-seq Data and SLAF Tags
In this study, a total of 38 million sequence data reads were obtained, each of 80 bp (×2). The majority of bases (80.55%) were of high quality (scores of ≥30; Table 3). The number of SLAF sequences were 99,180, and 98,812 for JY30 and EP6392 plants, respectively. Average SLAF sequence depths were 18.13-fold in the parental samples and 45.83-fold and 48.44fold in the R-and S-pools, respectively (Table 3). SLAF tags (n = 102,033) were divided into three categories: polymorphic, non-polymorphic and repeat, according to analysis of allele frequencies and gene sequences differences. Polymorphic SLAF tags were identified (n = 5,471), indicating a polymorphism rate of 5.36% (Table 4).

Association Analysis
A total of 1,961 tags from among the 5,471 polymorphic SLAF tags were selected for QTL identification based on criteria of sequence depth >5-fold in parents and the genotype of one allele being derived from JY30 and the other from EP6392 (Table 4). Peak regions above the threshold value were defined as those where Loess fitted values were greater than standard deviations above the genome-wide median in the (SNP_index) plot. One candidate region associated with cucumber aphid resistance spanning 0.31 Mb  on Chromosome 5 (Chr5: 26,684,572-26,994,642, cucumber line 9930 reference genome assembly, version 2) was identified with average (SNP_index) values above the threshold value of 0.2122 (Figure 1). The candidate region contained 43 genes and 20 SLAF makers. Based on functional annotation, 16 of the 43 genes were selected as good candidates for association with aphid resistance in cucumber ( Table 5). Details of the 20 SLAF makers are provided in Supplementary  Table S2.
Conversely, the expression levels of four genes, encoding stress-induced protein KIN2 (Csa5M643240.1 and Csa5M643260.1), probable LRR receptor-like serine/threonineprotein kinase At5g15730-like (Csa5M642150.1) and F-box family protein (Csa5M643280.1) were significantly higher in T1 and T2 treatments than in CK. Moreover, the expression levels of Csa5M643240.1 and Csa5M643260.1 were best significantly higher in T2 treatment than in T1 treatment. The expression level of the gene encoding F-box family protein (Csa5M643280.1) in T2 treatment was best significantly higher than in CK, and    was higher than in T1 treatment. The expression level of the gene encoding F-box/kelch-repeat protein (Csa5M643290.1) in T2 treatment was best significantly higher than in T1 treatment and in CK, and expression of this gene level was higher in T1 treatment than in CK; however, the difference was not statistically significant. Hence, these five genes were up-regulated after aphid infestation.
The expression of genes encoding probable E3 ubiquitinprotein ligase HERC4-like (Csa5M641620.1) and E3 ubiquitinprotein ligase At3g02290-like (Csa5M642120.1) did not appear to be regulated in response to aphid infestation. In addition, the expression levels of genes encoding protein NSP-interacting kinase 1-like (Csa5M642730.1) and Glycine-rich cell wall structural protein 1.0 (Csa5M643350.1) in T1 and T2 treatments were not significantly different from those in CK (Table 6; Figure 2).

DISCUSSION
Aphid is one of major pests affecting cucumber production. Severe aphid infestation can cause many visible symptoms, including leaf curling and yellowing, and plant wilting and stunting. In another aspect, the virus propagate with extreme rapidity by aphid. Understanding of the genetic mechanisms underlying aphid resistance is the most effective way to decrease aphid damage and produce improved quality cucumbers. SLAFseq technology is an effective and high-resolution technique for fine mapping of QTLs. The calculation of the SNP-index allows accurate quantitative evaluation of the frequencies of parental alleles, as well as the genomic contribution from the two parents, in F 2 individuals. The combination of SLAF-seq technology, the SNP-index method and BSA provides an efficient method to identify genomic regions related to traits. In this study, we used this approach to analyze two pooled F 2 population samples and detect a genomic region associated with aphid resistance in cucumber. A main effect QTL of 0.31 Mb on Chr5, including 43 genes, was obtained by association analysis. Sixteen of the 43  genes were identified as candidates for association with aphid resistance based on gene annotation. Plant hormones have important roles in signaling, and are vital in plant defenses against pathogens (Pieterse et al., 2009). The expression levels of a gene encoding the auxin transport protein BIG-like, was down-regulated in cucumber plants in response to aphid infestation. When green peach aphid (Myzus persicae Sulzer) infest Arabidopsis thaliana, the synthetics of Salicylic acid (SA) is promoted (Pegadaraju et al., 2005). Moreover, sugars are important in the defense of plants against pest infestation. Genes associated with sugar metabolism in Apium graveolens, Arabidopsis thaliana and Nicotiana attenuate are upregulated after aphid infestation (Moran et al., 2002;Heidel and Baldwin, 2004;Voelckel et al., 2004;Divol et al., 2005;Qubbaj et al., 2005;Thompson and Goggin, 2006). However, in this study, the expression level of the gene encoding endo-1,3(4)-beta-glucanase 1-like was down-regulated in cucumber plants in response to infestation, which differs to the findings of previous studies in other plant species. Genes encoding the transcription factors MYB59-like, the transcription factor HBP-1a-like, the F-box/kelch-repeat protein At5g15710-like, the betaglucan-binding protein and the proline-rich receptor-like protein kinase PERK10-like were also down-regulated in plants infested with different numbers of aphids. Our results indicate that aphid infestation resulted in repression of these seven genes and that they may be associated with aphid resistance in cucumber ( Table 6).
E3 ubiquitin ligase can function as a regulator of plant disease-resistance signaling in Arabidopsis thaliana (Stegmann et al., 2012), rice (Park et al., 2012), and tobacco (Gilroy et al., 2011). In a study revealing the mechanism of interaction between the rice blast fungus effector factors AvrPiz-t and relative R gene Piz-t (Park et al., 2012) found that APIP6, which interacts with the effector AvrPiz-t, is a RING type E3 ubiquitin ligase, which positively influences the immune response to rice blast fungus. Stegmann et al., (2012) found that E3 ubiquitin ligase can negatively regulate immune responses induced by pathogen associated molecular patterns (PAMPs). In this study, the expression levels of genes encoding the probable E3 ubiquitin-protein ligase HERC4-like and E3 ubiquitin-protein ligase At3g02290-like, were not obviously regulated in plants infested with different numbers of aphids, inconsistent with previous reports in other plant species.
R genes with LRR domains have major roles in regulation of the resistance of plants to pathogens and insects. Such genes confer resistance to the blue alfalfa aphid in Medicago truncatula Gaert (Klingler et al., 2005) and to soybean aphid in soybean (Kim et al., 2010a,b). The expression level of the gene Csa5M642150.1, encoding probable LRR receptor-like serine/threonine-protein kinase At5g15730-like, was up-regulated in cucumber plants in response to infestation with aphids. The expression levels of the genes encoding F-box family protein, F-box/kelch-repeat protein and the genes (Csa5M643240.1 and Csa5M643260.1) encoding stress-induced protein KIN2 were also up-regulated in plants in response to aphid infestation. These results indicate that aphid infestation enhanced expression of these five genes, suggesting that they could play an important role in aphid resistance in cucumber.

AUTHOR CONTRIBUTIONS
XC, FZ, XQ, and DL conceived the experiment and made the revision of the manuscript. DL and MC performed the research. DL and MC collected data. DL, MC, and QX analyzed the date and wrote the manuscript. All authors reviewed and approved this submission.