报春花科3种植物对青藏高原适应性进化的转录组学研究*

郝豆豆，张勇群，施静，拉多，雷鸣

(1.西藏自治区人民政府驻成都办事处医院，四川 成都 610041；2.西藏大学，西藏 拉萨 850001；3.青藏高原生物研究所，西藏 拉萨 850012)

青藏高原拥有高寒干旱、强辐射、低分压的环境现状，生活在青藏高原的生物面临着来自其恶劣环境的各种压力[1]。气候、地貌、生物间的相互作用等因素决定了适宜物种分布的地理范围[2]，青藏高原作为独特的地质-地理-生态单元，拥有独特的生物资源，成为了世界生物多样性研究的热点地区。报春花科(Primulaceae)植物在我国高原广泛分布，它与龙胆科(Gentianaceae)、杜鹃花科(Ericaceae)一起被称为高山三大植物类群。该科植物种类繁多，对于维持高寒生态系统稳定方面发挥着重要作用[3]，且很多种植物在藏药中具有较高药用价值[4]。报春花属(Primula)植物花色艳丽，富含黄酮类化合物，这些化合物有益于人体健康，是医药和化妆品行业的重要资源[5]。点地梅属(Androsace)是极端环境生境中的高山垫层植物，能够促进其他植物物种生长，改善土壤微生物群落，并增加植物、节肢动物和微生物的丰度和多样性[6]。白粉圆叶报春(Primulalittledalei)、柔小粉报春(P.pumilio)和禾叶点地梅(Androsacegraminifolia)是高海拔陆生多年生草本植物，其植株矮小、花色艳丽，大多分布在青藏高原海拔3 800～5 300 m的岩石缝隙中、草甸、碎石坡地，具有较高的观赏价值，在藏医药中全草及花可以入药[4]。

由于青藏高原恶劣的环境和复杂的古历史气候，近年来在生物如何适应高海拔环境和多样性形成机制研究方面已成为热点地区[7-8]。在相似的环境中，生物体经过选择压力进化出相似的适应性机制。转录组测序(RNA-seq)技术可以获取生物体在某个时间的所有RNA信息，并提供了一种快速有效的方法来获取大量的蛋白质编码基因，通过鉴定物种之间同源基因，可用于探索非模式生物适应生存环境的分子遗传机制[9-10]。

本研究使用Illumina的双末端RNA测序(RNA-seq)对白粉圆叶报春、柔小粉报春和禾叶点地梅进行全基因组表达研究，并且通过对3个物种的正选择单拷贝基因的分析，初步探讨调节这3个物种适应高海拔环境的生物学途径和分子机制，确定与高山环境适应相关的潜在候选基因。

1 材料与方法

1.1 试验材料

于2020年7月在中国西部的西藏自治区拉萨市夺底沟采集了白粉圆叶报春、柔小粉报春和禾叶点地梅的新鲜叶片，将150 mg的新鲜叶片剪碎，放入装有RNALater溶液的2 mL冻存管，并于-80 ℃的低温冰箱中保存，直到提取RNA。

1.2 转录组测序

用mirVana miRNA分离试剂盒(Ambion)提取高质量的总RNA，使用TruSeq链式mRNA LTSample Prep试剂盒生成测序文库，使用Illumina HiSeq X Ten对文库进行测序，文库构建和测序工作是在oebiotech(中国上海)公司完成，原始测序数据已提交到NCBI数据库。

1.3 序列组装和功能注释

3种植物的叶片转录组mRNA测序原始数据使用Trimmomatic程序处理[11]。去除长度小于50 bp且含有ploy-N的片段，获得cleanreads；使用Trinity程序[12]进行从头组装，去除<200 bp的转录本。使用Diamond程序与NCBI非冗余(NR)、SwissProt、 eggNOG、Pfam、GO分类和真核完整基因组数据库进行比对来注释基因的功能，在KEGG数据库中注释其潜在的代谢途径[13-14]。

1.4 三个物种直系同源基因的鉴定和选择压力分析

