• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Evaluation of the early defoliation trait and identification of resistance genes through a comprehensive transcriptome analysis in pears

    2023-01-06 10:57:38SHANYanfeiLlMengyanWANGRunzeLlXiaogangLlNJingLlJiamingZHAOKejiaoWUJun
    Journal of Integrative Agriculture 2023年1期

    SHAN Yan-fei ,Ll Meng-yan ,WANG Run-ze ,Ll Xiao-gang ,LlN Jing ,Ll Jia-ming ,ZHAO Ke-jiao,WU Jun

    1 College of Horticulture/State Key Laboratory of Crop Genetics and Germplasm Enhancement,Nanjing Agricultural University,Nanjing 210095,P.R.China

    2 Institute of Horticulture,Jiangsu Academy of Agricultural Sciences/Jiangsu Key Laboratory for Horticultural Crop Genetic Improvement,Nanjing 210095,P.R.China

    Abstract Early defoliation,which usually occurs during summer in pear trees,is gradually becoming a major problem that poses a serious threat to the pear industry in southern China.However,there is no system for evaluating the responses of different cultivars to early defoliation,and our knowledge of the potential molecular regulation of the genes underlying this phenomenon is still limited.In this study,we conducted field investigations of 155 pear accessions to assess their resistance or susceptibility to early defoliation.A total of 126 accessions were found to be susceptible to early defoliation,and only 29 accessions were resistant.Among them,19 resistant accessions belong to the sand pear species (Pyrus pyrifolia).To identify the resistance genes related to early defoliation, the healthy and diseased samples of two sand pear accessions,namely,the resistant early defoliation accession ‘Whasan’ and the susceptible early defoliation accession ‘Cuiguan’,were used to perform RNA sequencing.Compared with ‘Cuiguan’,a total of 444 genes were uniquely differentially expressed in‘Whasan’.Combined with GO and KEGG enrichment analyses,we found that early defoliation was closely related to the stress response.Furthermore,a weighted gene co-expression network analysis revealed a high correlation of WRKY and ethylene responsive factor (ERF) transcription factors with early defoliation resistance.This study provides useful resistant germplasm resources and new insights into potentially essential genes that respond to early defoliation in pears,which may facilitate a better understanding of the resistance mechanism and molecular breeding of resistant pear cultivars.

    Keywords: pear,early defoliation,resistance and susceptibility evaluation,RNA-seq,WRKY,ERF

    1.lntroduction

    Pear is a globally cultivated fruit that belongs to thePyrusgenus of the Amygdaloideae subfamily in the Rosaceae family (Wuet al.2013).The annual pear production and harvested area have reached 23.1 million tons and 1.3 million ha,respectively.China ranks first in pear production worldwide,with an annual pear production of 16.1 million tons (FAO 2020).Pears have been cultivated in East Asia,North America,and Europe for over 3 000 years (Baudryet al.1993).The breeding of preferred pear accessions usually takes around 13 to 15 years and requires a large population volume,labor force,acreage,and financial investment.Therefore,the long-term and high investment characteristics of pear cultivation underscore the importance of optimizing pear quality and yield.

    The genusPyrusis susceptible to various diseases and pests,which reduce pear production and fruit quality.Common diseases include pear scab,pear rust,black spot disease,and fire blight.Pear scab,caused byVenturia pirina(Aderh.),is the most significant pear disease,and it severely reduces pear yield and fruit quality (Jamaret al.2017).Pear rust is a foliar pathogen of fruiting trees that is induced by the fungusGymnosporangium sabinae(Dicks).It causes yellow to bright orange leaf spots on pear tree leaves and mummifies the young fruit,resulting in complete defoliation and crop loss (Kaliszet al.2015).Black spot disease,caused byAlternaria alternata(Fries),is one of the most harmful diseases in Japanese pear cultivation.It causes severe bitterness in mature fruit,as well as complete defoliation (Tanahashiet al.2016).Fire blight,caused by the bacteriumErwinia amylovora,is widespread throughout Europe and Asia.This disease is devastating to various pomaceous trees and can induce the death of plants and loss of fruit production (Mertogluet al.2020),as well as partial defoliation in midsummer(Reighardet al.2005).A major pest in European pear production areas is the pear psylla (Psylla pyriL.).The nymphs and adults of psylla feed on sap from young leaves,shoots,and buds,causing direct damage to plants in summer (Erler 2004).During feeding,the honeydew produced by psylla adults and nymphs accumulates on the surfaces of leaves and fruit,leading to russeting disease(Pehlevan and Kovanci 2018).In total,this pest causes a decline in fruit quality and the deterioration of pear plantations (Li and Tian 2020).

    In recent years,early defoliation has increasingly threatened the sustainable development of the pear industry in South China.Early defoliation usually happens during mid-June to late-August,so it occurs approximately three months earlier than the natural fall (Baudryet al.1993).Some results have shown that early defoliation is related to ascendant trends in disease severity (Challiolet al.2006;Bogoet al.2018).Entomosporiumis a widespread plant disease(Goncalveset al.2013) that rapidly develops during plant growth.It was shown to cause early defoliation,eventually leading to a pear yield reduction in Brazil (Bogoet al.2018).Pear leaf spot is a severe disease occurring in many areas of the world where pears are cultivated under warm and humid conditions,and it causes the early defoliation of pear trees in mid-summer (Bell and van der Zwet 2005).In addition,the diseases caused byAlternaria alternata,Colletotrichum acutatum,andColletotrichum fructicolaalso cause severe early defoliation of young leaves during summer (Fukaya 2004;Tanahashiet al.2004;Zhanget al.2015).Aside from disease infection,the honeydew excretion ofPseudococcus comstocki(an invasive pest) and development of sooty mold can also result in early defoliation (Pellizzariet al.2012).In addition to biotic stimuli,abiotic stresses also induce early defoliation (Linet al.2018;Gortariet al.2021).The long-term use of chemical fertilizers aggravates the occurrence of early defoliation in pear cultivation areas in southern China (Kanget al.2021).The high contents of NO3-N and NH4-N in soil can cause leaf blackening,and early defoliation occurs in‘Niitaka’ pear (Song 2002).Furthermore,the characteristics of pear accessions determine their resistance to early defoliation (Matsumotoet al.2008).The occurrence of early defoliation will induce flower and bud necrosis (Yimet al.2005;Wennekeret al.2011),as well as the deterioration of vegetative-reproductive variables and the soluble solid content in pears (Hudina and Stampar 1999;Cheolet al.2015).Early defoliation has resulted in pear production losses in most areas,especially in southern China (Huanget al.2010) and southern Brazil (Bogoet al.2018).It is a major concern among many fruit farmers due to its negative effects on pear tree vigor,fruit quality,and sustained high yield.However,there is no system for evaluating the resistance and susceptibility of different accessions to early defoliation,although this is the important information for both cultivation protection and future breeding in pear.In addition to pears,this unfavorable phenomenon is also widespread in poplar(van den Berget al.2007;Gortariet al.2018),kiwifruit vines(Cruz-Castilloet al.2010;Kwacket al.2014),and rubber trees (Linet al.2018).

    Significant efforts have been made to identify the genes that are potentially related to disease defense and stress responses in pears.The scab resistance geneRvp1,identified in European pear ‘Navara’,causes resistance to the fungusV.pirina(Bouvieret al.2012).PbrATG8c,an autophagy marker in other species,could respond to a variety of tested stressors inPyrus bretschneideriRehd.such as dehydration,salt,and pathogen infection (Sunet al.2021).In the scab-inoculated leaves of Japanese pear,the activity of phenylalanine ammonia lyase (PAL)was enhanced.Salicylic acid (SA) also accumulated,which indicated that PAL might be related to SA and play a role in enhancing resistance (Faizeet al.2004;Shiet al.2021).High-throughput sequencing has gradually become the mainstream method for determining candidate stress-response genes (Liu Det al.2020;Wang Y Fet al.2020;Yuanet al.2020;Zhanget al.2020;Zhuet al.2020).RNA-seq was conducted on sand pear leaves that were inoculated or mock-inoculated with pear black spot disease,and four genes were subsequently selected as candidates for resistance to pear black spot disease in sand pears(Yanget al.2015).The results of RNA-seq in Nakai (Pyrus pyrifolia) showed that the pear responded to the infection ofPenicillium expansumby signal transduction,including a reactive oxygen species burst,plant hormone signal transduction,and transaction factor regulation (Xuet al.2021).Many studies on pears have demonstrated that transcription factors (TFs),including NAC,MYB,bHLH,WRKY,and AP2/ERF,play key roles in the regulation of biotic and abiotic stresses (Huanget al.2015;Liet al.2015,2018;Gonget al.2019;Gaoet al.2020;Wang Z Yet al.2020;Donget al.2021;Wang Y Set al.2021).These TFs systematically regulate the expression of target genes by binding tocis-acting elements in the gene promoters.Although a massive amount of research has been conducted to reveal the genes related to certain diseases and stress in pears,research related to early defoliation in pear has rarely been conducted.

    In this study,field investigations were performed to evaluate the resistance or susceptibility to early defoliation among a total of 155 pear accessions,which will facilitate identifying the resistant pear accessions.To further identify the candidate genes that respond to early defoliation,two sand pear accessions,the resistant accession ‘Whasan’and the susceptible one ‘Cuiguan’,were selected for further transcriptome analysis.RNA-seq analysis and weighted gene co-expression network analysis (WGCNA) were conducted,and the differentially expressed genes related to early defoliation resistance were finally selected.This study provides selected resistant pear accessions and new insights into the potential genes underlying the early defoliation phenomenon in pears,which will provide resistant accessions and genetic resources for the breeding of new resistant cultivars.

    2.Materials and methods

    2.1.Field investigations and plant material collection

    The field investigations of 155 accessions were conducted at the pear germplasm resource nursery of Jiangsu Academy of Agricultural Sciences on September 2nd,2020.For each accession,a total of three to five perennial branches from two or three trees were assessed using the disease index.The proportion of diseased areas (S) on leaves was evaluated firstly by Image J Software (Zhouet al.2020).The disease severity response to early defoliation was classified into five levels according to the S value: level 0: S=0;level 1:0<S≤0.1;level 3: 0.1<S≤0.25;level 5: 0.25<S≤0.4;level 7:0.4<S≤0.65;and level 9: 0.65<S≤1,and some examples of leaves belonging to the different levels are shown in Fig.1.Then,the disease index (DI) was calculated according to the following formula (Niu and Huo 2015):

    Fig.1 The status of leaves assigned to the different levels of disease severity.S,the proportion of diseased areas on leaves.

    According to the calculated DI values,the pear accessions were classified as follows: disease-resistant:0<DI≤40%;and disease-susceptible: 40%<DI≤100% (Liet al.2017).

    Based on the evaluation results,resistant and susceptible leaves were collected from ‘Whasan’ (P. pyrifolia) and‘Cuiguan’ (P. pyrifolia),respectively.For each accession,eight spotted (diseased) leaves and eight spot-free (healthy)leaves were collected for comparisons and photographic documentation.Then,the leaves and the basal portions of the petioles were separately collected and stored at -80°C.

    2.2.RNA isolation and sequencing

    Total RNA was extracted from four samples of ‘Cuiguan’(i.e.,diseased leaves: ‘CG_DL’;diseased basal: ‘CG_DB’;healthy leaves: ‘CG_HL’;healthy basal: ‘CG_HB’)and four samples of ‘Whasan’ (i.e.,diseased leaves:‘WS_DL’;diseased basal: ‘WS_DB’;healthy leaves:‘WS_HL’;healthy basal: ‘WS_HB’).The RNA integrity number (RIN) was determined using 1% agarose gels and an Agilent 2100 bioanalyzer (Agilent Technologies,California,USA).The purity of RNA was determined using a NanoPhotometer spectrophotometer (Thermo Fisher Scientific,Massachusetts,USA).The eight samples described above were sequenced,and separate libraries for ‘CG_DB’,‘CG_DL’,‘WS_DB’,‘WS_DL’,‘CG_HB’,‘CG_HL’,‘WS_HB’,and ‘WS_HL’ were prepared using the NEBNext?UltraTM RNA Library Prep Kit (New England Biolabs,USA).The first four libraries were set as the controls.Three replications were performed for each library.Fragmented mRNA and oligonucleotides were used for the synthesis of first-strand cDNA.Then,secondstrand cDNA synthesis was completed using dNTPs and DNA polymerase I.The 3′ ends of the purified cDNA were adenylated,and polyA was added and connected to the sequencing connector.After screening out 250–300 bp cDNA by using AMPure XP beads (Beckman Coulter,Beverly,USA),PCR amplification was performed.The product of the PCR procedure was subsequently purified,and the libraries were obtained.Library quality was assessed by a Qubit2.0 Fluorometer (Thermo Fisher Scientific,Massachusetts,USA).Then,the insert size of the libraries was determined,and qRT-PCR was performed to ensure library quality.Sequencing was conducted on the Illumina HiSeq 2500 according to the efficient concentrations and target data volume.

    2.3.Sequencing data analysis and differentially expressed gene (DEG) identification

    The quality of sequencing reads was evaluated using FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/).Raw data were filtered using Trimmomatic software v0.39 (Bolgeret al.2014) to remove the read adapters and other low-quality reads.High-quality data were assembled and mapped to the reference genome of pear (Wuet al.2013) using HISAT2 v2.0.1 (Kimet al.2015).For expression quantification,FeatureCounts Software v2.0.0 was used to count the mapped reads in each gene (Liaoet al.2014),and the values of reads per kilobase per million (RPKM) mapped reads were calculated by an in-house Python script.Using “DESeq2 v1.26.0” in R,the DEGs were defined by |log2FoldChange|≥1 and aP-value of ≤0.01 (Mortazaviet al.2008;Dobinet al.2013).

    2.4.Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses

    PlantRegMap (http://planttfdb.cbi.pku.edu.cn/) was used for GO enrichment analysis of the DEGs,and the KEGG pathway enrichment of the DEGs was analyzed using KOBAS v3.0 (http://kobas.cbi.pku.edu.cn/ kobas3/genelist/).The results of the GO enrichment analysis were plotted as bar charts in Microsoft Excel 2016,and the results of the KEGG enrichment analysis were plotted as bubble charts using “ggplot2 v3.3.5” in R.

    2.5.ldentification of TFs and plant disease resistance genes related to disease and stress

    To predict the TFs,the protein sequences of the DEGs in the comparison groups were obtained from the pear protein database using TBtools v1.0 (Chenet al.2020),and subsequently submitted to the plant TF database(PlantTFDB v5.0;http://planttfdb.gao-lab.org/prediction.php).Eventually,the classified TF groups with an E-value of≤1e–10were obtained.To predict the putative plant disease resistance genes,we submitted the protein sequences of the DEGs to the Pathogen Receptor Genes database (PRGdb)v4.0 (http://prgdb.org/prgdb4/) (Calle Garciaet al.2021).Plant disease resistance genes were screened out with an E-value of ≥1e–6using the blast function.Moreover,the domains and classes of these genes were analyzed using DRAGO3 tools in PRGdb v4.0.

    2.6.Weighted gene co-expression network analysis(WGCNA)

    WGCNA was performed according to our previous study(Wang R Zet al.2021).First,we merged the gene expression matrices of all 24 samples.Then,we calculated the adjacencies in WGCNA with a best soft power of 19,and minimized the effects of noise and spurious associations using the Topological Overlap Matrix (TOM)(Yip and Horvath 2007).Subsequently,the co-expression relationships between the DEGs were obtained from TOM with a weight threshold of 0.02.Moreover,the coexpression relationships containing TFs and plant disease resistance genes,as well as stress-related genes,were extracted and visualized using Cytoscape v3.6.0 (Shannonet al.2003).

    2.7.Validation of gene expression by quantitative real-time PCR (qRT-PCR)

    The total RNA was extracted from the eight samples described above using a Plant Total RNA Isolation Kit(FORGENE,Chengdu,China),strictly following the manufacturer’s instructions.The concentrations were tested using Nanodrop and Agilent 2100.qRT-PCR was performed to verify the expression results of the transcriptome sequencing.The specific primers used to amplify the DEGs were designed using the Primer-Blast function in the National Center for Biotechnology Information (NCBI) Search database (https://www.ncbi.nlm.nih.gov/) and were procured from Tsingke Biotechnology Co.,Ltd.(Beijing,China).The primer sequences are shown in Appendix A.RNA samples were used to synthesize cDNA using TransScript?One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech,Beijing,China).qRT-PCR amplification was completed using Genious 2X SYBR Green Fast qPCR Mix (ABclonal Technology,Wuhan,China) on LightCycler 480 (Roche,Basel,Switzerland).The PCR procedure was as follows: pre-incubation at 95°C for 5 min;amplification with 45 cycles of 95°C for 3 s,60°C for 10 s,and 72°C for 3 s;and final cooling at 40°C for 30 s.Three biological samples and three technical repeats were evaluated to confirm the reliability of the results.PbrGADPHwas used as a reference gene.The relative expression level of each gene was calculated using the 2–ΔΔCtmethod (Livak and Schmittgen 2001).

    2.8.Data analyses

    Data processing was performed using the R package and Microsoft Excel 2016.The stacked bar diagram,Venn diagrams,and heatmaps were drawn using Origin 2019,Venn diagram in the R package,and TBtools v1.0,respectively.The visualized charts of RPKM values and relative expression levels were drawn using GraphPad Prism v6.01,and the correlation coefficient between them was analyzed using the CORREL function in Microsoft Excel 2016.

    3.Results

    3.1.Evaluation of the resistance and susceptibility to early defoliation in the different pear accessions

    Based on the standards of disease severity evaluation and the DI calculation,126 of the 155 accessions were found to be susceptible to early defoliation,and only 29 accessions were resistant based on the field investigations.The distribution of DI approximated a normal distribution,the values varied from 18.34 to 97.86%,and were mainly distributed between 40 and 80%,which indicated that most of the pear accessions were disease-susceptible(Fig.2-A;Appendix B).For the disease severity evaluation,the different accessions showed varied resistance levels.Among the five levels,19.35% pear accessions had leaves in level 0,56.77% had leaves in level 1,85.16% had leaves in level 3;and all of the pear accessions included leaves of level 5 (Table 1).Further statistical analysis revealed that pear accessions showed resistance when their leaves ranked as level 0 to level 3,accounting for more than 60%of the total leaves evaluated (Appendix B).

    Table 1 The number of accessions distributed among the different levels of disease severity

    Fig.2 The disease index investigations of 155 accessions genus (Pyrus L.).A,frequency distribution of disease index values of the 155 accessions.B,the number of disease-resistant (pink) and disease-susceptible (green) accessions in the different cultivated species.P,P.pyrifolia;B,P.brestschneideri;U,P.ussuriensis;C,P.communis;and I,interspecific hybrid.

    Furthermore,the 155 accessions were classified into five cultivated species.There were clearly more diseasesusceptible accessions than disease-resistant accessions in each pear species (Fig.2-B).Among the 29 resistant accessions,19 accessions ofP. pyrifoliawere resistant to the disease,accounting for 65.52% of all resistant accessions.In addition,four resistant accessions (‘Xueli’,‘Huangxiang’,‘Jinhua’,and ‘Yali’) belong toP.brestschneideri,and five resistant accessions (‘Jimi’,‘Eli No.2’,‘Zhongai No.2’,‘Eli No.1’,and ‘Tienong No.25’) are interspecific hybrids.Only one resistant accession (‘Gaolici’) belongs toP.communis.Furthermore,the most resistant accession was ‘Atago’(P. pyrifolia),while ‘Pugua’,‘Whasan’,‘Gaoyou-1’,‘Long No.10’,‘Jinshui No.3’,‘Gaoyou-2’,and ‘Songdao’,which all belong toP. pyrifolia,were the top ten resistant accessions(Appendix B).All of these resistant accessions could be recommended as parents in future breeding projects.In addition,some of the mainly cultivated sand pear cultivars,such as ‘Cuiguan’,‘Wonhwang’,‘Nijisserki’,and ‘Housui’,exhibited extremely high disease index values,which indicated their high susceptibility to early defoliation and unsuitability for resistance breeding.

    To explore the underlying mechanism of the resistance differences,the disease-resistant accession ‘Whasan’(DI=21.37%) and disease-susceptible accession ‘Cuiguan’(DI=86.11%) were selected as representatives for further analysis (Appendix B).The status evaluations of whole trees showed that many more leaves remained on ‘Whasan’than on ‘Cuiguan’ (Fig.3-A).Heathy leaves of ‘Whasan’and ‘Cuiguan’,on which no disease spots were found,were selected as the controls.Comparing the diseased leaves of these two accessions,it was obvious that the diseased area of ‘Cuiguan’ was much higher than that of ‘Whasan’.In addition,the lesion areas showed symptoms of black spots and yellowing (Fig.3-B).The results above indicated significant differences between the resistant accession‘Whasan’ and the susceptible accession ‘Cuiguan’.

    Fig.3 The phenotypes of leaves from ‘Cuiguan’ and ‘Whasan’.A,the conditions of whole trees.B,the symptoms of leaves.

    3.2.RNA-seq of two pear accessions with contrasting disease resistance

    The diseased leaves (DLs),the basal petioles of diseased leaves (DB),the healthy leaves (HLs),and the basal petioles of healthy leaves (HB) of ‘Cuiguan’ (CG) and ‘Whasan’(WS) were collected for sequencing.After RNA-seq,the quality of eight libraries from the 24 samples (three biological replicates were collected for each tissue in the two accessions) was assessed (Appendix C).The number of raw reads for these libraries ranged from 40.7 to 52.9 million.After trimming the reads,a total of 39.1 to 51.4 million high-quality reads and 58.7 to 77.2 billion bases per library were obtained.The average Q20 (sequencing base calls with an error rate of <4%) was more than 98%,and the average GC content (proportion of guanine (G) and cytosine(C) in the nucleotide sequence) was 46.82% (Appendix C).The amount of RNA sequencing data produced for each sample was at least six giga base pairs (Gbps).A total of 34 624 genes were obtained,excluding the genes with RPKM values of 0 in all samples (Appendix D).These results indicated that the sequencing quality was high enough for further analysis.

    To investigate the relationships between gene expression and resistance differences,the DEGs in ‘Whasan’ and‘Cuiguan’ were analyzed.The number of DEGs in ‘WS_DL’vs.‘WS_HL’ was the highest (6 262 DEGs).A total of 3 316 DEGs were up-regulated,and 2 946 DEGs were downregulated in this comparison group (Fig.4-A).In contrast,the number of DEGs was the least in ‘CG_DB’vs.‘CG_HB’(441 DEGs).Compared with the healthy basal group,114 and 327 genes were up-regulated and down-regulated,respectively,in the diseased basal sample of ‘Cuiguan’.To further identify the candidate DEGs related to early defoliation,an overlapping analysis was conducted.We speculated that DEGs related to resistance mainly existed in the comparison groups ‘WS_DBvs.WS_HB’ and ‘WS_DLvs.WS_HL’.Finally,444 DEGs were shared in ‘WS_DBvs.WS_HB’ and ‘WS_DLvs.WS_HL’ but not in ‘CG_DBvs.CG_HB’ or ‘CG_DLvs.CG_HL’ (Fig.4-B).To identify the expression patterns of the DEGs,a further cluster analysis was conducted.The expression patterns of DEGs between ‘WS_DB’ and ‘WS_HB’ and between ‘WS_DL’and ‘WS_HL’ showed significantly opposite expression trends (Fig.4-C),indicating the different expression profiles between the healthy and diseased tissues.In addition,most DEGs showed different expression trends between basal petioles and leaves,which suggests that these may be tissue-specific DEGs.Meanwhile,the statistical data showed that 126 DEGs with the highest expression level were centralized in ‘WS_HB’,while only 91 DEGs with the highest expression level were distributed in ‘WS_DB’(Appendix E).

    Fig.4 Characteristics of differentially expressed genes (DEGs) related to early defoliation.A,stacked bar diagram showing the number of DEGs among the four comparison groups.Four samples from ‘Cuiguan’ (diseased leaves: ‘CG_DL’;diseased basal:‘CG_DB’;healthy leaves: ‘CG_HL’;healthy basal: ‘CG_HB’) and ‘Whasan’ (diseased leaves: ‘WS_DL’;diseased basal: ‘WS_DB’;healthy leaves: ‘WS_HL’;healthy basal: ‘WS_HB’).All DEGs (pink),up-regulated DEGs (yellow) and down-regulated DEGs (green)are quantified.B,Venn diagram showing the overlap of DEGs among the comparison groups (‘diseased leaves/basal petiole’vs.‘healthy leaves/basal petiole’ for both ‘Cuiguan’ and ‘Whasan’).C,hierarchical clustering of 444 DEGs uniquely expressed in ‘Whasan’ (‘WS_DL’ vs.‘WS_HL’ and ‘WS_DB’ vs.‘WS_HB’).Red represents high expression levels of each gene in certain tissues,and blue represents low expression levels.

    3.3.GO enrichment analysis of DEGs in response to early defoliation

    The potential functions of the DEGs were further categorized into three groups and 45 enriched terms using GO enrichment analysis (Fig.5-A;Appendix F).In the biological processes,terms related to secondary metabolic,biosynthetic,and regulation processes,such as phenylpropanoid biosynthetic process and chalcone biosynthetic and metabolic process,were mainly enriched (Fig.5-A).Genes related to the reproductive processes (23 genes) and the multi-organism processes (22 genes) were mostly enriched.In contrast,the number of genes associated with chalcone metabolic and biosynthetic processes,regulation of anthocyanin biosynthetic processes,proanthocyanidin biosynthetic processes,and regulation of flavonoid biosynthetic processes was the lowest (two genes) (Appendix F).Moreover,genes related to the regulation of hormone level,response to wounding,and biotic stimulus were also significantly enriched (Fig.5-A),and these are all related to abiotic stress and disease incubation in leaves (Cheonget al.2002;Xuet al.2015).Only a few genes were enriched in protein ubiquitination and protein modification.For the molecular function,the catalytic activity,naringeninchalcone synthase activity,and transferase activity were minor categories.For cellular components,the cell wall,external encapsulating structure,cell periphery,membrane,and exocyst were significantly enriched terms (Fig.5-A).

    Fig.5 Gene Ontology (GO) functional enrichment analysis of 444 differentially expressed genes (DEGs) uniquely expressed in‘Whasan’.A,GO enrichment of DEGs.Terms marked by red boxes are further analyzed.B,expression profiles of DEGs selected from the red boxes in Fig.5-A.Hierarchical clustering was performed to classify DEGs with similar expression patterns based on normalized reads per kilobase per million (RPKM) values of the DEGs.The blue to red represents the expression levels from low to high.The genes marked in red are the focus.Four samples from ‘Whasan’ (diseased leaves: ‘WS_DL’;diseased basal:‘WS_DB’;healthy leaves: ‘WS_HL’;healthy basal: ‘WS_HB’).

    We also analyzed the expression levels of the genes involved in the regulation of hormone level,response to wounding,and biotic stimulus.It was obvious that half of these DEGs (12 DEGs) showed their highest expression levels in the diseased basal petiole of ‘Whasan’,which included DEGs related to spermine synthase,choline/ethanolamine kinase,chalcone synthase,and thaumatinlike protein,among others.In addition,genes encoding cytokinin dehydrogenase showed higher expression levels in healthy ‘Whasan’ leaves.Furthermore,two WRKY TFs,WRKY11 and WRKY22 (Pbr018160.1andPbr032698.1),were highly up-regulated in healthy basal petioles and leaves,respectively (Fig.5-B;Table 2),and they are widelyknown to be related to disease resistance and abiotic stress responses (Meiret al.2010;Aliet al.2018).

    Table 2 The basic information for the 37 candidate differentially expressed genes

    3.4.KEGG enrichment analysis of DEGs in response to early defoliation

    Subsequently,the KEGG enrichment analysis showed that 69 pathways were enriched.Among them,13 pathways were apparently enriched,with aP-value of ≤0.05 as the threshold (Fig.6).It is significant that flavonoid biosynthesis,the biosynthesis of secondary metabolites,and plant–pathogen interaction were the most abundant.Genes related to metabolic pathways were mostly enriched in 57 genes.Terms annotated as plant–pathogen interactions were significantly enriched in response to stress (Kahlon and Stam 2021).In addition,other pathways,including metabolic pathways,phenylpropanoid biosynthesis,circadian rhythmplant,zeatin biosynthesis,glutathione metabolism,sulfur metabolism,MAPK signaling pathway-plant,pentose and glucuronate interconversions,and other glycan degradation,were enriched,and circadian rhythm-plant was also enriched(Fig.6-A;Appendix G).

    To further clarify the cause of early defoliation,we focused on the significantly enriched pathway ‘plant–pathogen interaction’.In this enriched term,genes encoding calmodulin-like proteins (CMLs) were the most common(five genes),and they have been reported to have functions in stress response (Vadasseryet al.2012).Among them,the expression levels ofCML3(Pbr026767.2),CML41(Pbr006192.1),andCML42(Pbr028492.1) increased in the diseased leaves of ‘Whasan’.In contrast,the expression levels of two otherCMLgenes (Pbr021489.1andPbr013669.1) were higher in the healthy leaves of‘Whasan’.In addition,two stress-related TFs (Hsuet al.2013;Jeonet al.2016;Hussainet al.2018),ethyleneresponsive TF (Pbr031830.1) and WRKY TF (Pbr032698.1),were up-regulated in the healthy leaves and basal petioles,respectively (Fig.6-B).

    Fig.6 Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of 444 differentially expressed genes (DEGs)uniquely expressed in ‘Whasan’.A,KEGG pathway enrichment of DEGs.Red represents a higher significance and blue represents a lower significance.The size of the bubbles indicates the number of genes enriched in each pathway. Terms marked by a red box are further analyzed.B,expression profiles of DEGs selected from the red box in Fig.6-A.Hierarchical clustering was performed to classify the DEGs with similar expression patterns based on normalized reads per kilobase per million (RPKM) values of the DEGs.The blue to red represents the expression levels from low to high.The genes marked in red are on the focus.Four samples from ‘Whasan’ (diseased leaves: ‘WS_DL’;diseased basal: ‘WS_DB’;healthy leaves: ‘WS_HL’;healthy basal: ‘WS_HB’).

    3.5.ldentification of DEGs functioning in defense responses and construction of a co-expression network

    To comprehensively screen out TFs associated with early defoliation,DEGs from the GO and KEGG analyses were further submitted to PlantTFDB v5.0.The results showed that the genesPbr001424.1,Pbr018160.1,andPbr032698.1were identified as WRKY TFs,whilePbr031830.1was identified as an ERF TF,which was identical with the annotation (Table 2).Interestingly,WRKY and ERF TFs were previously reported to respond to several stressors through an intricate regulatory network in plants (Merajet al.2020).Therefore,to further understand the regulatory mechanism related to the stress responses in early defoliation,a coexpression network was then constructed centering on these four TFs,with the 444 DEGs as the background.Among the 444 DEGs,plant disease resistance genes were identified using PRGdb v4.0.Filtering out the DEGs with a low co-expression weight threshold of 0.02,266 mutual coexpression relationships were identified (Fig.7;Appendix I).A total of 175 DEGs were co-expressed with these four TFs (purple nodes),and 19 plant disease resistance genes (pink nodes) and 40 stress-related genes (orange nodes) were identified among them.It was observed that the node sizes ofWRKY11(Pbr018160.1) andWRKY22(Pbr032698.1) were the largest,indicating that they had the most connections with other DEGs.Most of the DEGs thatWRKY11andWRKY22correlated with were stress related (22 DEGs),and several DEGs were plant disease resistance genes (six DEGs).There was also a strong connection betweenWRKY11andWRKY22(red curved line) (weight value=0.34).Another WRKY TF,WRKY70(Pbr001424.1),was mostly correlated with plant disease resistance genes (13 DEGs) and several stress-related genes (two DEGs).Furthermore,an ERF family memberCRF4was mostly correlated with stress-related genes(20 DEGs).These results suggest that these four TFs are closely related to plant disease defense and stress responses and may play important roles in the response to early defoliation.In addition,53 plant disease resistance genes were predicted in the resistant accession ‘Whasan’(Appendix H),representing 11.94% of the unique ‘Whasan’DEGs.These resistance genes shared seven domains.Most of the genes (19 DEGs) contained transmembrane domains,and only one DEG belonged to the coiled-coil domain.Furthermore,the leucine-rich repeat domain (10DEGs),lectin-like motifs (seven DEGs),nucleotide-binding site (four DEGs) and kinase domain (three DEGs) were also identified (Table 3).According to their various domains,the DEGs were mainly divided into four classes,including the RLP (15 DEGs),KIN (eight DEGs),LECRK (seven DEGs),and RLK (seven DEGs) classes (Table 4).

    Fig.7 Co-expression network of differentially expressed genes centralized on the WRKY and ERF transcription factors.Cytoscape representation of the co-expressed genes with edge weights of >0.02.Purple nodes indicate transcription factors,pink nodes indicate plant disease resistance genes and orange nodes indicate stress-related genes.Transcription factors and plant diseased resistance genes,as well as stress-related genes,are connected via undirected lines.The red curved line indicates the interaction between WRKY11 and WRKY22.

    Table 3 Summary of predicted pear resistance genes

    Table 4 Summary of predicted resistance differentially expressed genes (DEGs)

    3.6.Validation of RNA-seq data by qRT-PCR

    As discussed above,our analyses identified some critical genes,including TFs,such as WRKY and ERF TFs,CMLs,and plant disease resistance genes.The expression levels of these genes were verified by qRTPCR.TFs play essential roles in regulating stress-related genes.The expression level ofWRKY70in ‘Cuiguan’was higher than that in ‘Whasan’,whileCRF4showed higher expression levels in ‘Whasan’.Additionally,in‘Whasan’,the expression level ofWRKY22in healthy leaves was higher than that in diseased leaves.CMLs are reported to be related to stress responses.In the present study,the expression levels ofCML44andCML42were higher in ‘Whasan’ compared to ‘Cuiguan’.In addition,compared with healthy leaves,CML3andCML41were up-regulated in the diseased leaves in ‘Whasan’.The expression levels of plant disease resistance genes were also verified.In ‘Whasan’,LeEIX2andRFO1were significantly up-regulated in the diseased leaves compared with the healthy leaves.Various expression levels of the abovementioned genes in different disease-resistance accessions and tissues verified the importance of those genes in regulating the response to stress,which was related to early defoliation.Moreover,we compared the results of qRT-PCR with the RPKM values of these genes.The correlation coefficient between them was calculated by applying the correlation function,and the value was defined as “r”.The results show a high correlation between the transcriptome data and qRT-PCR (average “r”=0.71),indicating the reliability of the transcriptome data (Fig.8).

    Fig.8 qRT-PCR validation of the nine selected differentially expressed genes from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.The blue bars represent relative expression level corresponding to the left Y-axis.The orange lines represent reads per kilobase per million (RPKM) values from RNA-Seq data.“r” represents the correlation coefficient between RPKM and relative expression data.Four samples are from ‘Cuiguan’ (diseased leaves: ‘CG_DL’;diseased basal:‘CG_DB’;healthy leaves: ‘CG_HL’;healthy basal: ‘CG_HB’) and ‘Whasan’ (diseased leaves: ‘WS_DL’;diseased basal: ‘WS_DB’;healthy leaves: ‘WS_HL’;healthy basal: ‘WS_HB’).The error bars represent the standard deviations (n=3).

    4.Discussion

    4.1.The assessment of resistance and susceptibility to early defoliation in pear accessions

    In the present study,we performed field investigations on 155 pear accessions,and the DI was calculated to evaluate the resistance of these accessions to early defoliation(Appendix B).The results clearly show that the diseasesusceptible pear accessions (126 accessions) outnumbered the disease-resistant accessions (29 accessions) (Fig.2).Among the resistant accessions,19 belonged toP.pyrifolia,which indicated that sand pear accessions might be more resistant to early defoliation,which was similar to previous research (Liet al.2017).Accession ‘Whasan’ exhibited strong resistance to early defoliation (Appendix B;Fig.3).Similarly,a previous study reported that the symptom of early defoliation in ‘Whasan’ was only mild even during the rainy season (Wanget al.2016).Notably,sand pear‘Atago’ showed the strongest resistance among the 155 accessions (Appendix B).Previous reports indicated that‘Atago’ showed high resistance to diseases such as black spot and pear leaf spot (Chen 2008;Denget al.2008),which are the common diseases causing early defoliation (Bell and van der Zwet 2005;Tanahashiet al.2016).Furthermore,several accessions were found to show tolerance to early defoliation (Appendix B),such as ‘Pugua’,‘Gaoyou-1’ and‘Jinshui No.3’.Interestingly,the female parent of ‘Kousui’is ‘Kikusui’,and the disease index values of these two accessions are similar (Appendix B).Previous studies indicated that the resistance is an inheritable trait in pears(Bueset al.2000;Ohet al.2021).Therefore,we speculate that the resistance of ‘Kikusui’ could be passed on to the progeny.For the other resistant accessions evaluatedviathis field investigation,we can use their parents for mining more resistant germplasm resources,and all of the 29 identified resistant accessions can provide useful breeding resources for early defoliation resistance.

    Among the disease-susceptible accessions,‘Cuiguan’was one of leading cultivars in South China (Gaoet al.2021),and the occurrence of early defoliation was widespread after exposure to stress (Liuet al.2011;Huanget al.2018;Yanget al.2018).Consistent with those findings,‘Cuiguan’ was also evaluated as a disease-susceptible accession in our study.Another sand pear accession,‘Wonhwang’,showed susceptibility to early defoliation,which was consistent with previous research (Wanget al.2016).Other accessions,such as ‘Housui’,‘Nijisserki’,and ‘Niitaka’ (Appendix B),which have not been reported as resistant to early defoliation,were also evaluated as susceptible accessions in this study.Therefore,we propose that these disease-susceptible accessions should not be recommended as parents for resistance breeding,and effective measures must be taken to alleviate early defoliation in these susceptible accessions.Recently,bioorganic fertilizer treatment (BIO,2%,w/w)or BIO combined with humic acid (HA) was reported to be effective for alleviating the early defoliation in pear(Kanget al.2021).In addition,the fact that the proportion of resistant accessions,i.e.,those with disease severity evaluated as levels 0,1 and 3,was more than 60% will be instructive for evaluating resistant or susceptible accessions more directly during field investigations.

    4.2.Stress response-related TFs were induced by early defoliation

    The GO and KEGG enrichment analyses and WGCNA revealed a correlation between TFs and the DEGs enriched in stress response,indicating that the occurrence of early defoliation might be closely related to the stress response.It is well known that differential responses to many stressors in plants are carried out by an intricate regulatory network of transcription repressors and activators,including the common TFs AP2/ERF,WRKY,bHLH,bZIP,MYB,and NAC(Merajet al.2020).Early defoliation is induced by complex stimuli,and plants will have corresponding reactions through their regulatory circuits.The roles of WRKY TFs in plant defense have been reported,and they occupy a central position mediating fast,positive,and negative regulation of disease resistance (Meiret al.2010) and play a profound role in abiotic stresses (Aliet al.2018).Notably,WRKY TFs were identified by GO and KEGG enrichment analyses and WGCNA in the present study,which indicated that they might play a role in the pear tree response to early defoliation.WRKY22is known to positively regulate resistance to pathogens inArabidopsisand pepper,and the loss of function ofWRKY22enhances susceptibility to disease in these two plants (Hsuet al.2013;Hussainet al.2018).Similarly,WRKY11positively enhances the salt and drought stress tolerance of tobacco and rice (Xuet al.2014;Leeet al.2018).It was also demonstrated thatMsWRKY22directly binds to the W-box (elements in the promoter region ofMsWRKY11,[(T)(T)TGAC(T/C)]) and activatesMsWRKY11expression in alfalfa (Wenet al.2021).In the present study,the expression patterns ofWRKY11andWRKY22were similar,and they showed higher expression levels in healthy tissues (Table 2).Additionally,a strong correlation was found betweenWRKY11andWRKY22(weight value=0.34) (Fig.7;Appendix I) in the co-expression network.We also found that a W-box (TTGACC) exists in the promoter ofWRKY11(Appendix J) (PlantCARE,http://bioinformatics.psb.ugent.be/webtools/ plantcare/html/),which was similar to previous findings (Wenet al.2021).Thus,we speculate thatWRKY22might activateWRKY11by binding to the W-box,thereby enhancing plant resistance to stress,and leading to less early defoliation.AnotherWRKYmember,WRKY70,was found to have a higher expression level in diseased tissues.This result was also found in the shoot tissues of susceptible chickpea,whereCaWRKY70was significantly accumulated under Foc1 (Fusarium oxysporum) infection (Chakrabortyet al.2020).The different expression patterns in the present study indicated that WRKY TFs might function redundantly in plant defense.

    As plant-specific TFs,AP2/ERF TFs are widely involved in various biological processes,especially in abiotic stimuli tolerance (Liu X Met al.2020).Cytokinin response factors(CRFs) are a subfamily of AP2/ERF that are involved in stress responses,including oxidative stress,drought stress,and disease defense (Guet al.2002;Shiet al.2014;Hugheset al.2021).AtCRF2,a homologous gene ofCRF4inArabidopsis,positively regulates plant immunity and has a strong correlation with the cold response (Jeonet al.2016).The overexpression ofCRF4leads to a stronger tolerance to freezing upon exposure to freezing temperatures,while mutants lackingCRF4expression are more sensitive to freezing (Zwacket al.2016).In the present study,CRF4showed a higher expression level in ‘Whasan’,which might be the reason why that accession showed stress resistance,contributing to its tolerance to early defoliation.In addition,WGCNA showed thatCRF4was co-expressed with stressrelated genes,indicating its possible role in causing leaf resistance.

    4.3.Resistance genes of pears in response to early defoliation

    Recent studies have shown that the signaling pathway regulates defoliation (Xuet al.2019;Liet al.2020).It is widely known that Ca2+is a mediator that plays an important role in plant signaling in response to biotic and abiotic stresses (Kudlaet al.2010).CMLs are major Ca2+sensors controlling various cellular functions (Chevalet al.2013),and they play a wide array of roles in stress perception.In the present study,CML44(Pbr021489.1) andCML42(Pbr013669.1) exhibited up-regulated trends in healthy leaves,which were consistent with those inArabidopsis(Chiassonet al.2005).AnotherCML42(Pbr028492.1)showed the opposite expression trends,and was highly expressed in diseased leaves.This difference could be explained by multiple functions under different stimuli,as was shown in a previous study in whichCML42was reported to function as a negative regulator in herbivory defense and as a positive regulator in the defense against UV-B stress(Vadasseryet al.2012).In addition,the co-expressed network showed various relationships betweenCMLs and TFs,including CRF4,WRKY70,and other TFs (Fig.7;Appendix I),which might lead to differentCMLfunctions for plant defense.A previous study showed that WRKY TFs were predicted asCMLtarget genes (Navarroet al.2004).Thus,it can be presumed that probable regulation or interaction relationships might exist between these two TF types and severalCMLs,resulting in different roles for theCMLs and subsequently influencing the disease resistance of plants.However,the specific correlations between other TFs andCMLs are unclear.Additionally,the direct connection between early defoliation and the Ca2+signaling pathways remains complex,and these issues are worthy of further exploration.

    Plants are usually attacked by a variety of pathogens,such as fungi,viruses,and bacteria.In this study,the “plant–pathogen interaction” term was significantly enrichedviaKEGG enrichment analysis,implying that early defoliation might be related to plant pathogens.In our study,53 plant disease resistance genes were uniquely identified in the resistant accession ‘Whasan’ using PRGdb v4.0 (Appendix H),which may be related to the greater resistance of ‘Whasan’.Resistance genes were developed during the long battle between plants and pathogens,and they present the plant with specific receptors and defense mechanisms (Andolfo and Ercolano 2015).These genes share some conserved domains and are categorized into various classes according to the arrangements of their domains (Calle Garciaet al.2021) (Tables 3 and 4).Leucine-rich repeats (LRRs),the major domain found in most classes of resistance genes,are essential for recognition specificity and resistance.For example,the polygalacturonase-inhibiting proteins (PGIPs)contain an LRR domain,which limits the growth of plant pathogens and induces the defense responses in plants (Diet al.2006).LeEix proteins belong to RLP (containing an LRR and a transmembrane domain),and initiate a defense response upon induction by ethylene-inducing xylanase(EIX) of non-pathogenicTrichodermain tomato (Kruijtet al.2005).However,there are receptors that contain domains other than LRRs.RFO1 is a receptor-like kinase with no extracellular LRR domain,and it encodes a new type of disease-resistance protein that imparts resistance to a broad spectrum ofFusariuminArabidopsis thaliana(Diener and Ausubel 2005).Pid2belongs to a new resistance class,LECRK,that contains lectin-like motifs instead of the LRR domain,and it elicits resistance to blast fungi in rice (Wanget al.2015).A previous study showed thatWRKYcould be involved in the process of plant defense by regulating resistance genes (Renet al.2021).In our study,these resistance genes were co-expressed with WRKY TFs(Fig.7).We predicted that the higher resistance of ‘Whasan’was not merely due to the existence of resistance genes but also due to the potential interactions between WRKY TFs and resistance genes.

    5.Conclusion

    In this study,we performed field investigations and assessed the early defoliation resistance of 155 pear accessions.A total of 126 accessions were susceptible to early defoliation and 29 accessions were resistant.Among them,we selected two sand pears (P. pyrifolia),a disease-resistant accession‘Whasan’ and a disease-susceptible accession ‘Cuiguan’,for a subsequent RNA-seq analysis.A total of 444 genes were uniquely differentially expressed in the resistant accession‘Whasan’,compared with the susceptible accession‘Cuiguan’.Further functional enrichment analysis of these DEGs confirmed that early defoliation was related to plant–pathogen interaction,wounding,biotic stimuli response,and hormone level regulation.Combined with WGCNA and enrichment analysis,some WRKY and ERF TFs were identified as potentially essential TFs in response to early defoliation.Meanwhile,the different expression levels of these genes,which were verified by qRT-PCR,indicated their importance in early defoliation.These results provide new insights into the potentially essential genes that respond to early defoliation in pear and lay a foundation for future research on the molecular breeding of resistant accessions.

    Acknowledgements

    This work was supported by the earmarked fund for Jiangsu Agricultural Industry Technology System,China (JATS[2021]453),the National Key Research and Development Program of China (2021YFD1200200),and the earmarked fund for China Agriculture Research System (CARS-28).We thank the State Key Laboratory of Crop Genetics and Germplasm Enhancement,Centre of Pear Engineering Technology Research,and the Bioinformatics Center of Nanjing Agricultural University,China for supporting this project.

    Declaration of competing interest

    The authors declare that they have no conflict of interest.

    Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm

    久久久色成人| 在线播放无遮挡| 久热久热在线精品观看| 国产成人a∨麻豆精品| 亚洲,一卡二卡三卡| 久久国产亚洲av麻豆专区| 亚洲色图综合在线观看| 亚洲精品久久午夜乱码| 深爱激情五月婷婷| 天美传媒精品一区二区| 少妇猛男粗大的猛烈进出视频| 91精品国产国语对白视频| 大片免费播放器 马上看| videossex国产| 久久国内精品自在自线图片| 久久久久网色| 麻豆精品久久久久久蜜桃| 亚洲成人一二三区av| 高清在线视频一区二区三区| 舔av片在线| 国产精品久久久久久久电影| 这个男人来自地球电影免费观看 | 18禁在线播放成人免费| av国产久精品久网站免费入址| 国产黄色视频一区二区在线观看| 熟女人妻精品中文字幕| 人妻少妇偷人精品九色| 少妇熟女欧美另类| 18禁动态无遮挡网站| 国产日韩欧美在线精品| 在线观看免费视频网站a站| 国产av一区二区精品久久 | 男女啪啪激烈高潮av片| 一区二区三区精品91| 哪个播放器可以免费观看大片| 国产一区二区在线观看日韩| 看非洲黑人一级黄片| 国产精品偷伦视频观看了| 欧美3d第一页| 精品一区二区三区视频在线| 午夜老司机福利剧场| 少妇精品久久久久久久| 日本一二三区视频观看| 成人影院久久| 亚洲精品色激情综合| 免费不卡的大黄色大毛片视频在线观看| 国产精品国产三级专区第一集| 国产免费一级a男人的天堂| 99视频精品全部免费 在线| 日韩成人伦理影院| 亚洲最大成人中文| 大陆偷拍与自拍| 2022亚洲国产成人精品| 街头女战士在线观看网站| 国产一区二区三区av在线| 五月玫瑰六月丁香| 欧美日韩视频精品一区| 成人毛片60女人毛片免费| 国产熟女欧美一区二区| 三级国产精品欧美在线观看| 久久久亚洲精品成人影院| 久久国产乱子免费精品| 国产成人免费观看mmmm| 一级a做视频免费观看| 欧美一区二区亚洲| 六月丁香七月| 美女cb高潮喷水在线观看| 成年免费大片在线观看| 欧美成人午夜免费资源| 丝袜喷水一区| 人妻夜夜爽99麻豆av| 国产一区亚洲一区在线观看| 亚洲av免费高清在线观看| 亚洲国产毛片av蜜桃av| 男女啪啪激烈高潮av片| 亚洲,一卡二卡三卡| 99热这里只有精品一区| 日本免费在线观看一区| 欧美国产精品一级二级三级 | 妹子高潮喷水视频| 人人妻人人爽人人添夜夜欢视频 | 一本—道久久a久久精品蜜桃钙片| 性色avwww在线观看| 涩涩av久久男人的天堂| 亚洲精品自拍成人| 深爱激情五月婷婷| 97在线人人人人妻| 永久免费av网站大全| 日韩av免费高清视频| 人妻少妇偷人精品九色| 成人高潮视频无遮挡免费网站| 黄片wwwwww| 国产免费视频播放在线视频| 国产 一区 欧美 日韩| 制服丝袜香蕉在线| 国产亚洲最大av| 韩国av在线不卡| 久久av网站| 日本-黄色视频高清免费观看| 亚洲自偷自拍三级| 麻豆乱淫一区二区| 欧美人与善性xxx| 免费观看a级毛片全部| 美女福利国产在线 | 哪个播放器可以免费观看大片| 久久6这里有精品| 国产无遮挡羞羞视频在线观看| 色综合色国产| 精品午夜福利在线看| 国产淫片久久久久久久久| 午夜免费观看性视频| 亚洲av国产av综合av卡| 亚洲欧美日韩卡通动漫| 中文字幕亚洲精品专区| 寂寞人妻少妇视频99o| 亚洲内射少妇av| 亚洲欧美一区二区三区国产| 少妇高潮的动态图| 国产亚洲欧美精品永久| 国产精品伦人一区二区| 国产在线免费精品| 啦啦啦在线观看免费高清www| a 毛片基地| 蜜桃在线观看..| 深夜a级毛片| 国产成人免费无遮挡视频| 日日啪夜夜爽| 大话2 男鬼变身卡| 日韩成人av中文字幕在线观看| 特大巨黑吊av在线直播| 免费观看在线日韩| 日本av手机在线免费观看| 亚洲精品乱码久久久v下载方式| 人体艺术视频欧美日本| 男女国产视频网站| 亚洲不卡免费看| 高清午夜精品一区二区三区| 舔av片在线| 纯流量卡能插随身wifi吗| av天堂中文字幕网| 男女国产视频网站| av黄色大香蕉| 亚洲真实伦在线观看| 嫩草影院新地址| 老熟女久久久| 婷婷色麻豆天堂久久| 黄色怎么调成土黄色| 内地一区二区视频在线| 91午夜精品亚洲一区二区三区| 久久99热这里只频精品6学生| 久久久久久久大尺度免费视频| 国产av精品麻豆| 国产欧美另类精品又又久久亚洲欧美| 欧美3d第一页| av在线蜜桃| 在线观看国产h片| 18禁动态无遮挡网站| 肉色欧美久久久久久久蜜桃| 国产精品女同一区二区软件| 国产日韩欧美在线精品| 欧美亚洲 丝袜 人妻 在线| 国产av国产精品国产| 久久精品国产鲁丝片午夜精品| 三级经典国产精品| 国产精品.久久久| 国产精品久久久久久精品电影小说 | 国产成人a区在线观看| 秋霞伦理黄片| 日本黄色日本黄色录像| 中文字幕av成人在线电影| 2018国产大陆天天弄谢| 91精品伊人久久大香线蕉| 中文字幕亚洲精品专区| 亚洲最大成人中文| 日韩av不卡免费在线播放| 国产成人91sexporn| 丰满乱子伦码专区| 嫩草影院入口| 国产精品蜜桃在线观看| 国产又色又爽无遮挡免| 最近中文字幕2019免费版| 国产精品成人在线| 免费av中文字幕在线| 亚洲精品国产av成人精品| 成年av动漫网址| 国产伦在线观看视频一区| 成人毛片a级毛片在线播放| 一本—道久久a久久精品蜜桃钙片| 久久韩国三级中文字幕| 少妇的逼水好多| 在线观看一区二区三区| 国产精品一区www在线观看| 99国产精品免费福利视频| 亚洲国产精品专区欧美| 香蕉精品网在线| 久久人妻熟女aⅴ| 免费大片18禁| 青春草视频在线免费观看| 亚洲美女黄色视频免费看| 亚洲va在线va天堂va国产| 欧美变态另类bdsm刘玥| 久久久久久伊人网av| 日韩一区二区三区影片| 日本黄大片高清| 日韩成人av中文字幕在线观看| 国产一区有黄有色的免费视频| 蜜臀久久99精品久久宅男| 亚洲高清免费不卡视频| 亚洲精品一二三| 菩萨蛮人人尽说江南好唐韦庄| 久久精品夜色国产| 久久女婷五月综合色啪小说| 日本av手机在线免费观看| 亚洲第一av免费看| 精品国产乱码久久久久久小说| 中文乱码字字幕精品一区二区三区| 日韩,欧美,国产一区二区三区| 久久久国产一区二区| 久久久久久久国产电影| 黄片无遮挡物在线观看| 成年人午夜在线观看视频| 欧美激情极品国产一区二区三区 | 亚洲精品一二三| 久久久久久伊人网av| 建设人人有责人人尽责人人享有的 | 欧美激情极品国产一区二区三区 | 日本av手机在线免费观看| 国产黄色视频一区二区在线观看| 欧美+日韩+精品| 少妇人妻一区二区三区视频| 国产综合精华液| 久久久久久久大尺度免费视频| 一级毛片黄色毛片免费观看视频| 中文精品一卡2卡3卡4更新| 日韩亚洲欧美综合| 高清毛片免费看| 国产成人午夜福利电影在线观看| av.在线天堂| 搡女人真爽免费视频火全软件| 观看美女的网站| 蜜桃亚洲精品一区二区三区| 久久亚洲国产成人精品v| 亚洲精品视频女| 亚洲精品国产色婷婷电影| 国产中年淑女户外野战色| 国产在线男女| 男女啪啪激烈高潮av片| 国产精品三级大全| 亚洲精品第二区| 日韩一区二区视频免费看| 最新中文字幕久久久久| 亚洲精品第二区| av.在线天堂| 欧美zozozo另类| 国产av一区二区精品久久 | 亚洲精品中文字幕在线视频 | 亚洲欧美精品专区久久| 亚洲精品成人av观看孕妇| 欧美区成人在线视频| 久久综合国产亚洲精品| 蜜桃在线观看..| 亚洲精品色激情综合| 狠狠精品人妻久久久久久综合| 亚洲av国产av综合av卡| 亚洲精品日本国产第一区| 亚洲中文av在线| 亚洲精品国产成人久久av| av在线老鸭窝| 十分钟在线观看高清视频www | av在线老鸭窝| 国产成人免费观看mmmm| 国产一级毛片在线| 亚洲av中文av极速乱| 一本久久精品| 久久女婷五月综合色啪小说| 亚洲精品久久午夜乱码| 亚洲久久久国产精品| 少妇的逼好多水| 中文字幕人妻熟人妻熟丝袜美| av播播在线观看一区| 人妻少妇偷人精品九色| 熟女电影av网| 性色avwww在线观看| 亚洲最大成人中文| 黄色怎么调成土黄色| 小蜜桃在线观看免费完整版高清| 国产亚洲欧美精品永久| 国产色婷婷99| 免费不卡的大黄色大毛片视频在线观看| 国产91av在线免费观看| 看非洲黑人一级黄片| 国产精品一二三区在线看| av在线蜜桃| 国语对白做爰xxxⅹ性视频网站| 成人高潮视频无遮挡免费网站| 搡老乐熟女国产| 男女无遮挡免费网站观看| 菩萨蛮人人尽说江南好唐韦庄| 国产v大片淫在线免费观看| 久久久亚洲精品成人影院| 久久国内精品自在自线图片| av福利片在线观看| 午夜福利视频精品| 欧美激情极品国产一区二区三区 | 熟女电影av网| 日韩电影二区| 日韩大片免费观看网站| 伊人久久国产一区二区| 亚洲综合色惰| 老师上课跳d突然被开到最大视频| 婷婷色综合大香蕉| 成人二区视频| 精华霜和精华液先用哪个| 青春草国产在线视频| 久久久亚洲精品成人影院| 亚洲精品日本国产第一区| 国产有黄有色有爽视频| 亚洲欧美清纯卡通| 人人妻人人澡人人爽人人夜夜| 一区二区三区四区激情视频| 伊人久久精品亚洲午夜| 国产av国产精品国产| 99re6热这里在线精品视频| 一级毛片电影观看| 中文字幕久久专区| www.色视频.com| 晚上一个人看的免费电影| 日产精品乱码卡一卡2卡三| 亚洲欧美一区二区三区国产| 久久热精品热| 99热6这里只有精品| 永久免费av网站大全| 午夜福利影视在线免费观看| 免费观看在线日韩| 日韩亚洲欧美综合| 交换朋友夫妻互换小说| 最近2019中文字幕mv第一页| 国产永久视频网站| 国产黄色视频一区二区在线观看| 国产爱豆传媒在线观看| 午夜激情福利司机影院| 免费人妻精品一区二区三区视频| 18禁在线无遮挡免费观看视频| 国产探花极品一区二区| 国产亚洲91精品色在线| 久久久久久久国产电影| 欧美变态另类bdsm刘玥| 国产成人a∨麻豆精品| 国产欧美日韩精品一区二区| 久久久久精品久久久久真实原创| 日韩欧美精品免费久久| 一个人看视频在线观看www免费| 久久精品久久久久久噜噜老黄| 亚洲,一卡二卡三卡| 亚洲丝袜综合中文字幕| 看非洲黑人一级黄片| 久久97久久精品| 国产成人一区二区在线| 中文精品一卡2卡3卡4更新| 久久精品国产a三级三级三级| 亚洲精品日韩av片在线观看| 女人十人毛片免费观看3o分钟| 久久久a久久爽久久v久久| 日韩一本色道免费dvd| 免费观看无遮挡的男女| 免费人成在线观看视频色| av免费观看日本| 女性被躁到高潮视频| www.av在线官网国产| 国国产精品蜜臀av免费| 丝瓜视频免费看黄片| 亚洲精品视频女| 97热精品久久久久久| 在线播放无遮挡| 各种免费的搞黄视频| 晚上一个人看的免费电影| 午夜福利网站1000一区二区三区| 久久6这里有精品| 男女无遮挡免费网站观看| 亚洲熟女精品中文字幕| 国内精品宾馆在线| 黄色一级大片看看| 免费久久久久久久精品成人欧美视频 | www.色视频.com| 国产淫片久久久久久久久| 国产探花极品一区二区| 高清日韩中文字幕在线| 国产 精品1| 99精国产麻豆久久婷婷| 国产综合精华液| 一级毛片aaaaaa免费看小| 亚洲av男天堂| 婷婷色av中文字幕| 国产爱豆传媒在线观看| 精品久久久久久久末码| 97超碰精品成人国产| 日韩中文字幕视频在线看片 | 99热国产这里只有精品6| 亚洲欧美一区二区三区国产| 欧美区成人在线视频| 嫩草影院新地址| 国产午夜精品一二区理论片| 精品午夜福利在线看| 欧美日韩亚洲高清精品| 国产成人91sexporn| av在线老鸭窝| 丰满乱子伦码专区| 免费看光身美女| 一区二区三区四区激情视频| 美女脱内裤让男人舔精品视频| 大香蕉97超碰在线| 亚洲一区二区三区欧美精品| 国产精品久久久久久久久免| 久久久久久久久大av| 国产日韩欧美在线精品| 精品一区在线观看国产| 亚洲久久久国产精品| 在线观看免费高清a一片| 成人亚洲欧美一区二区av| 一区二区三区乱码不卡18| 欧美三级亚洲精品| 亚洲欧美日韩另类电影网站 | 人妻少妇偷人精品九色| 中文欧美无线码| 亚洲精品国产av蜜桃| 亚洲av男天堂| 久热这里只有精品99| 狂野欧美激情性xxxx在线观看| 国产精品99久久久久久久久| 午夜福利视频精品| 男女免费视频国产| 国产精品国产三级国产专区5o| 日本午夜av视频| 亚洲国产精品国产精品| 亚洲av福利一区| 99久久精品一区二区三区| 国产精品爽爽va在线观看网站| 亚洲国产精品999| 自拍欧美九色日韩亚洲蝌蚪91 | 色视频www国产| 欧美区成人在线视频| 国产精品女同一区二区软件| 成人午夜精彩视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| 91精品国产九色| 久久精品人妻少妇| 日韩免费高清中文字幕av| 欧美日本视频| 高清av免费在线| 97在线人人人人妻| 欧美日韩精品成人综合77777| 建设人人有责人人尽责人人享有的 | 国内少妇人妻偷人精品xxx网站| 亚洲四区av| 欧美最新免费一区二区三区| 午夜免费鲁丝| 国产成人一区二区在线| 噜噜噜噜噜久久久久久91| av专区在线播放| 亚洲av福利一区| 一个人免费看片子| 草草在线视频免费看| 91精品国产九色| 欧美亚洲 丝袜 人妻 在线| 五月天丁香电影| av天堂中文字幕网| 最近中文字幕高清免费大全6| 国产色婷婷99| 免费观看无遮挡的男女| 麻豆成人av视频| 美女内射精品一级片tv| 91久久精品电影网| 啦啦啦啦在线视频资源| 亚洲国产av新网站| 精品视频人人做人人爽| 精品少妇久久久久久888优播| 亚洲经典国产精华液单| 丝袜脚勾引网站| 亚洲熟女精品中文字幕| 最近最新中文字幕大全电影3| 丰满乱子伦码专区| 国产白丝娇喘喷水9色精品| 亚洲精品国产av成人精品| 亚洲精品乱码久久久久久按摩| 国产无遮挡羞羞视频在线观看| 久久久久久久久大av| 99久久精品一区二区三区| 日本爱情动作片www.在线观看| 国产高清不卡午夜福利| 人妻少妇偷人精品九色| 黄片wwwwww| 免费少妇av软件| 美女xxoo啪啪120秒动态图| 女性生殖器流出的白浆| 亚洲国产精品成人久久小说| 五月伊人婷婷丁香| 夜夜爽夜夜爽视频| 久久99热这里只频精品6学生| 你懂的网址亚洲精品在线观看| 亚洲欧美一区二区三区国产| 最新中文字幕久久久久| 成年女人在线观看亚洲视频| 亚洲欧美日韩卡通动漫| 永久免费av网站大全| 蜜桃亚洲精品一区二区三区| av又黄又爽大尺度在线免费看| 蜜桃在线观看..| 亚洲高清免费不卡视频| 日韩av在线免费看完整版不卡| 亚洲精品国产av蜜桃| 国产亚洲91精品色在线| 日本黄色日本黄色录像| 亚洲精品日本国产第一区| 精品人妻视频免费看| 少妇精品久久久久久久| 在线观看美女被高潮喷水网站| 国产熟女欧美一区二区| 日本黄色片子视频| 久久久久久久久大av| 欧美高清成人免费视频www| 大片免费播放器 马上看| 在线观看av片永久免费下载| 欧美高清性xxxxhd video| 啦啦啦在线观看免费高清www| 内射极品少妇av片p| 久久午夜福利片| 久久久久久九九精品二区国产| 中文字幕亚洲精品专区| 一级a做视频免费观看| 亚洲中文av在线| 黄色视频在线播放观看不卡| 亚洲欧洲国产日韩| 久久这里有精品视频免费| 久久青草综合色| 午夜福利网站1000一区二区三区| 中文字幕亚洲精品专区| 精品久久国产蜜桃| 久久女婷五月综合色啪小说| 秋霞在线观看毛片| 中国三级夫妇交换| 久久精品国产亚洲av涩爱| 国产精品嫩草影院av在线观看| 日本色播在线视频| 成人一区二区视频在线观看| 91精品伊人久久大香线蕉| 亚洲精品亚洲一区二区| 97在线人人人人妻| 偷拍熟女少妇极品色| 九九在线视频观看精品| 99热这里只有是精品在线观看| 亚洲一级一片aⅴ在线观看| 欧美日本视频| 大香蕉97超碰在线| 久久精品国产自在天天线| 久久这里有精品视频免费| 99久久综合免费| 亚洲精品乱码久久久久久按摩| 伊人久久国产一区二区| 深爱激情五月婷婷| 三级国产精品欧美在线观看| 纯流量卡能插随身wifi吗| 免费av中文字幕在线| 黑丝袜美女国产一区| 色视频www国产| 亚洲av福利一区| 国产成人精品福利久久| 少妇 在线观看| 我的女老师完整版在线观看| 内射极品少妇av片p| 国产精品久久久久久精品电影小说 | 麻豆国产97在线/欧美| 国产av国产精品国产| 久久久久久久久久人人人人人人| 97热精品久久久久久| 亚洲国产高清在线一区二区三| 在线 av 中文字幕| 午夜日本视频在线| 久热久热在线精品观看| 少妇人妻久久综合中文| h视频一区二区三区| 午夜免费鲁丝| 大陆偷拍与自拍| 亚洲国产精品成人久久小说| 亚洲精品自拍成人| 干丝袜人妻中文字幕| 亚洲怡红院男人天堂| 亚洲av中文字字幕乱码综合| 天堂8中文在线网| 成人免费观看视频高清| 大香蕉97超碰在线| 日本与韩国留学比较| 熟妇人妻不卡中文字幕| 午夜激情久久久久久久| 国产精品不卡视频一区二区| 免费观看性生交大片5| 欧美最新免费一区二区三区| 美女cb高潮喷水在线观看| 亚洲av日韩在线播放| 欧美精品亚洲一区二区| 亚洲人成电影免费在线| av在线老鸭窝| 国产精品三级大全| 人妻一区二区av| 日日爽夜夜爽网站| 你懂的网址亚洲精品在线观看| 亚洲一区中文字幕在线| 天天躁夜夜躁狠狠久久av| 亚洲熟女毛片儿| 久9热在线精品视频| 新久久久久国产一级毛片| 亚洲精品国产一区二区精华液| 久久av网站| 99久久99久久久精品蜜桃| 男的添女的下面高潮视频| 999精品在线视频|