基于SWChen 系统预测靶向林麝(Moschus berezovskii) FASN、HMGCR 和CYP11A1基因miRNAs的研究*

谌颖莲，曾德军，赵贵军，封孝兰，张承露，吴佳勇，竭 航

(重庆市药物种植研究所，特色生物资源研究与利用川渝共建重点实验室，重庆 408435)

麝香是一种珍贵的动物源药材，具有较高的药用价值和香料价值[1-4]，但其形成过程复杂且未知，限制了麝香产业的发展和产香效率的提高。研究表明：麝香的化学组分包括脂肪酸类的烷烃物质或衍生的醇类、酮类和醛类 (如麝香酮、6-甲基庚烷-1,6-二醇和十六烯醛等)以及固醇类、甾类和芳香族类物质 (如胆固醇、间甲酚和二氢雌甾酮等)[2]，这些化学组分对麝香形成至关重要，同时也是麝香药用和香用价值的关键物质。动物的脂肪酸合酶(fatty acid synthase，FASN)、3-羟基-3-甲基-戊二酰辅酶A 还原酶(3-hydroxy-3-methylglutaryl-CoA reductase，HMGCR)和细胞色素P450 家族11 亚家族A1 (CYP11A1)是长链的脂肪酸类衍生物、大环碳链酮类以及芳香族甾类物质合成的关键酶[5-9]。基于麝香的化学组分和合成相关化学组成的关键酶，在林麝(Moschus berezovskii)中筛选调控FASN、HMGCR和CYP-11A1基因表达的调控因子，对探究麝香可能的生物合成机制具有重要意义。

MicroRNA (miRNA)是一种长度约为20 nt的非编码RNA，通过作用靶基因mRNA的3′-UTR 在转录后水平抑制靶基因的表达，对动物的生长和代谢等过程具有关键的调控作用[10-11]。筛选调控林麝FASN、HMGCR和CYP11A1基因表达的miRNAs 能够进一步了解麝香的生物合成过程，然而现在已知的林麝基因序列信息非常少，缺少林麝基因的mRNA 和miRNA 序列信息[1,12]。由于牛(Bos taurus)、山羊(Capra hircus)和绵羊(Ovis aries)与林麝编码基因的3′-UTR 序列并不完全一致，通过同源物种来推测调控林麝基因的miRNAs 获取信息并不完整，缺少适合预测林麝miRNAs的工具。本研究通过基于Smith-Waterman 局部序列比对算法[13]编写的SWChen系统对调控林麝FASN、HMGCR和CYP11A1基因表达的miRNAs 进行预测筛选，并以人(Homo sapiens)的FASN基因为例，与TargetScan、mi-RDB 和miRSytem 预测软件和miRTarBase 数据库比较，评价SWChen 系统对miRNAs 预测能力；利用代谢组学测序来确定麝香的化学成分，对与合成相关化学成分有关的FASN、HMGCR和CYP11A1基因进行miRNAs 预测，以林麝、牛、山羊和绵羊预测的结果进行集合论分析确定候选miRNAs，以期为麝香的生物合成机制研究提供参考。

1 材料与方法

1.1 试验动物

118 头健康成年雄性林麝 (5 岁) 饲养于重庆市药物种植研究所林麝养殖基地(N29°，E107°，海拔678 m)，饲养采用精饲料和青绿饲料[红薯藤 (Ipomoea batatas)和光叶海桐 (Pittosporum glabratum)嫩叶]相结合的方式，自由采食。精饲料包括 65%玉米、25%大豆、6%麦麸、0.4%食盐、1.5% CaCO3、0.8% CaHPO4、0.15%维生素和0.15%矿物质。每天17：30 混合精饲料与红薯藤等补充饲喂。

1.2 样品采集

在5—6 月泌香盛期，采用微创香腺采集法获取林麝香腺组织样本，立刻放入RNA 保存液中并转移至-80 ℃保存；在11—12 月麝香成熟阶段，通过“活麝取香”法采取林麝成熟麝香3 份，置于离心管内保存。

1.3 代谢组检测

将采集到的麝香样品送至上海百趣生物有限公司采用GC-MS 进行非靶标代谢组学检测，使用甲醇—氯仿(体积比为3∶1)溶液进行萃取，通过总离子流图信号获取麝香的主要化学成分及相对含量。

1.4 转录组测定及序列分析

