未来组项目文章 |异源四倍体棉花全长转录组测序分析

未来组参与的异源四倍体棉花全长转录组项目文章于9月11号见刊New Phytologist，华中农业大学研究团队基于PacBio测序，并整合Iso-Seq流程开发适用于区分2套亚基因组转录本的分析方法，克服了短读长测序在解析多倍体isoforms的技术瓶颈，揭示了纤维特异性的可变剪接事件，2套亚基因组中部分同源基因的isoforms差异，并在isoforms水平上揭示了miRNA对可变剪接事件的调控，为研究多倍体物种可变剪接提供了新的研究角度。本次研究采用PacBio SMRT测序技术，对异源四倍体棉花进行全长转录组测序分析，直接得到更为完善的全长isoforms，完善棉花转录组注释，解析了异源四倍体棉花组织和2套亚基因中AS复杂性，并整合多组学数据，揭示AS的调控机制。
已有大量研究基于高通量测序技术揭示了真核生物转录组中的可变剪接（AS）的广泛性和复杂性，以及在植物发育阶段或应激反应中对AS进行全基因组范围内研究。基于短读长的RNA-seq在准确拼接重构全长isoforms时充满挑战，无法避免假阳性AS事件，尤其是在面对多倍体物种时，情况更加复杂。

研究方法

1.采集异源四倍体棉花Gossypium barbadense L. cv 3-79的根、下胚轴、叶、花瓣、花粉和花柱6个组织样本，提取RNA，等量混合，反转录为全长cDNA，构建1-2kb，2-3kb，3-6kb文库，上机PacBio RSII，P6C4，共测15Cells。

2. 采集Gossypium barbadenseL. cv 3-79 6个发育时期（开花后0天即0 DPA、7 DPA、10 DPA、12DPA、20 DPA和30 DPA）的棉花纤维样本，提取RNA，等量混合，反转录为全长cDNA，构建1-2kb，2-3kb，3-6kb文库，上机PacBio RSII，P6C4，共测15Cells。

研究结果

1 全长转录组分析流程优化

