2022年7月，刘仁义团队（福建农林大学）和张兴坦团队（福建农林大学/中国农科院深圳农业基因组研究所）在Horticulture Research期刊（IF=7.291）发表了题为Pan-transcriptome assembly combined with multiple association analysis provides new insights into the regulatory network of specialized metabolites in the tea plant Camellia sinensis 的研究论文，这篇论文很有意思，不仅涉及宏转录组组装，还涉及GWAS、eQTL分析，是一篇很好的数据再利用的范文！今天我们一起来看一下，文章都有哪些内容。
茶树中的特异性代谢物对植物调节胁迫抗性和风味形成至关重要，而当前研究对其合成和调节途径并不完全清晰。本研究首先组装了茶树宏转录组，其次鉴定了大叶茶（CSA）纯系和小叶茶（CSS）纯系之间的遗传差异，并探究了代谢与遗传之间的关系，最后揭示了次级代谢物调控网络及其关键基因。
文章技术路线
1. 数据获取：从NCBI下载136份茶树RNA-seq数据；从EBI下载对应的136份茶树的代谢组数据，以上数据源于同一项目提交的转录组和代谢组数据。
2. SNP分析、构建进化树、群体结构分析。
3. 构建泛转录组、转录本表达定量、获取基因序列。
4. 鉴定CSA纯系和CSS纯系之间的差异表达基因、ePAVs。
5. GWAS分析代谢物相关SNP、线性回归模型鉴定代谢物相关QTTS、鉴定候选基因、鉴定eQTLs。
文章内容
1. 泛转录组拼接
利用聚类分析将转录组数据分为3个群体（Group1-3），并基于有参基因组拼接和从头拼接的办法对每个群体进行群体水平的转录本拼接，用于泛转录组拼接：
①下载RNA-seq数据并与“铁观音”基因组比对后，将筛选到的SNP用于进化分析和群体结构分析。结果表明，3个群体中的CSA纯系和CSS纯系均有聚类趋势。由于茶树存在自交不亲和现象，因此其余的品种可能是杂交种或与近缘种发生了基因渗入（图1A）。基于参考基因组拼接并组装了每个群体的转录组（图1B），同时过滤掉低表达基因。
②未比对上基因组的reads经从头拼接后，被鉴定为新基因。
至此，泛转录组拼接完成。
图1 茶品种的遗传分化。A 136个茶品种的进化分析。B Reads在参考基因组上的比对率
2. 基因表达与功能差异分析
分析发现共有21891个基因为核心基因（在至少95%的品种中表达），896个基因为品种特异基因（图2A），并且核心基因和特异基因主要参与的KEGG通路不同。此外，高表达基因和中等表达基因主要富集在植物基本生物学通路上，而低表达基因主要富集在胁迫适应途径上，这或许是植物的自我保护机制。
由于相关性和PCA分析结果不理想（图2C，D），作者便检测CSA和CSS之间的DEGs，发现这些DEGs可能决定物种分化。CSS的上调基因主要参与次级代谢物相关途径，而CSA的上调基因参与的生物学过程更广。
CSA和CSS分别存在210和154个特异表达基因（ePAVs），并且二者的参与的主要生物学过程存在差异，这与CSA和CSS的表型差异密切相关。
图2 群体间的基因表达差异。A基因表达分布。B不同表达水平的基因比例。C皮尔森相关性。D主成分分析
3. 次级代谢物积累差异
茶氨酸、可可碱、大多数原花青素类物质、C、EC主要积累在CSA，而ECG、EGCG主要积累在CSS（图3A，B）。此外，二者特异性积累的其他代谢物种类也有所差异。
图3 CSA和CSS的代谢差异。A氨基酸和多肽、生物碱、黄酮醇、原花青素是差异标志物。B儿茶素（C）、表儿茶素（EC）、表没食子儿茶素（ECG）、表没食子儿茶素没食子酸酯（EGCG）的差异
4. 次级代谢调节网络和候选基因
综合GWAS和QTLs分析结果，鉴定到78个候选基因（图4A，B）。在泛转录组水平上鉴定的Cis-eQTL和Trans-eQTL在基因组上的分布不均匀（图4C，D）。从eQTL-代谢物网络中鉴定到候选转录因子，如AP2/ERF-ERF、MYB、WD40、bHLH、C3H等。
图4 次级代谢物调节网络。A 3种GWAS分析结果。B QTTS和GWAS分析所得候选基因。C MLM和FarmCPU分析所得eQTLs。D基因密度分布。基因密度、Cis-eQTL和Trans-eQTL在基因组上的分布情况（a-c）。E调节网络。网络包含代谢物、候选基因、eQTLs及候选基因+eQTLs
5. CsTGY14G0001296 可能决定原花青素积累差异
GWAS和QTTS均定位到基因CsTGY14G0001296与32种代谢物相关，但该基因表达趋势与原花青素类物质的积累趋势之间无显著正相关性（图5A，B）。功能注释信息显示该基因是无色花色素双加氧酶（ANS）基因，负责原花青素合成，表明该基因可能在形成紫叶品种富含原花青素品种过程中发挥重要作用。
图5 CsTGY14G0001296与32种代谢物的相关性系数（A）与热图
6. 黄酮醇差异积累的机制
CsF3’5’H（CsTGY11G0002074）在CSS中高表达，这可能是造成儿茶素类化合物在CSA和CSS间积累特异性不同的原因（图6A）。GWAS和QTTS分析、基因启动子分析结果也支持该推论（图6B）。此外，CsF3’5’H的表达还可能受NAC等TF影响。
图6 黄酮醇积累差异。A 黄酮类化合物合成通路。B GWAS分析CsF3’5’H为候选基因
文章亮点
这篇文章真是数据再利用的典范，下载转录组、代谢组数据进行泛转录组组装和GWAS、eQTL分析也颇具新颖性，值的一看！下面我简单总结下该文章的几处亮点，供您参考!
1. 本文系统研究了茶树特异性代谢方面的表型、代谢、遗传之间的联系，填补了园艺经济作物在特异性代谢的分子调节机制方面的空白。
2. 单个参考基因组不能代表该物种全部的基因多样性，而泛基因组组装成本较高，因此考虑用采用泛转录组组装的方式，既能降低成本，又兼具检测基因组变异和基因表达差异的能力。本文为首次在茶树中应用泛转录组测组装策略。
3. 由于QTL定位和GWAS分析不足以解释所有表型差异，而表型差异与DEGs、ePAVs等密切相关，因此本研究使用泛转录组分析揭示了DEGs和ePAVs在CSA和CSS分化中的作用。
4. CSA的原花青素含量高于CSS，可能是由CsANS造成的。儿茶素类化合物在CSA和CSS间特异积累的机制可能与CsF3’5’H的表达差异和多种TF调节相关。这些因素会影响茶树品质。
未来，通过定向杂交CSA和CSS并鉴定DEGs，可能有助于选育抗病性更强、风味更佳的品种。
网友评论