通过总RNA 提取后，送样至上海生工生物科技有限公司采用Hi-seq 2500 测序。然后将测序数据进行Trinity 拼接、序列比对和功能注释等获得FASN、HMGCR和CYP11A1基因mRNA 序列并进行结构分析，确定3′-UTR 序列以及林麝miRNAs 序列数据；以牛、山羊和绵羊为参考，对3 个基因的3′-UTR 序列进行比对，并使用MEGA 6.0 及最大似然法构建进化树；通过RNA 二级结构分析(unafold.rna.albany.edu)了解各物种基因3′-UTR 序列的差异情况。

1.5 SWChen 系统对miRNAs 预测准确性的验证

由于目前常用的miRNAs 预测软件如miRSystem (mirsystem.cgm.ntu.edu)、TargetScan (www.targetscan.org)和miRDB (mirdb.org)仅预测部分常见物种的miRNAs，无法对调控林麝基因的miRNA 进行准确预测，因此，本研究在miRBase 数据库(www.mirbase.org)下载人的miRNAs序列，并在GenBank (www.ncbi.nlm.nih.gov/genbank/)查询人FASN基因3′-UTR 序列(NM_0041 04.5)，利用SWChen系统 (https://github.com/YLCHEN1992/SWChen)对调控FASN表达的miRNAs 进行预测，通过与在线软件miRSystem、TargetScan 和miRDB的预测结果比较，以miRTarbase 数据库(mirtarbase.cuhk.edu.cn)中调控FASN表达的全部miRNAs (miRTarget_All)和具有强力证据可调控人FASN基因表达的miRNAs(miRTarget_High)作为参考，检测SWChen 系统对miRNA 预测的准确性。

1.6 调控FASN、HMGCR 和CYP11A1 基因表达的miRNA 预测

将得到的林麝mRNA 3′-UTR 序列及获得的林麝miRNAs 数据库利用SWChen 系统进行分析，预测调控林麝FASN、HMGCR和CYP11A1基因表达的miRNAs；在GenBank 中获取牛、山羊和绵羊FASN、HMGCR和CYP11A1基因的mRNA序列(表1)。在miRBase 数据库获取牛、山羊和绵羊的miRNA 序列，由于绵羊miRNA 数量较少，将山羊和绵羊miRNA 进行整合得到羊的miRNA数据库，共同用于山羊和绵羊miRNA的预测分析。利用SWChen 系统预测调控牛、山羊和绵羊FASN、HMGCR和CYP11A1基因的候选miRNAs，并通过双阈值(种子序列匹配评分及亲和程度评分同时考虑)预测调控基因表达的miRNAs，以增加预测精度，最后将预测结果进行集合论分析以期获取麝香合成相关和候选miRNAs 分子。

表1 物种基因序列GenBank 编号Tab.1 GenBank No.of species gene sequence

2 结果与分析

2.1 SWChen 系统miRNA 功能预测结果

由图1 可知：SWChen 通过种子区匹配预测(single)调控人FASN的基因共有193 条，其中匹配miRTarget_All 数据8 条，匹配miRTarget_High数据3 条，而调高miRNA 与mRNA 亲和力判定阈值后(双阈值)与miRTarbase 数据库不匹配；TargetScan 在线软件预测有745 条，其中匹配miRTarget_All 数据15 条，匹配miRTarget_High数据3 条；miRSystem 在线软件预测有62 条，其中匹配miRTarget_All 数据4 条，匹配miRTarget_High 数据4 条；miRDB 在线软件预测有47条，其中匹配miRTarget_All 数据6 条，匹配miRTarget_High 数据2 条。4 个软件预测的miRNAs结果集合统计(图2)显示：TargetScan 出现530条独立miRNAs，miRSystem 和miRDB 分别出现10 条和1 条独立的miRNAs；SWChen 系统总共出现3 条独立的miRNAs，有186 条与Target-Scan 结果重叠，31 条与miRSystem 结果重叠，39 条与miRDB 结果重叠，双阈值SWChen 结果除miRSystem 均包含在其他软件预测结果中。

图1 4种miRNA 预测软件与miRTarget 数据库匹配情况Fig.1 Comparing prediction results of four miRNA prediction softwares by using miRTarget database match

图2 4种miRNA 预测软件预测结果集合统计Fig.2 The gather statistic to resuls of four miRNA prediction softwares

2.2 林麝成熟麝香代谢组测序分析结果

林麝麝香化学组分含量最高的前10 个物质列于表2。结果显示：成熟麝香中的主要成分为长碳链的酮和醛等物质、特殊脂肪酸碳链的环酮类物质以及芳香族和甾类物质。

表2 成熟麝香主要化学组分Tab.2 The main chemical components of matured musk

2.3 物种基因序列进化树构建及RNA 二级结构分析

由图3 可知：同一基因之间3′-UTR 序列存在差异较小，其中林麝FASN和HMGCR基因3′-UTR 序列与牛和羊被分为两类，而林麝CYP11A1基因3′-UTR 序列与山羊和绵羊被归为一类。进一步对各基因的3′-UTR 序列进行RNA 二级结构进行预测，结果(图4)显示：林麝、牛、山羊和绵羊FASN基因3′-UTR 序列RNA 二级结构的最小自由能分别为-222.80、-234.60、-229.10 和-234.70 kcal/mol，HMGCR基因3 ′-UTR序列RNA 二级结构的最小自由能分别为-466.00、-477.60、-267.20 和-1 417.30 kcal/mol，CYP11A1基因3′-UTR 序列RNA 二级结构的最小自由能分别为-110.10、-112.90、-105.50 和-112.00 kcal/mol，其中林麝和牛FASN基因3′-UTR 与山羊和绵羊相比出现较多的单链空泡。

图3 基因3′-UTR 序列进化树比对Fig.3 Sequence alignment of gene 3′-UTR by evolution tree construction

图4 各物种基因3′-UTR RNA 二级结构预测Fig.4 Predicted RNA secondary structures of gene 3′-UTR sequences in respective species

2.4 SWChen 系统分析各物种调控FASN、HMGCR和CYP11A1 基因表达的miRNAs

通过miRNA 转录组序列保守分析发现：林麝miRNAs 数据总有1 391 条，miRBase 数据库中查询牛的miRNAs 数据为1 030 条，山羊和绵羊的数据为519 条。SWChen 系统的种子匹配预测结果显示：调控林麝、牛、山羊和绵羊FASN基因表达的miRNAs 在相应数据库中的占比分别为5.25%、5.53%、1.70%和2.38%；调控HMGCR基因表达的miRNAs 占比分别为12.55%、13.11%、10.88%和6.12%；调控CYP11A1基因表达的mi-RNAs 占比分别为3.54%、3.69%、4.76%和4.22%；在林麝和牛中，调控FASN和HMGCR基因表达的miRNAs 占比高于山羊和绵羊，而调控CYP11A1基因表达的miRNAs 占比低于山羊和绵羊。集合论分析结果(图5 和表3)显示：在4 个物种中，调控FASN基因表达的miRNAs 交集为miR-24-3p 和miR-30b-3p；调控HMGCR基因表达的mi-RNAs 交集为miR-29b、miR-146a 和miR-2284d；调控CYP11A1基因表达的miRNAs 交集为let-7家族。

表3 调控基因的miRNA 交集Tab.3 miRNA intersection of targeting genes

图5 4 个物种中调控FASN、HMGCR 和CYP11A1 基因表达的miRNAs 集合分析维恩图Fig.5 Venn diagram of miRNAs regulating FASN, HMGCR and CYP11A1 genes expression intersection analysis in four species

双阈值分析结果显示(表4)：调控FASN基因表达的miRNAs 仅在山羊和绵羊中存在miR-874-5p 交集；调控HMGCR基因表达的miRNAs在林麝、牛和山羊中存在miR-545-3p、miR-429和miR-181b-3p 交集，调控CYP11A1基因表达的miRNAs 仅在林麝、山羊和绵羊中存在miR-532-3p的交集。

表4 在各物种中预测调控基因的miRNAs (双阈值)Tab.4 Predicted miRNAs targeting genes in species respectively (dual thresholds)

3 讨论

SWChen 系统是基于Smith-Waterman 局部序列比对算法和碱基互补配对原理[11,13]，预测2 条序列之间的相互作用能力和相似性情况的程序，主要基于R 语言软件编写，用于RNA 之间的作用情况预测[14-17]。计分规则参数使得序列比对过程中出现连续3～4 个的相互匹配才能抵消1 个错配，符合miRNA 序列种子匹配规则，通过矩阵最大值进行重新回数算分即为SWChen 系统的最终评分，选择适合的评分参数对于提高相互作用预测准确性具有重要意义。

Smith-Waterman 算法是评价2 条序列之间功能域和保守域的比对算法，相比BLAST 算法[18]和马尔可夫模型[19]推算，Smith-Waterman 算法更加准确，能够满足预测需求，但是这种算法的速度最慢，且会忽略第2 评分结合位置。miRTarBase数据库包含了现有已验证调控功能的miRNAs，而在所有物种的研究中，人FASN研究证据较多，因此以人种为对象评估SWChen 系统的准确性具有更高的可信度。程序评估结果显示：4种预测工具的结果差异并不大，其中SWChen 相比miRSystem 和miRDB 工具能预测出较多的候选miRNAs，匹配miRTarBase 数据库的miRNAs条数更多；而TargetScan 预测的miRNAs 条数约为SWChen的4 倍，但匹配数据库的miRNAs 约为SWChen的2 倍，提示SWChen的背景噪音相对TargetScan 较小。所有的miRNA 预测工具均基于miRNA 对靶基因的识别机制[11]，其中TargetScan 在计分过程中考虑结合位置的保守性和miRNA 亲和程度，再去除物种序列保守性判断之后预测的miRNAs 会增加，而SWChen 并不考虑序列保守性，同时将评分控制为80 分，即种子区域的完全匹配，占据识别机制中的大部分识别情况；而TargetScan的算法包括了miRNA 第2～8 位和第3～8 位的匹配规律[11]，因此预测出较多的候选miRNAs，但却增加了错误miRNAs的噪音，这是TargetScan 和SWChen 之间预测结果差异的原因。miRSystem 系统为集成多种预测工具的miRNA 预测系统，在预测结果中获得强有力证据的miRNAs 匹配条数最多，背景噪音相对较小，但匹配miRTarget_All 数据的数量也减少，在miRNA的靶基因研究中，miRNA 预测主要基于多工具的交集筛选，因此集成的miRSystem 会得到更多的miRTarget_High 匹配结果；双阈值预测能够得到较多的重叠miRNAs，暗示双阈值的准确性可能更高。SWChen 系统能够满足miRNA 靶基因预测的需求，能够得到相对其他工具同等价值的结果，且离线预测、数据库及参数的可变动使其更加灵活，能够预测林麝miRNAs 对靶基因的调控情况，可为特种动物、昆虫、植物和真菌等关键miRNA的筛选提供便利。

已有研究表明：麝香的形成过程主要有2种假设[1-4,12]。第一，林麝香腺分泌初香后，初香物质储存于香囊中经微生物的发酵熟化后形成成熟麝香；第二，林麝香腺直接分泌麝香的主要化学成分储存于香囊中积累，而微生物不是麝香形成的决定因素。无论林麝通过哪种方式形成成熟麝香，林麝合成和分泌的化学组分对麝香的形成至关重要。通过代谢组学的手段对麝香进行检测，使用不同的萃取试剂得到的化学组分和含量有一定差异，在麝香的甲醇和氯仿提取物中主要成分为脂肪酸类的烷烃类、醇类、酮类、醛类、固醇类、甾类及芳香族类物质，与其他文献报道[2]一致。脂肪酸衍生类、环酮类、甾类和固醇类为麝香的主要成分，而FASN、HMGCR和CYP11A1是合成上述物质的关键基因。

FASN基因对于长链脂肪酸的合成和活化过程具有重要作用[5,9,20-22]，如碳链的延伸、末端酰基的活化和双键的还原等，特别是与麝香酮(3-甲基环十五酮)碳链骨架的形成具有较强的关联性；同时FASN 合成的亚基酶对碳链的转运、去饱和、脱羟基和碳链活化等具有重要作用，如酰基载体蛋白β-酮酰ACP 还原酶和β-酮酰ACP水解酶等。对于固醇类和甾类物质而言，HMGCR是其合成的限速酶[8-9]，它参与甲羟戊酸的合成，而甲羟戊酸作为基础物质进一步形成异戊烯基焦磷酸、鲨烯、羊毛固醇和胆固醇等产物，对于类固醇物质的合成和分泌具有重要作用。CYP11A1是胆固醇的侧链裂解的关键酶[4,8,23]，对于孕烯酮、雄烯二酮和脱氢雄甾酮等物质合成至关重要，而在麝香检测到大量的雄性激素前体物质，如反式-脱氢雄甾酮，其合成同时需要HMGCR和CYP11A1的参与。因此，FASN、HMGCR和CYP11A1对于麝香的生物合成有较强的关联性，具有较大的研究价值和潜力。miRNAs 是基因转录后调控的关键因子[10-11]，对生物代谢组分的调控具有关键作用。1 个基因可受多条miRNAs调控，同时1 条miRNAs 可调控多个基因的表达，预测共同靶向林麝多个麝香合成相关基因表达的miRNAs 对探究麝香生物合成机制具有重要意义，同时也为牛和羊的肉质和繁殖性状相关研究提供参考。

林麝、牛、山羊和绵羊的FASN、HMGCR和CYP11A1基因3′-UTR 序列差异较小，同源物种之间基因的3′-UTR 有保守性，保守区域可能为共同的miRNAs 作用位点，同时也验证了序列拼接的正确性[12]；较高的自由能及更多的单链空泡结构将增加miRNA 与mRNA的结合概率[14-15]。通过SWChen 系统分别对4 个物种中调控FASN、HMGCR和CYP11A1基因的miRNAs 进行预测，在林麝和牛中，FASN基因受控的miRNAs 在各自miRNA 数据库中的占比相对山羊和绵羊较高，暗示林麝和牛FASN基因受miRNA 调控的潜力强于山羊和绵羊，与RNA 结构分析的结果一致；林麝HMGCR基因受miRNA 调控的数量占比较高，可能与其3′-UTR 长度和序列复杂程度有关；而CYP11A1基因在4 个物种之间预测的候选miRNA 数量差异不大，且RNA 二级结构空泡结构及最小自由能差距不明显，暗示4 个物种CYP11A1基因受miRNA 调控的潜力一致。候选miRNAs 数量在相应物种总体数据库中的占比和基因3′-UTR的RNA 二级结构和最小自由能预测，可为分析miRNA 对基因的调控潜能提供参考。

集合论分析显示：4 个物种FASN基因共同受到miR-24-3p 和miR-30b-3p 调控，miR-24-3p 可通过调控锌指同源基因ZHX2(zinc fingers and homeoboxes 2)影响FASN基因的表达[24]，miR-30c-5p 可调控下调FASN的表达改善肝性脂肪沉积[22]，预测miR-24-3p 和miR-30b-3p 有很大的潜力调控4 个物种的FASN基因表达。预测调控HMGCR基因表达的miRNAs 结果显示：miR-29b、miR-146a 和miR-2284d 为4 个物种miRNA的交集部分，miR-29 能够抑制小鼠HMGCR的表达而影响在肝脏的胆固醇沉积[25]，miR-29b 作为miR-29的同源miRNAs，具有较大的潜力调控4 个物种HMGCR的表达，miR-146a 与免疫炎症、心血管疾病、糖尿病及总胆固醇含量等具有一定的相关性[26-29]，这些因素均与固醇类和甾类物质合成相关，miR-146a 存在很大的潜能调控HMGCR基因的表达影响固醇类物质的合成导致相应疾病发生；通过双阈值分析可知：调控HMGCR基因表达的miRNA 有较大可能为miR-181b-3p 和miR-429 家族，miR-181d-5p 能够显著下调猪的HMGCR基因[30]，而鸡肝脏中的miR-429 与HMGCR的表达具有较强的相关性[31]，因此 miR-29b、miR-146a、miR-181b-3p 和 miR-429 都具有调控HMGCR表达的潜力。调控CYP11A1基因表达的miRNAs 集合分析显示均属于let-7 家族成员，let-7 家族是细胞中含量最为丰富的miRNA 家族，参与动物免疫、繁殖和生长发育等过程[32-34]；CYP11A1属于雄性激素的合成限速酶[8-9]，在各物种中均有较高的保守性。结合UTR 进化树、二级结构和miRNA 占比可知：4 个物种CYP11A1基因的相关参数基本接近一致，暗示let-7 对保守性较强的CYP11A1有较强的调控潜能；进一步的双阈值预测结果显示：miR-532-3p 也能够调控CYP11A1基因的表达，而miR-532-3p 被报道涉及女性不孕不育[35]，miR-532-3p 可能通过调控CYP11A1基因导致相关生殖类疾病。而在林麝中未发现共同调控FASN、HMRCR和CYP11A1的miRNAs 分子，但miR-205 和miR-143 家族存在交集，有参与麝香合成调控通路的可能。研究表明：miR-205主要参与AKT 和Wnt 信号通路[36-40]，miR-143主要参与NF-kappaB 信号通路[41-43]，暗示麝香的生物合成机可能涉及miR-205 和miR-143 参与的AKT、Wnt 和NF-kappaB 信号通路。

4 结论

本研究通过自主开发的SWChen 系统预测发现：调控林麝FASN基因表达的miRNA 包括miR-24-3p 和miR-30b-3p，调控林麝HMGCR基因表达的miRNA 有miR-29b、miR-146a、miR-181b-3p 和miR-429 家族，调控林麝CYP11A1基因表达的miRNA 有miR-532-3p 和let-7 家族，共同调控林麝3 个基因表达的miRNAs 主要为miR-205 家族和miR-143 家族。研究结果可为研究林麝、牛、山羊和绵羊FASN、HMGCR和CYP11A1的转录后调控提供参考，同时对解析林麝的麝香合成机制具有重要意义。