在异源四倍体棉花全长转录组分析中，针对四倍体物种开发了整合性Iso-Seq数据分析流程(https://github.com/Nextomics/pipeline-for-isoseq)，其中包含了数据质控、转录本分类、isoforms聚类及转录组后续分析（Figure 1），使用Samtools phase对来源于2套序列相似性极高的亚基因组的转录本进行区分（Figure 2）。

Figure 1 异源四倍体棉花的Iso-Seq数据分析流程

Figure 2 At和Dt 亚基因组上的转录本比对

2 Iso-Seq解析了多倍体棉花广泛的可变剪接、可变多聚腺苷酸、融合基因、新LncRNA等形式

（1）下机数据分析，共得到全长转录本1,096,932（ca. 43.2%）（Figure 3A），经mapping、phasing、clustering、consensus后总共得到44,968个基因的176,849个isoforms，其中全长isoforms平均2,175bp，比参考序列转录本平均1,462bp的长度长（Figure 3C），并且Iso-Seq可得到更多的多外显子isoforms（Figure 3D）。

（2）通过与参考序列比对，更新了18,008个基因，确定了222个融合基因，在新的转录本中确定了2,447个LncRNA，与LncRNA_V1比较，确定了365个新的LncRNA（Figure 3E）。在Iso-Seq中检测的44,968个基因中，基因上平均polyA位点数目为2.82（Figure 3F），同时分析了polyA位点侧翼核苷酸序列，其表现出核苷酸偏向性特征，在3’UTR的polyA剪切点的上游富含U碱基，在下游富含A碱基（Figure 3G），同时在polyA剪切点上游确定了2个保守的motifs（AAUAAA和UGUA）。

Figure 3 经Iso-Seq得到的棉花转录组图示

1. 转录本分类 B.全长isoforms在文库中mapping汇总 C. 参考数据和Iso-Seq数据中转录本长度分布 D. 在参考数据和Iso-Seq数据中isoforms中外显子数目 E. Iso-Seq数据及中检测的全长LncRNA与LncRNA_v1的Venn图 F. 基因上polyA位点数目分布 G. polyA剪切点（-50 ~ +50）的核苷酸相对频率

（3）经Iso-Seq检测133,229个AS事件，发现63.8%基因的AS事件为内含子保留（IR）（Figure 4A），平均每个基因对应3.93个isoforms，大概是参考注释中的2.9倍（Figure 4B），另外，随机挑选了5个基因，根据其预测转录本设计引物，经RT-PCR来验证Iso-Seq中AS事件的检测，发现扩增片段与预测片段相符，同时，也发现了isoforms的表达量具组织特异性（Figure 4C）。

Figure 4 AS事件特征图示及全长isoforms的RT-PCR验证

A. AS事件分类及相对应的基因和AS事件数目

B. 参考数据和Iso-Seq数据中基因对应isoforms数据

C. 5个基因的AS事件的RT-PCR验证

3 组织特异性isoforms的结构和表达量差异分析

结合来自纤维样本和非纤维样本的Iso-Seq数据对比，检测到来自15,871个基因的66,652个两者共有AS事件（Figure 5A），经RNA-Seq数据对isoforms的表达量进行分析，发现组织特异性isoforms与isoforms组织特异性表达表现出一致（Figure 5B），挑选PB.1316进行验证，PB.1316可以转录为T1和T2两个isoforms，实验发现T1有些在非纤维组织中表达，如根和叶，在纤维发育过程中，T1表达量不断降低，而T2出现高表达（Figure 5C）。

Figure 5组织特异性AS事件和isoforms图示

A.组织特异性AS事件和相应基因的Venn图

B. 组织特异性AS事件类型图示

C. 最大组织特异性数值分布

D.PB.1316的2个isoforms（T1和T2）的转录模式 E. 对7个样本中PB.1316进行RT-PCR验证

4亚基因组的同源基因表现出isoforms结构差异

（1）结合已研究的16,077同源基因对和Iso-Seq数据，得到6,202对同源基因，可以分为3大类，其中group I中1,605对基因中At基因组isoforms数目呈收缩，group II中3,017对基因两者没有差异，而group III中1,580对基因中Dt基因组isoforms数目呈收缩，并对三个group进行GO分析（Figure 6A）。

（2）挑选来自不同group的三个基因对，group I的PAP10在At亚基因组中只转录1个isoforms，而在Dt亚基因组中转录6个不同isoforms；group II的ERD3在At和Dt亚基因组中都转录6个不同isoforms；group III的CPN60A在At基因组中可转录16个不同isoforms，而在Dt基因组中仅转录1个isoforms（Figure 6C-D）。

Figure 6  At和Dt亚基因组中同源基因isoforms数目

A. At和Dt亚基因组中同源基因isoforms log2比值及GO注释

B. At和Dt亚基因组同源基因PAP10的isoforms

C. At和Dt亚基因组同源基因ERD3的isoforms

D. At和Dt亚基因组同源基因CPN60A的isoforms

5 isforms多层面的调控——miRNAs和表观修饰

（1）对6个组织进行小RNA测序并整合已有的研究中测序的miRNAs，发现miRNAs结合AS在isoforms水平调控基因表达量。如PB.42410有5个isoforms，其中2个是miR399的靶序列； PB.18525可转录2个isoforms，其中1个isoform是miR397的靶序列，因为ES等事件PB.42410或PB.18525转录的其他isoforms而缺失了miRNA靶位点；PB.44799和PB.44311因为IR事件而获得miR7484和miR8634的靶位点；PB.27256和PB.2778转录的isoforms在 3’UTR和5’UTR分别获得miR396和miR827的靶位点（Figure 7A-F）。

Figure 7 miRNAs对全长isoforms靶向确定

（2）结合已有棉花纤维发育阶段的表观修饰分析数据与此次研究中的纤维样本的AS数据，在分析不同AS中核小体占位密度和DNA甲基化水平中，发现其可能都在定义外显子中起重要作用（Figure 8A和B）；基于全长isoforms的UTRs注释，将其与CDSs的DNA甲基化水平对比，发现CG，CHG，CHH甲基化水平差异较小（Figure 8C,D,E,F,G和H），另外，经RdDM通路参与DNA甲基化的24-ntsiRNA在UTRs区域明显增加，可能与3’UTRs区域的CHG和CHH甲基化水平增加相关，这些都给AS和基因表达的甲基化调控提供了新思路。

Figure 8 isoforms中核小体占位和DNA甲基化分析图示

本次研究中经Iso-Seq测序，拓展了多倍体棉花的转录组研究，并从组织特异性和亚基因组角度确定了其中可变剪接的复杂性，结合多组学数据，分析了可变剪接事件的调控机制，而这些研究的基础关键在于需先获得全长isoforms。同时，文章也给予我们在转录组研究一些启发，在表型差异的转录组研究中，不仅需关注与其相关基因的表达量相关，也要考虑到可变剪接中的不同isoforms作用。

参考文献：

Wang,M., Wang, P., Liang, F., Ye, Z., Li, J., Shen, C., Pei, L., Wang, F., Hu, J.,Tu, L., Lindsey, K., He, D. and Zhang, X. (2017),A global survey ofalternative splicing in allopolyploid cotton: landscape, complexity andregulation. New Phytol. doi:10.1111/nph.14762