使用OrthoMCL软件[15]鉴定3个物种的编码基因的同源基因簇(正交群orthogroups)，排除编码长度小于50个氨基酸的基因，然后筛选出一对一的单拷贝直系同源基因用于后续分析。使用KaKs_Calculator(http://code.google.com/p/kaks-calculator/wiki/KaKs_Calculatoref)软件包[16]计算非同义突变率(Ka)，同义突变率(Ks)和Ka/Ks值。Ka/Ks值可用于确定是否有选择压力作用于蛋白质编码基因，如果Ka/Ks>1，则考虑为正选择效应，该基因称为正选择基因；如果Ka/Ks=1，则认为存在中性选择；如果Ka/Ks<1，则认为该基因受到纯化选择/负选择。 Ka/Ks值越大，序列中发生的非同义突变越多，这样的基因是快速分化的基因，这对于物种的进化非常重要；与此相反，当Ka/Ks为0或接近于0时，在序列中主要发生同义突变，这些基因相对保守，并且它们可能在进化中受到强大的选择约束[17]。Fisher’s exact test用于验证每个直系同源物的Ka和Ks值，正选择基因的Ka/Ks值必须大于1，P值<0.05[18]。最后对筛选出的正选择基因进行GO和KEGG功能富集分析，注释基因功能和代谢途径。

2 结果与分析

2.1 转录组测序及质控结果

对白粉圆叶报春、柔小粉报春和禾叶点地梅3个物种的新鲜叶片进行转录组测序，对原始数据进行质量控制后，总共获得20.94 Gb的数据和145 269 464 clean reads，所有reads的GC含量平均值为45.16%，Q30为93.85%,见表1。

表1 报春花科3种植物测序数据概览

转录组测序结果符合后续组装分析的要求。原始数据已提交到NCBI(SRA)数据库中，白粉圆叶报春、柔小粉报春和禾叶点地梅3个物种对应生成的数据查询ID分别为SRR9110566、SRR9110567和SRR9112929。

2.2 转录组组装和unigene功能注释

从头组装质控后的核酸片段，3个物种中一共得到了51 666个unigene。Unigene的平均长度为920bp，ContigN50的长度为1 158bp(N50 = 1 158bp)(表2)，说明这3个物种的转录组测序得到的序列组装结果好，测序质量高，满足后续的分析要求。将3个物种的所有unigene比对到蛋白质数据库(包括NCBI非冗余(NR)、SwissProt、eggNOG、Pfam、GO分类和KEGG代谢通路数据库)。据每个数据库中成功比对注释的基因数目和占比(表3)，NR数据库注释到的Unigenes的数量最多。3个物种分别从蛋白质数据库中预测得到的编码序列(CDS)数目分别为15 213、16 781和16 940 个(表4)。在数据库中未注释到的unigene数目分别为649、793和1 290个，所占百分比分别为4%、4.5%和7.1%整体的功能注释成功率较高。

表2 3种植物转录组组装信息汇总

表3 基于7个蛋白质数据库的基因功能注释

表4 编码序列数量的统计

2.3 正交群鉴定和正选择分析

用OrthoMCL软件分析鉴定3个物种编码基因的同源基因簇(正交群)，3个物种共有的正交群(orthogroups)为7 030个，其中白粉圆叶报春(989个)和柔小粉报春(679个)特有的正交群数量相近，并且显著少于禾叶点地梅(1 919个)特有的正交群数量。这表明白粉圆叶报春和柔小粉报春的亲缘关系更近(图1) 。从3个物种共有的7 030个正交群中鉴定出2 151个包含单拷贝直系同源基因的正交群，将这些直系同源基因用于非同义(Ka)和同义(Ks)替换进化分析，其中有827个单拷贝直系同源基因的 Ka/Ks(ω)大于1、Fisher testP<0.05，表明这些基因受到了正选择作用；而且，其中有553个单拷贝直系同源基因的Ka/Ks(ω)大于1、Fisher testP<0.01。

图1 3个物种共有和特有的同源基因簇

2.4 3个物种的正选择单拷贝直系同源基因功能注释

对827个受到正选择的单拷贝直系同源基因进行了GO富集分析，包括细胞成分(CC)，分子功能(MF)和生物过程(BP)(图2)。

图2 正选择基因的GO功能注释

三大功能的生物学过程中包含的基因数目最多，分子功能包含的最少。在不同的功能群中发现了一些基因表现出快速进化模式，具有强选择(Ka/Ks>1)的迹象，很可能与环境适应和应急反应有关，因为其中大多数基因参与编码植物代谢产物生物合成途径中的载体(GO0009921，GO0005458)或酶(GO1903600谷氨酰胺酶复合物，GO0000276 ATP酶合成复合物，GO0004076生物素合成酶，GO0033984、GO0050589花青素合成酶，GO0003998酰基磷酸酶，GO0003842脱氢酶)，还有很多受到正选的基因与遗传信息的加工有关(GO:0001105，GO:0001158)。

KEGG富集分析结果包括四类：细胞过程(CP)，环境信息处理(EIP)，遗传信息处理(GIP)，代谢(MB)(图3)。共有684个正选择基因在KEGG数据库中成功注释，涉及23个代谢途径。受到正选择的单拷贝同源基因的功能主要富集在运输、分解代谢、信号转导(信号肽酶复合物亚基、信号肽酶复合物催化亚基SEC11、Sec依赖性信号转导途径)、遗传信息加工(复制和修复、折叠、转录、翻译)、能量代谢(核苷酸代谢、碳水化合物代谢、辅因子和维生素代谢、脂质代谢)(图4)。参与了嘌呤和嘧啶代谢的3个受到正选择的基因——肌苷5’-单磷酸脱氢酶(guaB)、次黄嘌呤磷酸核糖基转移酶(HPRT1)和胞苷脱氨酶(cdd)——也参与了药物代谢。同时，鉴定出8个受到正选的基因与过氧化物酶体代谢通路相关，包括过氧化物酶体膜蛋白输入受体(PEX3、PEX19、PEX13、PEX14)，过氧化物酶体膜蛋白(MPV17)、对脂质运输和代谢的过氧化物酶体植酸酰辅酶A羟化酶(PHYH)和脂肪酰基辅酶A还原酶(FAR)、参与视黄醇代谢的脱氢酶/还原酶SDR家族成员4(DHRS4)，以及3种与类黄酮生物合成相关的酶——查尔酮异构酶、黄酮醇合酶和花青素合酶。将富集到KEGG结果中的PSG基因的Ka/Ks值进行排序，表5中展示了受到正选择压力最强的前30个基因。

图3 正选择基因的KEGG通路分类

图4 正选择基因在KEGG富集分析结果中的前20项注释结果

表5 KEGG富集结果中根据Ka/Ks值排序后的前30个基因

续表5

3 讨论与结论

在自然界中，植物会面临很多的胁迫方式，例如盐碱、重金属、干旱、高温、寒冷、强辐射等等，对其生存造成威胁[19]。在有害胁迫的刺激下，植物体内会产生响应胁迫机制，如：生物碱类、精氨酸和脯氨酸代谢在猴樟(Cinnamomumbodinieri)响应盐碱胁迫机制中起主要作用[20]；毛白杨(Populustomentosa)悬浮细胞受到氧化胁迫，抑制了磷酸酶活性导致线粒体功能异常，于激素相关的基因表达下调，进而影响生长发育[21]。通过直系同源基因的选择压力分析，可以获得植物在胁迫环境下的受到不同选择作用的基因表达模式[22]。

白粉圆叶报春、柔小粉报春和禾叶点地梅3个物种中很多受到正选择作用的单拷贝直系同源基因与恶劣的高原环境有关，并参与了多种生物过程。其中28个受到正选择的基因(positively selected genes,PSG)参与能量代谢(如碳水化合物、运输)，这对于3种植物应对青藏高原上的寒冷的气候至关重要。受到环境胁迫时，细胞内部产生氧化胁迫，因此，激素调节、抗氧化等相关通路被激活，从而使植物适应生长环境[23]。过氧化物酶体的功能复杂，很多重要的代谢活动都发生在过氧化物酶体中，例如活性氧(ROS)和活性氮(RNS)的基础代谢，以及脂质代谢，激素生物合成、代谢的生理过程[24]。本研究鉴定出了8个受到正选择的基因(SOD1、DECR2、MPV17、SOD2、MLYCD、PEX19、PXMP2、MPV17)参与过氧化物酶体反应。过氧化物酶体具有很强的增殖能力，并能根据周围的细胞状态进行生化适应，从而帮助植物适应环境变化。类黄酮参与植物对刺激的响应，鉴定出3个受到正选择的基因编码与类黄酮生物合成相关的酶，黄酮类是植物的主要次生代谢物，是花青素生物合成和代谢的重要前体[25]，还具有抗氧化剂活性[26,28]，并提供针对生物(如昆虫、食草动物、病原体)和非生物胁迫(如紫外线辐射、多余光能的消散)的保护。

此外，受到正选择的单拷贝同源基因(如RBMX、HNRNPG、SNRPF、SMFRP-L27、MRPL27、rpmA、RP-L11、MRPL11、rplK、RP-S29e、RPS29、RP-L36、MRPL36、rpmJ等)中有很多与遗传信息的加工相关，如编码RNA聚合酶蛋白成分(RPB4、RPB9、RPABC3、RPC11，RPA49)，参与RNA加工和修饰、mRNA剪接/加工因子CDC5/CEF1、38B和SPF27，在青藏高原各种恶劣环境因素造成的基因损伤修复过程中发挥重要的作用；3个受到正选择的基因(LHCA1，LHCA2和LHCB3)与光收集复合体(LHC)同源，构成了光系统I(PSI)和PSII的天线系统[27]，这些正选择基因可能在响应强烈的太阳辐射和青藏高原的长时间日照方面发挥重要作用。

生长素是通过控制细胞分裂和分化来调节植物生长和发育，生长素在组织内的分布差异是通过极性生长素转运(PAT)机制介导的，PIN1样生长素转运蛋白的基因在受到了正选择作用(Ka/Ks> 1、P=2.05E-07)，该基因在控制植株高度中扮演重要的角色，白粉圆叶报春、柔小粉报春和禾叶点地梅植株矮小，有利于抵抗青藏高原的各种环境胁迫(例如强紫外线、低温、强风、食草动物等)。

综上所述，这3种植物为了适应复杂严酷的高海拔生境而进化出共同的适应机制——受到正选择的单拷贝同源基因参与了碳水化合物运输和代谢、氨基酸代谢、氧化磷酸化、脂肪酸生物合成和代谢、嘌呤和嘧啶代谢、次生代谢(类黄酮生物合成)和遗传信息加工及翻译后修饰。这些生物过程对于这3个物种应对青藏高原恶劣环境(紫外线辐射、气温低、强风、贫瘠的土壤)至关重要。本研究从分子生物学层面揭示了植物对高海拔环境的适应机制，为高山植物遗传进化的分子机制研究提供了思路。