河北省是枣树主产区，主栽品种有赞皇大枣、婆枣、金丝小枣等。枣果成熟恰逢河北省秋季多雨时期，遇雨裂果严重，致使枣果丧失商品价值，效益低下，严重制约了枣产业的发展[1]。因此，对枣裂果机制进行研究，为枣裂果的防治和改善提供基础资料显得尤为迫切和重要。
导致裂果的因素复杂，品种、成熟时期、果实解剖结构、矿质元素含量和分布、果实表皮机械性、生长环境等均对裂果有一定的影响[2]。关于水分与裂果的关系，研究认为，裂果主要是可溶性糖含量下降、枣果不同组织水势不同以及果肉吸水率不同等造成的[3-4]。但是关于水分胁迫下裂果的分子机制还未见报道。
近年来，随着组学技术的发展，转录组越来越多的应用于基础生理调控的研究，目前，仅有少量对枣裂果相关基因的报道。对中秋酥脆枣枣果研究发现，AACT、PG、SHP、HXK1和LOX基因可能与枣裂果相关[5]；对自然裂果的伏脆蜜研究表明，裂果相关的基因可能有JAZ、ARF、SAUR和CML等[6]。然而基于2代平台的短读长RNAseq方法在识别复杂转录本isoform方面存在固有限制，因为它们不能测序全长转录本[7]。3代测序平台的全长转录本测序的出现打破了这一限制，例如Oxoford Nanopore Technology(ONT)平台全长转录组测序技术能够测得全长转录本，同时以较低的数据量即可实现定量，无须2代测序reads的辅助[8]。本研究以金丝新4号为试验材料，人工模拟降雨促使枣裂果，基于ONT平台对裂果前后的果实进行全长转录本的测序，从转录本水平阐明裂果相关机制，为完善枣基因组的注释和筛选裂果相关基因奠定基础。
1 材料和方法
1.1 试验材料与试验设计
以石家庄果树研究所内种植的5 a的金丝新4号品种为试验材料，在脆熟期选取生长一致的金丝新4号枣树3株。在每株东南西北4个方向外围摘取生长一致的果实20个。取样后以固定流速喷淋果树，每隔5 h观察裂果发生，至裂果发生初期(50 h)按照整树裂果比例(10%)再次选取果实20个。每次取样后削取果皮，液氮速冻，然后-80 ℃冰箱保存用于测定。每个样品采用20个枣果取样混合为1次重复，每处理重复3次。
1.2 RNA的提取、文库构建与测序
采用天根RNAprep Pure植物试剂盒提取样品中的总RNA，用琼脂糖凝胶电泳、Nanodrop、Agilent 2100、Qubit 等检测合格后，进行等量混匀。将mRNA反转录为cDNA，加上switch oligo PCR接头，合成互补链，进行14轮带LongAmp标签的PCR扩增。然后用T4 DNA连接酶(NEB)将PCR产物和ONT接头连接并用Agencourt XP磁珠进行DNA纯化，最后按照Ligation Sequencing Kit(SQK-LSK109)说明书进行文库的制备，在Nanopore PromethION平台进行测序。
1.3 全长转录组测序数据分析
1.3.1 数据过滤与比对 从原始下机序列中过滤低质量(长度小于500 bp，Q值小于7)序列和核糖体RNA序列，并根据序列两端是否存在引物得到全长序列。将获得的全长序列进行Polish得到一致性序列，然后对一致性序列用Minimap 2软件与冬枣基因组[9-10]或构建的Contig比对结果进行比对和去冗余。用整合基因组浏览器进行可视化，利用Gffcompare[11]进行基因结构的优化。
1.3.2 基因结构分析 通过Astalavista软件获取每个样品存在的可变剪接类型[12]；采用TAPIS pipeline来识别可变多聚腺苷酸化(APA)[13]，并利用MEME分析polyA剪切位点上游50 nt的顺式作用元件；从去冗余后的转录本中筛选500 bp以上的转录本，利用MISA软件做SSR分析[14]；应用CPC、CNCI、CPAT、pfam蛋白结构域分析的4种方法对新发现的转录本进行lncRNA的预测，并取其交集作为最终预测结果。
1.3.3 转录本表达定量 转录本表达定量采用CPM(Counts per million)作为衡量转录本或基因表达水平的指标，CPM计算公式如下：
CPM =(比对到某一转录本上的reads/比对到参考转录组的片段总数)× 1 000 000
1.3.4 差异转录本筛选 使用DESEQ对组间差异基因进行筛选，以差异倍数(FC)>2且错误发现率(FDR)<0.01作为筛选标准。
1.3.5 功能注释与富集 使用 Blast [15]软件将所有序列和差异表达序列与NR、Swiss-Prot、GO、COG、KOG、Pfam和KEGG数据库进行序列比对和功能注释。使用KOBAS[16]软件以KEGG数据库中通路为单位，应用超几何检验，找出与整个基因背景相比，在差异表达转录本中显著性富集的通路。
1.4 qRT-PCR验证
挑选裂果相关的10个转录本，设计引物后送至上海生工公司合成，将提取的RNA反转录成cDNA后，使用Roche的SYBR Green qPCR试剂盒在Bio-Rad iQ 5.0 qPCR仪上进行扩增。相对基因表达量使用2-ΔΔCt法计算；使用内参基因标准化。每个基因的所有测定均为一式三份在相同条件下进行。
2 结果与分析
2.1 测序数据统计
对喷淋前金丝新4号(J-0-1、J-0-2、J-0-3)和喷淋后样品(J-50-1、J-50-2、J-50-3)共6 个样品进行全长转录组测序，共获得17.47 Gb Clean data，平均每个样品测序达到2.92 Gb，每个样本数据产出如表1所示。与冬枣基因组进行序列比对的结果表明，平均比对效率达83.78%。所有一致性转录本序列通过软件minimap 2比对到参考基因组后进行去冗余分析，最终获得13 368条转录本序列。
表1 金丝新4号果皮转录组测序数据及与参考基因组的比对统计
Tab.1 Statistics of sequence output of peel transcriptome and alignment with the reference genome
样品名称Sample IDClean reads数Clean reads全长序列reads数Full-length reads全长序列比例/%Full-length percentage (FL)比对到基因组的reads数Mapped reads比对率/%Mapped rates J-0-13 106 0502 193 25370.611 840 97283.94J-0-23 244 9212 296 33470.771 929 00484.00J-0-33 268 9342 286 64869.951 907 93883.44J-50-13 308 0412 318 90470.101 943 73683.82J-50-23 147 1812 208 19870.161 846 01183.60J-50-32 334 9471 689 69772.371 417 41283.89
2.2 基因结构分析
2.2.1 可变剪接 从Astalavista软件分析结果中，对转录本发生5种可变剪接(Alternative Splicing，AS)事件：外显子跳跃、可变3′剪接位点、可变外显子、可变5′剪接位点、内含子保留进行统计，喷淋前J-0组和喷淋后J-50组中预测的可变剪接事件数量比例统计见图1。其中J-0组239个可变剪接事件，而J-50组有189个可变剪接事件。喷淋后可变3′剪接位点类型比例有所增加，内含子保留则有所减少。
J-0.喷淋前；J-50.喷淋后；5种可变剪接包括外显子跳跃、可变3′端剪接位点、可变外显子、可变5′端剪接位点、内含子保留。
J-0.Peels without spraying；J-50.Peels with spraying；The alternative splicing includes 5 types：exon skipping，alternative 3′ slice site，mutually exclusive exon，alternative 5′ slice site and intron retention.
图1 喷淋前后果皮可变剪接事件统计
Fig.1 Statistics of alternative splicing in groups with and without spraying
2.2.2 可变多聚腺苷酸化 利用3代Iso-Seq技术检测与分析可变多聚腺苷酸化(Alternative Polyadenylation，APA)，结果如图2所示。APA事件中鉴定到9 178个PolyA位点，其中6 743个基因有2个或2个以上PolyA位点，占表达基因的70.48%；J-0组与J-50组PolyA位点较接近，J-50组有2个或2个以上PolyA位点的基因有所减少。利用MEME共找到显著富集的PolyA信号(PAS)motif 78个。
2.2.3 简单重复序列 用MISA从11 100条500 bp以上的转录本中共鉴定到7 321个简单重复序列(Simple Single Repeat，SSR)，各种类型的SSR的数目如表2所示。其中包含1个以上SSR的转录本有1 713个。
J-0.喷淋前；J-50.喷淋后。
J-0.Peels without spraying；J-50.Peels with spraying.
图2 可变多聚腺苷酸化位点统计
Fig.2 Statistics of alternative polyadenylation
表2 SSR分析结果统计
Tab.2 Statistics of SSR
类型Searching items数量Numbers序列总数Total number of sequences examined11 100包含SSR的序列数Number of SSR containing sequences4 799包含1个以上SSR序列数Number of sequences containing more than 1 SSR1 713单碱基Mono nucleotide3 433双碱基Di nucleotide2 317三碱基Tri nucleotide1 449四碱基Tetra nucleotide67五碱基Penta nucleotide23六碱基Hexa nucleotide32
2.2.4 长链非编码RNA预测 因长链非编码RNA(long non-coding RNA，lncRNA)不编码蛋白，通过对转录本进行编码潜能筛选，判断其是否具有编码潜能，从而可以判定该转录本是否为lncRNA。用4种分析方法预测lncRNA得到的结果取交集，结果如图3所示，共预测到114个lncRNA。
图3 根据4种筛选方法预测lncRNA的维恩图
Fig.3 Venn diagram of predicted lncRNA by four kinds of methods
2.3 转录本/基因表达定量和差异转录本/基因筛选分析
转录组测序J-0组和J-50组分别测得21 350，19 298个转录本，对转录本进行合并后J-0组和J-50组分别得到14 873，13 795个表达基因。其中2组共有转录本18 130个，共有基因13 098个。对转录本的表达水平CPM绘图结果如图4所示。
图4 样品转录本CPM箱线图
Fig.4 Boxplot of the transcript CPM of different samples
如图5所示，以FC>2和FDR<0.01作为筛选标准共筛选到248个差异转录本，其中179个上调，69个下调；对筛选到的差异转录本进行注释，注释结果如表3所示。
每个点代表一个转录本；横坐标.基因表达量变化倍数取2的对数；纵坐标.错误发现率取10的对数的负数；红色.显著上调的转录本；绿色.显著下调的转录本；黑色.无显著变化的正常转录本。
Every dot. One transcript；X-axis. The log2(fold change)of the gene expression；Y-axis. The lg(false discovery rate)；Red dots. The significantly up-regulated transcripts；Green dots. The significantly down-regulated transcripts；Black dots. The transcripts neither up-regulated nor down-regulated.
图5 差异表达转录本火山图
Fig.5 Volcano plot of differentially exoressed transcripts
分别对上调和下调转录本进行GO富集，选取生物学过程、细胞组分和分子功能中显著富集前3的条目(表4，5)，其中上下调转录本富集到的条目较为相似，都包括翻译、有丝分裂中心体分离正向调节、细胞质、催化前剪接体、胞质大核糖体亚基、核糖体结构成分、3-氧代酰基-[酰基载体蛋白]合酶活性等。上调的转录本还富集到了木葡聚糖代谢过程、木葡聚糖：木葡糖基转移酶活性等GO条目，下调转录本还富集到了冷适应、葡聚糖内切1，3-β-D-葡萄糖苷酶活性等GO条目。
对上下调转录本进行KEGG富集，结果如图6所示，上调转录本显著富集(q<0.05)的通路包括内质网蛋白质加工、吞噬体、氨基糖和核糖代谢以及淀粉和蔗糖代谢通路；下调转录本显著富集的通路包括内质网蛋白加工、氮代谢和光合作用-天线蛋白通路。
表3 注释的差异表达转录本或基因数量统计
Tab.3 Statistics of differentially expressed transcripts with annotation
数据库Database总数TotalCOGGOKEGGKOGNRPfamSwiss-ProteggNOG差异转录本DETs234105172102102233190202228
表4 上调表达转录本GO富集
Tab.4 GO enrichment of up-regulated transcripts
类别DomainGO代号GO_ID描述TermKS值KS生物学过程GO:0006412翻译4.90E-10BiologicalGO:0010411木葡聚糖代谢过程4.00E-07processGO:0046604有丝分裂中心体分离正向调节1.00E-06细胞组分GO:0005829细胞质3.10E-07Cellular GO:0071011催化前剪接体4.50E-07componentGO:0022625胞质大核糖体亚基3.00E-06分子功能GO:0003735核糖体结构成分4.30E-13Molecular GO:0016762木葡聚糖:木葡糖基转移酶活性4.30E-07functionGO:00043153-氧代酰基-[酰基载体蛋白]合酶活性2.30E-06
表5 下调表达转录本GO富集
Tab.5 GO enrichment of down-regulated transcripts
类别DomainGO代号GO_ID描述TermKS值KS生物学过程GO:0006412翻译1.20E-09BiologicalGO:0046604有丝分裂中心体分离正向调节9.20E-07processGO:0009631冷适应3.10E-06细胞组分GO:0071011催化前剪接体3.90E-07Cellular componentGO:0005829细胞质5.90E-07GO:0022625胞质大核糖体亚基2.40E-06分子功能GO:0003735核糖体结构成分2.40E-13MolecularGO:00043153-氧代酰基-[酰基载体蛋白]合酶活性2.20E-06functionGO:0042973葡聚糖内切1,3-β-D-葡萄糖苷酶活性6.70E-06
纵坐标.富集到的通路；富集因子.转录本富集程度；q值.富集的显著性。
Y-axis.The enriched pathways；Rich factor.The degree of transcripts enrichment；q value.The significance of enrichment.
图6 上调(A)和下调(B)转录本KEGG富集散点图
Fig.6 Scatter plot of KEGG enrichment of up-(A)and down-regulated(B)transcripts
进一步对富集的通路分析发现，内质网蛋白加工和吞噬通路中，编码钙网蛋白CNX(rna5911)和钙联蛋白CRT(rna4420.4，rna4420.7，ONT.19253)的转录本显著上调；内质网蛋白加工通路中热激蛋白转录本HSPA5(ONT.5890.1，ONT.5890.2)、HSP70(ONT.1837.2，ONT.1837.3，ONT.1837.4，ONT.6218.2，rna7964，rna7970)及部分HSP20(rna8760，rna8845)转录本显著上调，HSP20(ONT.1142.3，ONT.3347.1，rna15259，rna20226，rna35213，rna5102，rna9546)转录本显著下调；吞噬体通路中的微管蛋白αTUBA转录本(ONT.7131.2，rna33179)和微管蛋白β转录本TUBB(rna9061)显著上调。氨基糖和核糖代谢中显著上调的转录本还包括果糖激酶scrK(ONT.3895.1)、α-1，4-半乳糖醛酸转移酶GAUT(ONT.728.2，rna40026)、UDP葡萄糖6-脱氢酶UGD6(rna10718，rna10719)、UDP-葡萄糖醛酸4-表异构酶(ONT.3958.2)、UDP-芹菜糖/木糖合酶AXS(ONT.1060.1)等。根据富集结果构建的蛋白质网络互作图，如图7所示。
每个圆.一个蛋白；连线.相关作用关系；圆的大小.蛋白的连通性；圆的颜色.聚集系数；红色越多聚集系数越高。
Every circle.One protein；The line.Interactions between proteins；The size of circle.The connectivity of the protein；The color of circle.The clustering coefficient of the protein；Red indicates higher clustering coefficient than green.
图7 差异基因编码蛋白网络互作图
Fig.7 Protein-protein interaction of differentially expressed transcripts
2.4 qRT-PCR验证
选取10个候选转录本进行实时荧光定量对其表达量进行验证，引物序列如表6所示，结果如图8所示，10个转录本的表达水平和变化趋势与3代转录组的结果相一致，证明了3代测序结果的准确性和可靠度。
表6 实时荧光定量引物
Tab.6 qRT-PCR primers
转录本IDTranscript ID转录本名称Transcript description正向引物(5'-3')Forward primer反向引物(5'-3')Reverse primerrna5911CALRGTGAAGACAGCAAGACTACCTAACGATGAACCAGAAONT.5890.1BIPCCGAGAAGAGGACTACGATCTGCTGATTCACCACTGONT.4220.4CANXCGATAATGAGTTGCGAATGAATTGTCTTGGAAGGAATONT.1837.2HSP70CACCATTACCAACGATAAGCCTCATCTTCTGACTTGTrna8845HSP22CGTGTTCCTACTGTTGATGTCTCCGACTCCTCTCTTGATAONT.7131.2TUBAGAAGGATTACGAAGAAGTTGGAAATAGCAAGTCATTrna9061TUBBCTCTTACATCTCGTGGTTCATCATATTCTTGGCATCCONT.728.2GAUTTTAGGACTCATCGTGGTAGGAAGGTGGTTCAGGATAGAONT.1142.3HSP18.5GAGGAAATTCAGGCTACCAACAGTAACGGTGAGAACrna11228SCRKCAGCCATGAAGAGGTAGGTGAATCTAAGTAATGACAG-CAGCGAGTT
图8 候选转录本qRT-PCR验证
Fig.8 qRT-PCR validation of candidate transcripts
3 讨论与结论
本研究利用ONT平台对喷淋前后的金丝新4号进行测序，共得到17.47 Gb Clean data，与基因组平均比对效率达83.78%，证明了测序数据的准确性和有效性。通过数据过滤与去冗余，6个样本共得到13 368条转录本。由于读长长，无须拼接组装，3代ONT测序在定量的同时还能够对转录本的结构进行分析[8]。本研究对转录本中存在的AS、APA、SSR进行了统计，并对lncRNA进行预测，对照组和处理组分别鉴定到239，189个可变剪接事件，APA事件中鉴定到9 178个PolyA位点，并找到78个显著富集的PAS信号，预测到7 321个SSR和114个lncRNA，完善了基因组的注释。
可变剪接是真核生物基因表达的普遍调节机制，mRNA通过可变剪接对mRNA的稳定性、定位或翻译进行调控，增加了基因的多样性，增强了生物体对环境的适应性[17]。本研究结果表明，喷淋后裂果的J-50组转录本和可变剪接数目较J-0组减少，这与之前报道的棉花在盐胁迫下可变剪接数目增加的情况并不相同[18]，这可能是由于裂果对果实造成了不可逆的损伤，导致表达转录本和可变剪接的数量变少。从可变剪接的各类型所占比例来看，本研究裂果后可变外显子类型比例增加，内含子保留比例下降，与之前报道的玉米穗部低氮胁迫下内含子保留比例上升也不一致[19]，可能是由于物种的差别或是裂果的不可逆性。
钙网蛋白(Calnexin，CNX)和钙联蛋白(Calreticulin，CRT)是同源蛋白，具有相似的结构和功能。二者都是内质网中的钙结合蛋白，具有分子伴侣的生物学功能，能够调节Ca2+的稳态、细胞黏附和基因的表达[20]。Kim等[21]报道了CRT与拟南芥对水分胁迫的敏感性有关；此前关于番茄和文旦柚的研究中均报道了CNX基因和蛋白的差异表达与裂果有关[22-23]。本研究中CNX和CRT均显著上调，与上述结果相一致。通常认为，裂果的发生与水分胁迫下果皮与果肉存在水势和渗透压差别，果实吸水膨胀超过果皮承受能力而发生裂果[24]。在水分胁迫下细胞吸水，细胞中Ca2+浓度失衡，CRT和CNX上调调节Ca2+浓度及Ca2+参与的相关通路。同时注意到，本研究中分子伴侣BIP(HSPA5)转录本也显著上调表达，CRT 与BIP形成CRT-BIP复合体，促进了蛋白折叠[25]。
Koumanov[26]认为裂果与长时间降雨果皮温度下降造成冷应激，果皮皱缩开裂有关。冷应激则与热激蛋白等密切相关。HSP70主要参与蛋白的折叠和运输，并阻止蛋白聚合，阻止胁迫下的细胞程序性死亡，增强植物对非生物学胁迫的抗性[27]。小热激蛋白sHSP通过非极性方式与异常蛋白结合，在冷应激条件下可提高植物抗性[28]。本研究中，由于喷淋造成枣果冷应激，热激蛋白HSP70和部分sHSP转录本显著上调，GO富集也富集到了冷适应相关的条目(GO：0009631)，为Koumanov[25]的理论提供了一定的证据。
微管蛋白α 和微管蛋白β是细胞骨架的重要组成部分，对于维持细胞形态，参与物质运输等具有重要作用。本研究中，微管蛋白转录本显著上调，与之前报道的易裂品种与抗裂品种枣中相比，微管蛋白显著上调表达相一致[21]。
之前的研究中发现，裂果的发生与胞壁结构和果皮物质以及相应的细胞壁代谢酶等密切相关[6，29]。在本研究中，果皮中果糖激酶scrK、UDP葡萄糖6-脱氢酶UGD6、α-1，4-半乳糖醛酸转移酶GAUT、UDP-葡萄糖醛酸4-表异构酶、UDP-芹菜糖/木糖合酶AXS等与细胞壁组分相关的酶的转录本显著上调。其中，UGD6和UDP-葡萄糖醛酸4-表异构酶等在形成细胞壁多糖及纤维素、半纤维素前体过程中起着非常重要的作用[30]，GUAT是果胶合成中研究最多和最关键的酶[31]。由此可见，水分胁迫下，枣裂果后细胞壁合成相关的转录本显著上调，以应对裂果的发生。
ONT 3代全长转录组克服了2代转录组的缺点，以较低的数据量完成了全长转录本的测序和定量，完善了基因组的注释。水分胁迫下的裂果机制可能包括可变剪接事件的减少、钙网蛋白和钙调蛋白对钙离子的调节、持续喷淋后冷应激造成的热激蛋白的上调及细胞壁组分合成相关转录本的上调等。本研究首次从转录本水平对水分胁迫下枣裂果相关的机制进行了研究，完善了枣基因组的注释，并根据研究结果构建了蛋白互作网络，未来的研究可以针对筛选出的基因通过过表达或敲除试验验证其功能，通过酵母杂交试验验证其存在的互作关系，为改善水分胁迫下的裂果提供理论基础。
[1] 宗晓鸿.连阴雨对沾化冬枣的不利影响及应对措施[J].落叶果树，2019，51(3)：63-64.doi：10.13855/j.cnki.lygs.2019.03.022.
Zong X H.Impact of continuous rainy disaster on Zhanhua winter jujube and countermeasures[J].Deciduous Fruits, 2019, 51(3): 63-64.
[2] 马雯彦，庞晓明，续九如，李颖岳.果实裂果影响因子研究进展[J].华中农业大学学报，2010，29(6)：798-804.
Ma W X, Pang X M, Xu J R, Li Y Y. Advances in research on the factors influencing fruit cracking[J]. Journal of Huazhong Aricultural University, 2010，29(6):798-804.
[3] 卢艳清．枣抗裂种质筛选及其抗裂机理初步研究[D].保定：河北农业大学，2008.doi：10.7666/d.y1306932.
Lu Y Q.Screening of germplasm with high risistance to fruit cracking and primary study on their resistant mechanism in Chinese jujube[D].Baoding：Hebei Agricultural University，2008.
[4] 李克志，高中山.枣裂果机理的初步研究[J].果树学报，1990，7(4)：221-226.
Li K Z，Gao Z S.Preliminary study on the mechanism of fruit cracking in Chinese jujube[J].Journal of Fruit Science，1990，7(4)：221-226.
[5] 魏琦琦.基于转录组测序的枣成花和果实品质重要基因的研究[D].长沙：中南林业科技大学，2017.doi：10.7666/d.Y3269462.
Wei Q Q.The study of key genes on flowering and fruit quality from Ziziphus jujuba based on transcriptome sequencing[D].Changsha：Central South University of Forestry and Technology,2017.
[6] 王建宇.枣裂果生理特征及转录组和蛋白组学分析[D].阿拉尔：塔里木大学，2020.doi：10.27708/d.cnki.gtlmd.2020.000008.
Wang J Y.Physiological characteristics，transcriptome and proteomic analysis of jujube fruit cracking[D].Alaer：Tarim University，2020.
[7] Sharon D，Tilgner H，Grubert F，Snyder M.A single-molecule long-read survey of the human transcriptome[J].Nature Biotechnology，2013，31：1009-1014.doi：10.1038/nbt.2705.
[8] Sessegolo C, Cruaud C, Da Silva C, Cologne A, Dubarry M, Derrien T, Lacroix V, Aury J M. Transcriptome profiling of mouse samples using nanopore sequencing of cDNA and RNA molecules[J].Scientific Reports，2019，9(1)：1-12.doi：10.1038/s41598-019-51470-9.
[9] Li H.Minimap2：pairwise alignment for nucleotide sequences[J].Bioinformatics，2018，34(18)：3094-3100.doi：10.1093/bioinformatics/bty191.
[10] Liu M J，Zhao J，Cai Q L，Liu G C，Wang J R，Zhao Z H，Liu P，Dai L，Yan G J，Wang W J，Li X S，Chen Y，Sun Y D，Liu Z G，Lin M J，Xiao J，Chen Y Y，Li X F，Wu B，Ma Y，Jian J B，Yang W，Yuan Z，Sun X C，Wei Y L，Yu L L，Zhang C，Liao S G，He R J，Guang X M，Wang Z，Zhang Y Y，Luo L H.The complex jujube genome provides insights into fruit tree biology[J].Nature Communications，2014，5(1)：1-12.doi：10.1038/ncomms6315.
[11] Pertea G，Pertea M.GFF Utilities：GffRead and GffCompare[J].F1000 Research，2020，9：304.doi：10.12688/f1000research.23297.1.
[12] Song Q A，Catlin N S，Brad B W，Li S.Computational analysis of alternative splicing in plant genomes[J].Gene，2019，685：186-195.doi：10.1016/j.gene.2018.10.026.
[13] Abdel-ghany S E, Hamilton M, Jacobi J L, Ngam P,DevittN,Schilkey F,Ben-Hur A, Reddy. A survey of the sorghum transcriptome using single-molecule long reads[J].Nature Communications,2016,7:11706.doi:10.1038/ncomms11706.
[14] Beier S, Thiel T, Münch T, Scholz U, Mascher M. MISA-web：a web server for microsatellite prediction[J].Bioinformatics，2017，33(16)：2583-2585.doi: 10.1093/bioinformatics/btx198.
[15] Hung J H，Weng Z.Sequence alignment and homology search with Blast and ClustalW[J].Cold Spring Harbor Protocols，2016，2016(11)：1016-1021.doi：10.1101/pdb.prot093088.
[16] Kanehisa M，Sato Y，Kawashima M，Furumichi M，Tanabe M.KEGG as a reference resource for gene and protein annotation[J].Nucleic Acids Research，2016，44(D1)：457-462.doi：10.1093/nar/gkv1070.
[17] Francisco E B，Jimena G.Alternative splicing as a regulator of development and tissue identity[J].Nature Review of Molecular Cell Biology，2017，18(7)：437-451.doi：10.1038/nrm.2017.27.
[18] Zhu G，Li W，Zhang F，Guo W.RNA-seq analysis reveals alternative splicing under salt stress in cotton，Gossypium davidsonii[J].BMC Genomics，2018，19(1)：73.doi：10.1186/s12864-018-4449-8.
[19] 杨海龙.玉米穗部低氮响应可变剪接的鉴定[D].北京：中国农业科学院，2015.
Yang H L.Characterization of alternative splicing in response to low nitrogen in maize ears[D].Beijing：Chinese Academy of Agricultural Sciences，2015.
[20] Jia X Y，He L H，Jing R L，Li R Z.Calreticulin：conserved protein and diverse functions in plants[J].Physiologia Plantarum，2009，136(2)：127-138.doi：10.1111/j.1399-3054.2009.01223.x.
[21] Kim J H, Nguyen N H, Nguyen N T, Hong, S W, Lee H.Loss of all three calreticulins，CRT1，CRT2 and CRT3，causes enhanced sensitivity to water stress in Arabidopsis[J].Plant Cell Reports，2013，32(12)：1843-1853.doi：10.1007/s00299-013-1497-z.
[22] Wi
syk A，Lirski M，Fogtman A，Zagórski-Ostoja W, Góra-Sochacka A. Differences in gene expression profiles at the early stage of Solanum lycopersicum infection with mild and severe variants of potato spindle tuber viroid[J].Virus Research，2020，286：198090.doi: 10.1016/j.virusres.2020.198090.
[23] 李淑一.不同授粉的度尾文旦柚裂果相关差异蛋白的研究[D].福州：福建农林大学，2009.doi：10.7666/d.y1515084.
Li S Y.Study on the differential proteins of fruit dehiscent in different pollinated Duweiwendan pomelo[D].Fuzhou：Fujian Agriculture and Forestry University，2009.
[24] Measham P F, Bound S A, Gracie A J, Wilson S J. Incidence and type of cracking in sweet cherry(Prunus avium L.)are affected by genotype and season[J].Crop and Pasture Science，2009，60(10)：1002-1008.doi：10.1071/CP08410.
[25] Crofts A J，Leborgne-Castel N，Pesca M，VitaleA，Denecke J.BiP and calreticulin form an abundant complex that is independent of endoplasmic reticulum stress[J].The Plant Cell，1998，10(5)：813-823.doi：10.2307/3870667.
[26] Koumanov K S.On the mechanisms of the sweet cherry(Prunus avium L.)fruit cracking：Swelling or shrinking?[J].Scientia Horticulturae，2015，184：169-170.doi：10.1016/j.scienta.2015.01.002.
[27] Song A P, Zhu X R, Chen F D, Gao H S, Jiang J F, Chen S M.A Chrysanthemum heat shock protein confers tolerance to abiotic stress[J].International Journal Molecular Science，2014，15(3)：5063-5078.doi：10.3390/ijms15035063.
[28] Maimbo M, Ohnishi K, Hikichi Y, Yoshioka H, Kiba A. Induction of a small heat shock protein and its functional roles in Nicotiana plants in the defense response against Ralstonia solanacearum[J].Plant Physiology，2007，145(4)：1588-1599.doi：10.1104/pp.107.105353.
[29] 李伟明, 陈晶晶, 段雅婕, 胡会刚, 庞振才, 胡玉林.番荔枝果实后熟过程多糖代谢与果实软化和采后裂果的关系[J].植物生理学报，2018，54(11)：1727-1736.doi：10.13592/j.cnki.ppj.2018.0154
Li W M，Chen J J，Duan Y J，Hu H G，Pang Z C，Hu Y L.Effects of anti-ethylene treatments on polysaccharide metabolism，fruit softening and their relationship with postharvest cracking of atemoya fruit[J].Plant Physiology Journal，2018，54(11)：1727-1736.
[30] Han J，Pan Y X，Wang X F，Zhang Y，Ma Z Y.Antisense expression of Gossypium barbadense UGD6 in Arabidopsis thaliana significantly alters cell wall composition[J].Science China Life Sciences，2016，59：213-218.doi：10.1007/s11427-016-5004-y.
[31] 陈凯莉，许轲，张贤聪，王亚楠，汪志辉，王迅.果实中果胶代谢相关酶基因的研究进展[J].园艺学报，2017，44(10)：2008-2014.doi：10.16420/j.issn.0513-353x.2016-0846.
Chen K L, Xu K, Zhang X C, Wang Y N, Wang Z H, Wang X.Advances in genes information involved in pectin metabolism in fruit[J]. Acta Hoticulturae Sinica, 2017, 44(10): 2008-2014.