|
编译：卡德加，编辑：十九、江舜尧。 原创微文，欢迎转发转载。 棉花是一种重要的天然纤维作物，然而其基因图谱缺乏综合性和高分辨率。在本文中，研究者整合了四种互补的高通量技术，包括Pacbio long readIso-seq、链特异RNA-seq、CAGE-seq和PolyA-seq，系统地探索棉花不同组织或不同器官类型的转录景观。研究者设计了一个名为IGIA的计算流程，从整合的数据中重建精确的基因结构。研究结果揭示了棉花动态多样的转录图谱：组织特异性基因表达、TSSs和聚腺苷酸化位点的可变使用、可变剪接热点和转录读通(read-through)。这些被调控的事件在各个方面影响着许多基因，如功能RNA基序和蛋白质结构域的获得或丢失、DNA结合活性的微调(fine-tuning)以及同一复合体或通路中基因的协同调控等。这些方法和发现为进一步的功能基因组研究提供了宝贵的资源，如了解植物群落的天然SNP变异。 论文ID 原名：Multi-strategicRNA-seq analysis reveals a high-resolution transcriptional landscape in cotton 译名：棉花多策略RNA-seq揭示了高分辨率的转录组景观 期刊：Nature Communications IF：11.878 发表时间：2019.10 通讯作者：朱玉贤 通讯作者单位：武汉大学生物科学学院 DOI号：10.1038/s41467-019-12575-x 实验设计 结果 多策略RNA-seq用于高分辨率RNA景观 在之前的研究中，研究者已组装出亚洲棉基因组。为系统地识别完整和准确的基因模型，包括5 '和3 '的UTRs和捕获剪接isoform，研究者从花药、柱头、花瓣、苞片、萼片和全花等16个组织中提取RNA，在开花期后0天进行测序(DPA)，韧皮部、叶片、幼苗根、幼苗茎、子叶、种子、胚珠则为四个发育阶段(0、5、10、20 DPA)。Pacbio测序共得到92.8 Gb数据，ssRNA-seq测序数据为455.1Gb，CAGE-seq平均测序量为每个样83.1 million cleanreads，PolyA-seq平均测序量为68.8 million clean reads。 为了利用每种高通量技术的优势，研究者开发了一个集成分析流程来构建转录模型。对于支持Pacbio长read的区域，研究者使用了开发的IGIA(整合基因)工具，对于没有长read的区域，采用了TACO流程(图1a)。在IGIA方法中，包括TSSs、TESs和内含子在内的基因元件，分别从CAGE-seq、PolyA-seq和ssRNA-seq序列中鉴定。然后，根据所识别的元件，对打磨过(polish)的长片段进行分割。研究者在短read的基础上对每个元件的每次长read进行了验证，并在尽可能的情况下纠正了错误并完成了部分isoforms，并将它们划分为不同的isoforms类别(图1b)。在IGIA棉花基因中，56.7%只有一个转录本，17101个基因有两个或两个以上的转录本(图1c)。研究者的IGIA基因组与之前CGP和Cottongen注释结果有较大比例的重叠，并包含超过3400个新基因(图1d)。对于未在IGIA中标注的基因，它们的表达极低或长度过短(图1e)。 为了评估IGIA中的5'和3'注释，研究者计算了不同流程中注释的TSSs和TESs的距离，其峰值分别由CAGE-seq和PolyA-seq标识(图1f)。由于IGIA利用了CAGE-seq和PolyA-seq的信息，因此IGIA基因集的偏差相对于其峰值值最小。 接下来，研究者通过比较四组junction位点(JSs)来评估本研究中注释的剪接连接，包括来自ToFU、TACO、CGP和Cottongen的注释。在所有JS结果中，IGIA的特异性JSs数量最少(图1g)，表明较少的假阳性。基于IGIA基因结构注释，研究者对多个基因组特征进行了表征和比较，如亚洲棉中的UTRs、CDS、外显子和内含子，以及其他七个蛋白质编码基因的基因组。与动物基因组相比，亚洲棉的总体长度特征更接近于其他植物。然而，亚洲棉基因组仍具有独特的特征，包括比拟南芥和水稻更长的5'UTR、更宽的5'和3'UTR长度范围，以及比拟南芥更长的内含子。有趣的是，与其他两种棉花相比，亚洲棉基因的5'和3'UTR长度范围更广(图1h)，尽管这种差异可能是由于其他基因组的5'UTR和3'UTR注释不完整造成的。 为了探索这些功能元件的自然序列变化，研究者分析了来自230个亚洲棉家系的17,883,108个SNPs在IGIA基因中的分布，以及包括5'UTR、CDS、3'UTR和内含子在内的不同种类的元件。与预期一样，基因区域的SNP突变明显少于基因间区域，而CDS区域的自然变异频率最低(1.65 SNP/kb)。值得注意的是，TSSs和TES周围观察到SNP密度的急剧转变，在基因水平，外显子-内含子边界，以及UTR-CDS在mRNA水平的边界(图1i)。TSS上游30-40 nt观察到明显的山谷状分布(图1i上)，其中TATA box基序得到了富集。这些结果反映了IGIA对TSSs、TESs和JSs拼接的准确性。 图1 针对亚洲棉高分辨率RNA景观的多策略RNA-seq。A：整合基因亚型注释(IGIA)的实验设计和分析工作流程。B：IGIA识别精确isoform策略示意图。C：每个基因isoform数目的分布。D：CGP、Cottongen和IGIA基因注释比较的维恩图。E：FPKM分布，基因长度和低表达基因在基因亚群中的数量。F：TSS(左侧，CAGE-seq)和TES(右侧，PolyA-seq)与IGIA和其他方法组装的TSS(左侧，CAGE-seq)和TES(右侧，PolyA-seq)的峰值差异。G：仅由ToFU、CGP、TACO、Cottongen和IGIA支持的唯一剪切junction数。H：IGIA注释中棉花5'UTR、CDS、3'UTR、外显子和内含子的长度分布与其他7个品种的比较。I：SNP在IGIA基因复合体(上)和外显子(下)上的分布。 多种TSSs和可变启动子使用 基于改进的基因注释结果，研究者试图为亚洲棉建立时空动态转录组图谱。研究者使用ssRNA-seq数据量化了16个组织的基因表达，包括在两个生物重复中处于四个发育阶段(0、5、10和20 DPA)的胚珠。根据表达相关矩阵进行三组组织之间的比较分析：柱头、花药、花瓣三种组织；胚珠和种子四个发育阶段；其他营养组织和整朵花。这些结果与预期非常吻合：胚珠在不同发育阶段的基因表达最为相似，而营养组织和生殖组织的基因表达存在显著差异，因此被分别分组。 CAGE-seq能严格地选择5 '完整的RNA分子，用于在单碱基分辨率下识别TSSs。基于TPM > 0.5的结果，研究者确定44,728个TSS簇，对应22,863个基因位点，其中38.4%有两个或两个以上的TSS聚类簇(图2a)。为了估计CAGE-seq数据的准确性，研究者将簇分辨率定义为占簇中总信号80%以上的最小间隔。研究者评估了所有16个组织的CAGE信号的累积分布(图2b)，发现14.5%的TSS簇具有单碱基分辨率，71.6%的TSS簇具有10碱基分辨率。为了验证TSSs，研究者随机选择了55个基因，在16个组织的混合RNA样本中进行5 ' RACE试验。PCR结果显示47个基因(阳性率为85.5%)具有预期的条带，显示了较高的CAGE测序质量。 在一个基因位点上的多个TSSs将通过使用可变的TSS创造重新连接转录的潜力。接下来，研究者研究了全基因组TSS在棉花中的使用模式。将一个基因中的多TSSs分为三种类型：远端TSS (DT)、近端TSS (PT)和编码区TSS (CT)(图2c左)。研究者发现DT更容易被选择，表达水平明显高于PT和CT。 研究者的数据显示，选择性启动子使用可分别改变5888和2800个基因的5’UTR或蛋白N端长度(图2d上)。在可变5’UTR区域，检测到潜在功能RNA元素的富集，包括8387个miRNA结合位点，6377个U-rich元素，上游开放阅读框(uORFs) 5370个，RNA G-quadruplexes129个(RG4)，每一种都可能影响mRNA的稳定性或翻译效率(图2d下)。miRNA结合、U-rich和uORF对识别的元件的贡献最大，表明这些5’UTR的可变区域可能在基因表达的转录后调控中具有重要的调控作用。 接下来，研究者分析了2800个具有可变TSSs的编码蛋白的影响。大多数可变TSSs在不引起移码突变的情况下，会导致长度不超过500个氨基酸的末端截断(图2d)，这将导致基于预测的N端亚细胞定位信号或蛋白域的丢失。进一步分析了16个组织中TSS的使用动态，发现了6207个在组织间具有启动子开关的基因。可变TSS可能控制了参与的基因的结构域的存在或丢失、ABA信号转导(蛋白磷酸酶PP2C, Ga06g00624)，生长素转运(ABC转运蛋白，Ga14g01845)，硝酸盐转运(NRT1.2, Ga12g01164)，以及膜信号转导(富亮氨酸重复受体样蛋白激酶，Ga11g00886)。图2e是其中的一个典型例子，Ga12g01164，它是AtNRT1.2编码硝酸盐转运蛋白的同源基因。ssRNA-seq和CAGE-seq与0和5 DPA胚珠主要利用CT TSS表达短mRNA NRT-S的可能性一致，而晚期胚珠(10和20 DPA)等组织则转向DT TSS表达较长的mRNA NRT-L(图2e)。 研究者的5'RACE检测结果证实了TSS开关的现象(图2)。NRT1.2有12个跨膜(TM)结构域，但由于TSS从远端向近端变化，NRT-S编码蛋白N端截断的基础上，失去4个TM结构域，可能导致TM簇的蛋白构象由紧致状态变为松弛状态(图2g)。比较它们在NO3上的活性，在HEK-293细胞中，研究者对NRT-L和NRT-S进行了异位表达，测试同位素15NO3-的吸收。2和8小时的测定结果进一步证实了它们在NO3活性上的差异，NRT-S的吸收活性明显低于NRT-L(图2h)。这些结果表明，在发育早期(0和5 DPA)，胚珠可能具有独特的氮转运和代谢，不同于其他组织。上述结果表明，可变TSSs是棉花mRNA的一个共同特征，它通常在mRNA或蛋白中产生选择性的N端来调节组织的规范和发育。 图2 亚洲棉多重转录起始和可变启动子的使用。A：每个基因转录起始位点(TSS)的统计。B：CAGE-seq高分辨率识别TSS。C：TSS使用信号在多个TSSs基因中的统计，包括远端(DT)、近端(PT)和编码TSS(CT)。D：由可变TSS引起的5 '端变化的统计和CDS长度变化的分布。5'UTR区域的潜在功能RNA元件，包括RBP结合、U-rich、RG4、二级结构和uORFs。E：一个长(L)和短(S)isoform基因的例子，由于可变TSS的出现(箭头)。F:TSS引起变化的5'RACE验证。G：长、短isoform的三维蛋白质结构预测。H：15NO3-吸收NRT-L、NRT-S和空载体细胞系之间的活性比较。 发育调节的TES选择 选择性TES选择或APA产生的mRNA具有不同的3 '端，其编码序列或3 ' UTR各不相同，通过调节mRNA的稳定性、定位和翻译效率来增加转录组的复杂性。在棉花方面，APA还没有得到系统的研究，基于PolyA-seq的3 ' end信息，研究者对16个组织中所有表达的基因进行了全基因组TESs分析。从比对到棉花基因组的1.1×109个PolyA-seq reads中，共产生了70,502个TES簇，其中43,237个TES簇的平均TPM > 0.5，对应23,736个基因位点。在棉花所有表达基因中，40.2%至少有2个TESs(图3a)，平均每个基因有1.56个polyA位点，低于拟南芥的74.9%和大米中的47.9%。与上述TSS分析类似，研究者接下来评估了所有16个组织的polyA位点分辨率的累积分布，结果显示，36.4%的TES集群具有单碱基分辨率，76.9%的TES集群具有10碱基分辨率(图3b)。 研究者分析了TES周围序列的核苷酸(nt)组成，在PolyA位点发现了典型的A/U富集模式，与之前报道的拟南芥和水稻的结果相似，但在对照序列中没有发现(图3c)。这些结果表明，研究者的TES注释是可靠的，并且具有较高的分辨率。 为了研究APA在棉花组织中的动态变化，研究者从PolyA-seq数据中量化了TES用量并从10681个基因中鉴定出APA。为了区分TES的使用模式，研究者将单个基因中的多个TESs分为三种类型：远端TES (DT)、近端TES(PT)和编码区(CT)的TES。与转录起始的情况不同，远端TESs (DT)的使用频率低于PT，但DT和PT的表达水平远高于CT TESs(图3 d)。这与之前在水稻中的研究高度相似，但与小鼠的研究不同，在小鼠中，远端TESs的使用频率更高，表达水平高于PT。这一有趣的现象表明，植物和动物对选择性TES的偏好不同。 3 '端变异大部分(87.5%)局限于3 'UTR，只有8.3%的变异发生在编码区域和3 'UTR之间；3.4%会通过打破ORF框将编码RNA改变为lncRNA(图3e)。为了进一步分析改变3'UTR长度的影响，研究者分析了可变3'UTR区域的组成，发现了RNA功能基序的富集(图3f)，包括miRNA结合位点、U-rich元件和RG4。 接下来，研究中分析了3 '端缩短或延长的模式。UTRs通过16个组织，发现平均长度为3'UTRs在不同组织间差异较大(图3g)。为了从全基因组数据中深入了解APA的差异，研究者鉴定了具有3 'UTR的基因在两两组织之间的变换模式(图3h)。花药、柱头和花瓣(组织ID 1-3)相对较少，种子和胚珠(组织ID 4-8)的APA转换事件比营养组织(组织ID 9 - 11,14 - 16)的数目更多。大多数棉花胚珠和种子的APA交换可能在纤维发育和种子成熟过程中起重要作用。 在胚珠(20dpa)和种子中表达的1,899个mRNA中，有189个基因表现出APA转换，且差异显著(图3i)。在这些情况下，更多的基因使用种子中的远端polyA位点，如Ga07g00163，编码组氨酸激酶，这可能在跨细胞膜的信号转导中发挥作用。总的来说，根据结果，在棉花中建立了一个高分辨率PolyA位点景观，揭示了它们的序列特征、发育过程中的动态规律和组织规范。 图3亚洲棉中多重转录终止和可变3'UTR的使用。A：每个基因的转录末端位点(TES)的统计。B：PolyA-seq中TES簇信号的碱基分辨率的累积曲线。C：在TES周围的核苷酸组成。D：多TES基因中TES使用情况的统计，包括远端TES(DT)、近端TES(PT)和编码TES(CT)。E：对RNA影响的统计(饼图)和由于使用可变TES而导致CDS长度变化的分布。F：潜在的功能RNA元件，包括RBP结合，U-rich和RG4基序。G：组织间平均3'UTR长度。H：在两两组织间有使用可变TES的基因数量。I：APA在胚珠(20dpa)和种子(左)之间发生变化的基因，以及两个有代表性的基因例子(右)。 棉花动态剪接开关和微外显子 可变剪接分析已在雷蒙德棉(DD)和陆地棉(AADD)中开展过，但在亚洲棉(AA)基因组中还未有报道。基于IGIA的注释结果，研究者对亚洲棉的23,451个多外显子基因进行了系统分析。总的来说，在42.1%的多外显子基因中鉴定出23,756个事件，在AS基因中约有2.4个事件(图4a)。 结果显示内含子保留(RI)在亚洲棉中占大多数(62.2%)，在所有报道的植物中是最高的。这一结果也在另外两种棉属植物中被发现。其他三种AS事件模式(A3SS、A5SS和SE)的百分比分别为18.8%、11.3%和7.7%，如图4a所示。研究者进一步描述了AS事件对蛋白编码序列的影响，特别是对各种功能域的影响，发现AS对棉花Pfam和CDD蛋白结构域的完整性有很大影响，对具有无序区域、短线性基序(Slim)和跨膜(TM)特征的序列有很大影响(图4b)。随机选择不同AS事件的基因，用RT-PCR检测AS注释的准确性。ssRNA-seq在不同组织间的高动态剪接开关与RT-PCR结果一致。基于ssRNA-seq数据，对所有IGIA基因的AS定量显示了在组织内的全面动态概况。 微外显子(MEs)是一类特殊的外显子，其长度通常只有3 nt。本文对亚洲棉的MEs进行了总结。与人类结果相似，内外显子(internal exon)的长度分布在50 nt左右急剧下降(图4c)，研究者定义了MEs的长度阈值为51 nt，共从IGIA转录本中鉴定出7784个MEs。图4d显示Ga05g01029中的45 nt ME编码的乙烯应答转录因子RAP2-7蛋白，是一种AP2域的转录因子。这ME编码β-sheet的一部分，α-螺旋的AP2域。与裸子植物卷柏和水稻、拟南芥的直系同源结构比较，这个ME的长度和序列高度保守，尽管这个基因的其它外显子在不同基因组中差异很大(图4d)。RT-PCR检测结果显示，ME的AS在棉花组织中具有明显的组织特异性，这从组织间的外显子包合比例可以看出(图4e)。图4f显示在相同的蛋白浓度下，AP2-I可以同时与DNA探针结合，AP2-S则几乎没有结合，说明转录因子RAP2-7可能通过调控微外显子的剪接，改变其在不同组织中与基因组的结合活性，这是动物和植物共有的一种机制。 图4 亚洲棉中可变剪接调控和热点。A：利用IGIA注释对16个组织中4类AS事件及其相关基因进行统计。B：AS事件对蛋白质功能域的影响。C：内部外显子的大小分布。D：一个位于AP2 DNA结合域内的微外显子(ME)，在植物中保守。E：使用RT-PCR进行ME验证。F：EMSA分析显示，两种isoform w/o ME、I和S具有不同的DNA结合亲和力。G:具有代表性的AS热点(浅绿色突出显示)和RT-PCR验证的基因来自“冷”区域(CR)和热点(HS)。H：具有代表性的动植物基因组中的AS热点模拟。I：来自HSs和CRs的外显子和内含子的长度比较。 此外，某些基因的几个区域显示出高度富集的AS事件，称之为AS热点(图4g，亮绿色框)。基于AS热点的统计模型，研究者对棉花基因组进行了严格的扫描，发现了135个含有至少一个AS热点的基因。GO分析表明，这些基因的分子功能主要富集于ATP和RNA结合蛋白及蛋白激酶。进一步分析AS热点的蛋白特征，发现热点主要影响一个保守的蛋白结构域。例如，在蛋白激酶类shkb蛋白中，AS事件往往集中在编码MAPK激酶域的热点上。利用RT-PCR进一步确认该基因在“冷”区和热点的不同频率的AS(图4g，右)。利用相同的计算模型，研究者发现热点也存在于其他植物中，如水稻和拟南芥，但不存在于人类和小鼠中(图4h)。这表明AS热点可能是一种植物特异性剪接现象，可能与植物特异性顺式元件和反式因子有关。有趣的是，通过对热点地区和寒冷地区基因结构的比较，研究者发现炎热地区比寒冷地区的外显子更短，而内含子则不一样(图4i)。 综上所述，研究者的分析提供了棉铃虫选择性剪接景观的综合信息，对阐明亚洲棉选择性剪接的调控机制和作用具有重要价值。 多顺反子的发现及其基因组特征 多顺反子转录(Polycistronic transcription)是由多个ORF或单个启动子基因共同转录形成的多顺反子mRNA，在原核生物和真菌中普遍存在，在真核生物中少见。在植物中，一个罕见的例子是番茄，它是谷氨酰胺激酶和谷氨酰胺磷酸还原酶位点的双顺反子单位，类似于原核多顺反子操纵子。对亚洲棉全长转录本的IGIA注释分析表明，亚洲棉基因组中许多相邻的基因座具有跨两个或多个ORF的多顺反子转录，这些ORF在亚洲棉和CGP中被注释为独立的单顺反子基因(图5a)。 RT-PCR和Sanger测序显示，多顺反子转录本确实有两个长度为4.2和2.8kb的亚型，由于另一个内含子保留事件。TSS和来自混合组织的TES信号表明，这三个基因可能有自己的转录单位，与多顺反子表现出不同的表达模式，这得到了ssRNA-seq和RT-PCR结果的支持(图5b)。多顺反子的组织特异性表达提示它可能在组织特异性信号的控制下被激活。在严格的筛选标准下，在亚洲棉中共鉴定出1115个多顺反子(1081个双顺反子和34个三顺反子)(图5c)，占本研究所有注释位点的5%。 为了了解多顺反子的潜在机制和功能，研究者分析了这些多顺反子和控制基因的几个基因组特征。多顺反子内的ORFs(第一个ORF的终止密码子与第二个ORF的开始密码子之间的距离)明显小于单顺反子基因之间的距离(图5d)。来自同一多顺反子的基因在组织间的共表达相关性比随机对照更强(图5e)。此外，在相同的多顺反子转录本中的基因往往具有相似的GO注释，位于相同的KEGG通路，并相互作用(图5f-h)，这可能是由于它们的协同转录具有更好的协调性。为了研究多顺反子的起源，研究者将多顺反子与其他植物基因组进行共线性分析。结果表明，在大量植物中，许多多顺反子与基因位点存在共线性关系。多顺反子共线性存在于40多种物种中，多顺反子的保护表明多顺反子可能起着重要的调节作用(图5i)。 图5亚洲棉中多顺反子的鉴定及基因组特征。A：多顺反子的一个例子，Pacbio长read支持的三个基因。多顺反子中的基因可以独立转录，具有独立的TSSs和TESs。检测多顺反子表达的RT-PCR探针标记在IGIA转录本下面。B：多顺反子转录的RT-PCR验证，有四对探针。C：多顺反子中CDS对之间距离的统计与IGIA基因中CDS与其最近的CDS之间距离的比较。D：利用IGIA基因鉴定了具有两个和三个CDSs的多顺反子的数目。E：多顺反子中CDS对与随机CDS对的表达相关性。F-H：GO相似度的比较。F：相同KEGG通路的比例。G：可能的蛋白-蛋白互作。H：在多顺反子的CDS对和随机的CDS对之间。I：棉花多顺反子在54个植物基因组中的共线分析。 讨论 研究者精确绘制了亚洲棉TSSs、TESs、剪接位点和完整的isoform收集的基因图谱，以及它们在多个组织中的动态图谱。本研究综合了多种RNA-seq技术，将为棉花和植物研究者提供重要的参考。 新兴的RNA测序技术需要不断更新分析工具，进一步提高转录组的注释准确性。受GRIT算法的启发，IGIA将CAGE-seq和PolyA-seq进行了整合，分别确定了基因的起始和终止边界。对于基因区域，GRIT仅使用来自RNA-seq的短读数据来进行算法重建isoform结构，而IGIA则使用Pacbio长read为isoform骨架，短read使用RNA-seq用于纠正长read时的局部错误，从而避免了短read时易出错的重构。综合四种模式的RNA-seq数据，IGIA在本研究中得到了多种实验验证的支持，在基因注释方面表现良好。在本研究中，IGIA倾向于检测TSSs或TESs上具有特异性差异且具有可靠表达信号的isoform。IGIA要求对剪接连接提供多种证据，目的是识别可靠的基因结构，同时对isoform进行潜在估计。IGIA将进一步优化，以减少计算时间和内存需求，提高长read时的纠错能力，提高其灵敏度。 在以往的全基因组关联研究(GWAS)中，棉花中已经发现了许多与农艺性状相关的位点，但非CDS多态性(non-CDS polymorphisms)很难解释或优先研究。提供功能元件的完整目录，包括最新的5 ' UTRs和3 ' UTRs，研究人员可以提出新的假设，进一步分析CDS之外的重要序列变异。在本研究中，研究者整合了亚洲棉AA基因组上的932个GWAS位点，并观察到许多GWAS位点定位于调控区域，包括转录起始和终止位点、选择性5' /3' UTRs、组成型或可变剪接(AS)位点。这为进一步从这些GWAS位点解码棉花纤维长度、品质等重要农艺性状的分子机制提供了重要的基础。 据报道，人类启动子使用异常会导致多种疾病，包括癌症。然而，在植物中，只有在玉米和拟南芥中有有限的关于确定精确TSSs的报道。在本研究中，研究者为26,763个基因位点鉴定了91,707个TSS簇，其中近40%具有多个启动子，表明棉花基因组中广泛使用了可变启动子。在本研究中，2,800个基因具有不同的蛋白质isoform(CDS变化)是由亚洲棉的可变TSS选择引起(图2d)。此外，可变TSS对蛋白长度的影响可能不仅限于N端亚细胞靶向信号序列，还可能改变蛋白结构和功能(如硝酸盐转运体中TM结构的缺失，NRT1.2可能会减弱NO3-的摄取活性)(图2f-h)。这些结果拓宽了对植物可变TSS的认识，为植物学家提供了新的研究视角。 更多推荐 1 科研 | PNAS：转录组学揭示急性和慢性饮酒对肝脏昼夜新陈代谢有不同的影响
|
|
|
