1. 引言
橡胶草(Taraxacum kok-saghyz L. Rodin, TKS)也称为俄罗斯蒲公英，属于菊科植物，蒲公英属，橡胶草种，原产于新疆天山边界附近的特克斯河盆地与哈萨克斯坦的天山河谷接壤处 [1]。染色体数为2N = 2X = 16条 [2] [3]，二倍体的多年生草本植物 [4]。橡胶草的基因组大小约1.29 Gb，共包含46,000多个基因编码蛋白质 [5]。TKS的根中含有天然橡胶(Natural rubber, NR)，其中NR占总干重的3%~28% (多年生) [6] [7]，并且橡胶分子大小和性能与巴西橡胶树(Hevea brasiliensis)中的天然橡胶极为相似 [8] [9]。由于近缘物种短角蒲公英Taraxacum brevicorniculatum (TB)与TKS在生境和形态上较为类似且也能产天然橡胶 [10] [11]。TKS有着能够在寒冷和温带地区广泛生长等特点，生命期相对较短，含胶量和质量高，易转基因和收获方法较简单等优点 [12] [13]，橡胶草已被作为研究天然橡胶生物合成的理想模式植物，也正在成为天然橡胶替代生产最有发展潜力的作物。TKS根部除了含NR外，还含有约28%的次生代谢产物菊粉，菊粉是食品工业和生物乙醇制造中重要原料 [12] [13]，其副产品产值更大。然而，橡胶草在成为商业上可行的作物之前，仍然需要在橡胶产量和农艺性状等方面的改善。
早期由于橡胶草基础数据的缺乏，以基因组学、转录组学和蛋白组学等研究方法，挖掘橡胶草农艺性状相关的关键基因越来越受到关注 [14]。Luo等利用转录组挖掘橡胶产量相关的SNPs标记，为解析橡胶合成相关遗传性状和辅助选择育种基础 [15]。Cao等利用MeJA诱导橡胶草根的转录组分析发现，与橡胶形成相关的TFs与NR合成基因之间相互作用，共同调节次生代谢产物并影响橡胶的生物合成 [16]。Nowicki等利用基因组微卫星基因座的研究，gSSR和eSSR标记结果表明，TKS具有低至中等的遗传多样性，而该物种缺乏种群结构，这为菊科相关生物的多样性分析提供了一种具有分子标记功能性的工具，揭示了TKS遗传多样性与菊科物种密切相关 [17]。对TKS根部进行丙酮萃取，根据其提取物中三萜类化合物的组成分析，发现氧化角鲨烯环化酶TkOSC1和TkCYP716A263在橡胶草的乳胶中过表达后，会产生大量被修饰的五环三萜类化合物 [18]。橡胶草乳胶中敲除OSC基因，可以减少根中天然橡胶中三萜的含量。而在RNAi-OSC植株中并没有表现出明显的表型变化，含胶量也未受影响。这为TKS育种计划以及开发具有稳定低三萜烯含量的橡胶草植物提供坚实的基础 [19]。在橡胶草中使用CRISPR/Cas9敲除TkRALFL1基因后的植株根系生物量更高，菊粉和天然橡胶的产量也更高 [20]。橡胶草乳胶中敲除TkCPTL1基因，致使合成顺式-1,4-异戊二烯(IPP)减少，但TKS根中三萜和菊粉的含量却升高了，而三醇的含量未受影响。对TkCPTL1-RNAi植物的胶乳进行分析发现，具有与橡胶粒子相似的粒子将三萜烯包裹在由小橡胶粒子蛋白SRPP和稳定的磷脂壳中。相反，TkCPTL2基因表达下调，没有改变橡胶草的表型，这表明TkCPTL2蛋白在橡胶草中的功能有待进一步探究 [21]。TkCPTL1-RNAi植物和野生型橡胶草乳胶蛋白质组的比较发现了假定蛋白，这些蛋白与类异戊二烯生物合成途径，代谢产物以及合成橡胶的橡胶粒子膜的成分有关 [21]。Xie等对橡胶草成熟根的全蛋白鉴定建立了高质量的数据集，并提出橡胶延伸因子(REF)和SRPP在橡胶合成过程中同等重要 [22]。
Benninghausl等利用蛋白质组学和代谢组学的方法，在产胶和非产胶的橡胶草根中，Branase诱导乳胶RNA降解，对TKS根部代谢产物和蛋白质组成的影响，并鉴定出多种乳胶特异性蛋白，还包括迄今为止尚未在植物中进行讨论的类脂蛋白。Barnase的表达使有胶橡胶草显出非乳胶的表型，此结果有望作为植物与环境进一步相互作用过程中，对植物产乳胶进行前瞻性分析最有价值的基础，以阐明整个植物界中乳胶的进化和具体分布。此外，不同代谢途径的分子相互作用，也反应了橡胶合成网络关系有了更深入的了解 [23]。Zhang等利用间接超声结合尺寸排阻色谱法改进橡胶草橡胶提取结果显示，比常规提取方法所需的提取时间更短，效率更高，且该方法不影响TKS的橡胶分子量分布 [24]。Panara等比较非胶、低含胶量(Low Rubber, LR)和高含胶量(High Rubber, HR)三种橡胶草转录组差异基因表达，发现HR胶乳中橡胶含量较高，是由于许多直接参与橡胶合成基因表达呈正相关的结果，并表明天然橡胶的产生在转录水平上受到高度调控。而LR中较低胶含量可能与天然橡胶生物合成竞争的其它次生代谢产物合成中涉及的基因的较高表达有关 [25]。
橡胶草的叶片内含有一定量的花青素和菊糖，叶片的副产物在饲料、食品营养和药理等领域有着潜在的应用价值 [8] [26]。并且橡胶草叶片中的菊苣酸有潜在的抗氧化作用，可作为酚类化合物的抗氧化来源，也是中草药和生产健康食品/饲料的营养成分 [27]。综上所述发现，目前关于新疆橡胶草幼嫩与成熟叶片转录组的比较分析研究的鲜有报道，在本研究中，对新疆橡胶草幼叶及成熟叶片的转录组测序后的序列进行了组装、Unigene功能注释、参考基因组、差异表达基因分析、并在功能上注释差异表达基因分类和次级代谢产物合成分析等研究。橡胶草的转录数据也扩大并丰富了可利用的橡胶草基因组资源，并可能加速橡胶草的种质资源和育种改良。
2. 材料与方法
2.1. 材料与试剂
新疆橡胶草幼嫩叶片与成熟叶片取自石河子大学生命科学学院，新疆植物医学资源与利用教育部重点实验室种质资源圃内种植的新疆橡胶草。分别取种植30 d后的幼苗期(leaf_1)和种植120 d后的成熟期(leaf_2)长势良好，将健康的橡胶草叶片组织，迅速置于液氮中并保存至−80℃冰箱，以备试验所用。
2.2. 橡胶草叶片总RNA提取
使用改良的CTAB法提取橡胶草叶片总RNA，在CTAB方法的基础上，加1%的PVPP，加入β-巯基乙醇的SL 500 μL，两次氯仿抽提，并用RNA级的酚:氯仿:异戊醇 = 25:24:1进一步抽提，最后晾干沉淀，滴加RNase-Free ddH2O 30 μL，得到提取好的RNA溶液，标记备份−80℃保存。
2.3. 橡胶草叶片转录组测序
以新疆橡胶草幼苗期与成熟期种植苗的叶片为材料，分别提取叶片总RNA后，使用Nanodrop (Agilent 2100 Bioanalyzer)检测总RNA的浓度、总量、OD260/280、OD260/230和RIN值，RNA纯度使用NanoDropTM的紫外分光光度计检测，结果为A级。将合格总RNA片段化后，逆转录为cDNA，在cDNA的5'-和-3'末端加上接头，再进行PCR扩增，将扩增产物变性环化，最终获得单链环状的DNA文库。上机BGISEQ-500进行转录组测序，测序于华大基因公司(深圳)完成。
2.4. 原始测序数据过滤
测序的原始数据包含接头污染、低质量以及未知碱基N含量(>5%)过高的reads，数据分析之前需要去除这些reads后得到clean reads，以保证结果的可靠性。得到clean reads之后，使用HISAT将clean reads比对到参考基因组序列，使用Bowtie 2将clean reads比对到参考序列，之后再使用RSEM计算基因和转录本的表达水平。
2.5. 差异表达基因检测及生物信息学分析
根据样品组之间的差异表达基因(DEG)，使用DEseq 2和PossionDis方法进行差异基因检测，为了更直观地展示每个样品在不同FPKM区间的基因数目，对FPKM (FPKM ≤ 1、1 < FPKM < 10、FPKM ≥ 10)的三种情况进行了基因数目的统计，同时使用MA plot展示DEG的分布。然后根据差异基因检测结果，对其Gene Ontology (GO)功能以及KEGG (Kyoto Encyclopedia of Genes and Genomes)生物通路分类富集分析。
2.6. 差异表达基因的蛋白互作分析
为了研究leaf_1和leaf_2优先积累的蛋白质网络，使用了在线数据库STRING (https://string-db.org/)进行了差异表达基因的蛋白互作分析(Protein-protein interaction, PPI)。通过PPI分析，具有相互作用的DEG通常具有相似的功能。根据STRING蛋白互作数据库，对每组差异表达基因进行蛋白互作分析。
3. 结果与分析
3.1. 转录组测序与组装分析
对新疆橡胶草幼嫩和成熟叶片分别进行转录组测序，两个组织分别获得了214.16 Mb和211.67 Mb的Raw reads量，经过过滤后reads质量，最终获得幼嫩叶片165.94 Mb和成熟叶片166.46 Mb的clean reads (表1)。
Table 1. Quality statistics of filtered reads
表1. 过滤后的reads质量统计
Total Raw Reads (Mb)：过滤前的reads数；Total Clean Reads (Mb)：过滤后的reads数；Total Clean Bases (Gb)：过滤后的碱基总数；Clean Reads Q20 (%)：过滤后的reads中质量值大于20的碱基数占总碱基数的百分比；Clean Reads Q30 (%)：过滤后的reads中质量值大于30的碱基数占总碱基数的百分比；Clean Reads Ratio (%)：过滤后的reads的比例。
Table 1. Reference genome comparison result statistics
续表1. 参考基因组比对结果统计
aSample：样品名；bTotal Clean Reads：过滤后的reads总数；cTotal Mapping Ratio：比对上参考基因组的clean reads比例；dUniquely Mapping Ratio：唯一比对上参考基因组某一位置的clean reads比例。
过滤后的reads中质量值大于20和30的碱基数分别占总碱基数的百分比的96.36%和88.90%。过滤后的reads的比例占78%。转录本的质量指标测定，其GC值分别为40.81%和40.76%。
3.2. 橡胶草叶片转录组与参考基因组的比对分析
测序结果可用于后进行参考基因组对比，共获得332,396,368条clean reads，其中种植30 d幼苗叶片(leaf_1)共有165,935,798条，种植120 d后的成熟叶片(leaf_2)共有166,460,570条，与对比参考基因组分别占74.73%和67.66%。在Uniquely匹配时其中匹配率分别为37.76%和33.63%。参考基因组中总基因数为47,643条基因，转录测定发现有新的转录本分别为(leaf_1) 35,663条和(leaf_2) 34,119条，还预测了新基因数分别为14,890条和14,009条(表2)。
Table 2. Statistics of the number of genes and transcripts
表2. 基因、转录本数目统计表
使用Trinity对clean reads进行组装并获得转录本(Transcript)，leaf_1与leaf_2的Transcript总数分别为66,082个和65,187个，平均长度分别为737 nt和785 nt，其中200~1000 nt范围的转录本分别有83,461个和94,831个，占总量的75.65%和72.98%；1000~2000 nt的Transcript分别有19,547个和25,608个，占总量的17.72%和19.71%；大于2000 nt的Transcript分别有7316个和9506个，占总量的6.63%和7.32%。
对转录本进行聚类去冗余得到leaf_1和leaf_2的Unigene，leaf_1与leaf_2 Unigene的总数分别为66,749个和76,456个，平均长度分别为917 nt和955 nt，其中200~1000 nt的Unigene分别有44,108个和48,959个，占总量的66.08%和63.04%；1000~2000 nt的Unigene分别有16,193个和19,677个，占总量的24.26%和25.74%；大于2000 nt的Unigene分别有6448个和7820个，占总量的9.66%和10.23%。GC (%)含量分别为40.92%和40.74%。本实验对转录组数据的组装质量、Transcript和Unigene的长度分布情况进行了统计和分析(表1，表2)，结果发现为300~1000 nt长度的Unigene所占比重较大，说明测序质量较高可用于后续Unigene的功能注释。
3.3. 橡胶草叶片转录组All-Unigene的CDS预测
对新疆橡胶草叶片的所有Unigene中的CDS通过进行Blast比对Swiss-Prot数据库和Hmmscan，搜索Pfam蛋白同源序列，从而预测编码区域。从序列长度分布中我们共获得了55,120条CDS序列片段(图1)，总长度为52,588,077 nt，其中100~1000 nt的有36,043条占65.40%，1000~2000 nt的有14,799条占26.85%，2000~3000 nt的有3110条占5.64%，大于3000 nt的有1164条占2.11%。
Figure 1. CDS data length distribution for transcriptome of T. koksghz leaves
图1. 橡胶草叶片转录组CDS数据长度分布
3.4. 橡胶草叶片转录组Unigene功能注释
对测序后组装得到的橡胶草叶片转录组All-Unigene分别注释到七大功能数据库NR、NT、GO、Swiss-Prot、InterPro、KOG、KEGG的数据库(图2(A))，对每个数据库注释的Unigene数目进行统计，共有99,253条Unigene有对应的功能信息，其中，在NR中有58,402条，占总的58.84%；NT有29,578条，占总的29.80%；GO有19,809条，占总的19.96%；Swiss-Prot有38,614条，占总的38.90%；InterPro有47,437条，占总的47.79%；KOG有45,893条，占总的46.24%；KEGG有43,383条，占总的43.71%。
橡胶草叶片转录组Unigene的NR功能分类
已获得橡胶草叶片组装的All-Unigene注释到NR数据库，统计共有58,402条Unigene被注释结果。在菜蓟(Cynara cardunculus var. scolymus)、葡萄(Vitis vinifera)、芝麻(Sesamum indicum)、中果咖啡(Coffea canephora)和其它物种中都有同源序列分布。其中，Unigene与菜蓟相似序列占有61.35%；与葡萄相似序列占4.06%；与芝麻相似序列占有2.42%；与中粒咖啡相似的则有2.08%；还有30.09%的Unigene属于其他序列，可能还包含了橡胶草自身特有的基因序列(图2(B))。
Figure 2. All-Unigene species distribution in the database. A: All-Unigene annotated to the seven functional databases including NR, NT, GO, Swiss-Prot, InterPro, KOG, KEGG database and functional classification, respectively. B: Distribution of unique gene sequences of All-Unigene and cardoon, grape, sesame, medium grain coffee and T. koksghz
图2. All-Unigene在数据库中的物种分布。A：All-Unigene分别注释到七大功能数据库NR、NT、GO、Swiss-Prot、InterPro、KOG、KEGG的数据库和功能分类。B：All-Unigene与菜蓟、葡萄、芝麻、中粒咖啡和橡胶草自身特有的基因序列分布情况
3.5. 橡胶草叶片差异基因表达量分析
将橡胶草叶片leaf_1与leaf_2进行差异基因表达基因数目比对统计，发现leaf_1与leaf_2对比中12,637个表达基因上调，8103个基因表达下调，33,509个基因没有差异表达(图3)。
Figure 3. Statistical map of differentially expressed genes
图3. 差异表达基因数量统计图
3.6. 橡胶草叶片差异基因的功能注释
为进一步了解橡胶草幼叶片及成熟叶片响应的相关分子机制，对Differentially expressed gene (DEGs)进行Gene Ontology (GO)富集分析。GO主要分为生物过程(biological process)、分子功能(molecular function)和细胞组分(cellular component)三大类别。leaf_1与leaf_2的差异基因表达经GO分析分别分布的三大类别中，共有52个类型，其中生物过程中包括代谢过程(2750)、细胞过程(2381)和单一生物过程(1525)；细胞组成包括细胞部分(2067)、膜组成(1445)和组织部分(1264)；分子功能包括催化活性(3131)、结合(2538)和运输活性(2538)等。另外，大多数差异基因与生物过程相关，而参与分子功能的差异基因相对较少。在细胞组分中，以参与细胞、细胞组分、膜、膜组分、细胞器构成的富集最明显。参与生物过程的差异基因主要富集在代谢过程和细胞过程这两个类别；而分子功能中催化活性和结合是leaf_1和leaf_2分子功能差异基因富集的主要类别(图4)。
3.7. 差异表达基因的Pathway功能分析
为了进一步研究差异基因的生物学功能，进行了KEGG通路分析。分析差异基因KEGG代谢通路结果显示，差异表达基因能映射到新陈代谢(Metabolism)、遗传信息处理(Genetic Information Processing)、环境信息处理(Environmental Information Processing)、细胞进程(Cellular Processes)和有机系统(Organismal Systems)的5大类19个一级通路主要途径，其中涉及代谢通路的基因最多，共有10,434个，占到注释KEGG中的基因比例的61.78%。其次是遗传信息处理，占KEGG中的基因比例的23.23% (图5)。一级通路中包括翻译过程有777个基因，碳水化合物代谢有1572个基因，折叠基因1375个，环境适应662个基因，整体代谢和概观代谢基因3761个，氨基酸代谢948个基因，转运和分解代谢基因786个，脂质代谢基因724个等(图5)。这些结果表明，橡胶草幼嫩与成熟叶片差异表达基因主要与折叠与降解，翻译和碳水化合物代谢有关。
Figure 4. GO functional classification of the leaf_1-VS-leaf_2 differentially expressed genes
图4. Leaf_1-VS-leaf_2差异基因GO功能分类图
除此之外，还分析了前20个KEGG亚通路(图5(B))，结果显示，其中包含3525种基因的最丰富的途径是碳代谢。第二个包含2032个具有次级代谢产物合成的基因。第三子途径具有用于氨基酸的生物合成的675种基因，其次是用于嘧啶和嘌呤代谢的642和641种基因，以及内质网蛋白合成中的579种基因。这些KEGG亚途径中还富含涉及苯丙素和蔗糖代谢，糖酵解，苯丙素类生物合成以及氨基糖和核苷酸糖代谢。这些基因对于许多植物，特别是次级代谢产物和天然橡胶生产植物中类异戊二烯的生物合成至关重要。
Figure 5. KEGG Pathway classification chart of differentially expressed genes. A, the X-axis represents the proportion of the gene, the Y-axis represents the KEGG functional classification. B, top 20 of sub-pathway enriched. The X-axis represents the enrichment factor value (Rich Factor), the Y-axis represents the path name. The color of the dots in the figure represents the Q-value, and the lighter the color, the greater the value, the deeper the smaller, the latter represents the more significant the enrichment results; the size of the dots represents the number of DEGs, and the larger the number represents, the smaller the number represents
图5. 差异基因KEGG Pathway分类。A，leaf_1-VS-leaf_2轴表示基因所占比例，Y轴表示KEGG功能分类。B，前20个亚KEGG通路富集。X轴代表富集因子值(Rich Factor)，Y轴代表通路名称，图中圆点的颜色代表Q-value，颜色越浅值越大，越深值越小，值越小代表富集。结果越显著；圆点的大小代表DEGs数目，点越大代表数目越多，越小代表数目越少
3.8. 差异表达基因的蛋白互作分析
橡胶草叶片的蛋白质–蛋白质相互作用(PPI)网络分析可以反映蛋白质数据集的关键要素。为了研究橡胶草幼嫩叶片leaf_1和成熟叶片leaf_2优先积累的蛋白质的网络，使用了在线数据库STRING。对于leaf_1和leaf_2优先积累的蛋白质，可以观察到18个网络簇，其中8个网络簇由不同的颜色代表上调与下调，不同大小表示关系紧密程(图6)。其中最大的簇中的这些基因的蛋白质主要由油菜素类固醇相关蛋白组成。尤其是小亚单位核糖体蛋白和大亚单位核糖体蛋白，以及天冬氨酸氨基甲酰转移酶具有相互作用网络的大多数关键节点。
4. 讨论
早在1931年，人们就已经意识到俄罗斯蒲公英具有产胶、食品及其它副产物等应用价值 [14]。橡胶草叶内含次生代谢产物和花青素，能够改善关节柔韧性，增强血管强度和抑制炎症等药用价值 [9] [28]。随着基因组、转录组、蛋白组和代谢组测序技术不断更新合完善，组学研究为挖掘相关调控机制关键基因/蛋白提供了重要技术方法。其中近期快速发展的转录组测序(RNA-Seq)是用于研究功能基因组学分析技术。目前橡胶草尤其是原产地新疆橡胶草叶片的转录组信息较少，利用茉莉酸刺激橡胶草产胶转录组分析发现，IPP是三萜和橡胶合成的关键中间物质，三萜类物质和橡胶是通过IPP的转换互相合成 [16] [23]。比较LR和HR橡胶草根转录组，也说明次生代谢产物与橡胶合成呈负相关 [25]。
Figure 6. Network diagram of proteins interaction. Red represents up-regulated proteins, blue represents down-regulated proteins, and the size of the circle represents the number of interaction relationships. The larger the circle, the denser the relationship
图6. 蛋白互作网络图。红色表示上调，蓝色表示下调，圆圈大小表示相互作用的关系的个数，圆的大小表示表关系的强度，圆圈越大表示关系越密集
AtPAP1在短角蒲公英中的异源表达导致营养组织显示出红色/紫色色素沉着表型，并积累了大量的花色苷。该表型与在拟南芥中过表达AtPAP1所描述的表型一致，并且与在烟草，甘蓝型油菜和番茄中的异源表达相类似 [29] [30] [31]。次生代谢产物是植物中药用活性成分的重要组成部分 [32]，通过构建橡胶草不同组织器官的转录组文库，能更全面地研究TKS的次生代谢产物生物合成途径中关键基因的具体功能机制，定向提升次生代谢产物要用价值 [33]。但次生代谢产物与橡胶合成的具体调控机制未明确阐述。由于橡胶草转基因相关研究数据缺乏，橡胶合成关键基因与橡胶合成精确调控还需多方面验证。本研究通过RNA-Seq测序共获得了33,205 Mb的clean reads，将其与参考基因组对比后，匹配率70%以上的基因在基因组中匹配到，再通过幼嫩叶片与成熟叶片差异基因表达分析，并注释到KEGG的代谢通路中，其中花青素与萜类生物合成相关较多。蛋白互作预测分析油菜素类固醇合成核糖体蛋白合成途径成为主要互作网络关键点。成功地建立了TKS叶片转录组高质量的数据，以精确地区分幼嫩叶片和成熟叶片在转录水平的理论分析，获取更多的转录组数据将有助于从数据集中挖掘的一些实例，将为通过生物技术手段，提高橡胶草叶片次生代谢产物的产量和有效成分含量以及活性成分的合成途径以及更深入地了解橡胶草叶片的生理作用和生物学理论研究提供理论基础。
基金项目
国家自然科学基金地区基金项目(32060072)；石河子大学高层次人才科研启动项目(XJ2020001202)。
NOTES
*通讯作者。