选取云南省农业科学院茶叶研究所试验基地10年树龄的福鼎大白茶花(父本)、10年树龄的佛香2号花(母本)以及他们杂交后代5年树龄的不育花, 9月25号采摘它们的花芽、花蕾、花, 利用液氮迅速固样。

选用CTAB-LiC1方法^{[14, 15]}, 分别提取花芽、花蕾、花总RNA, 采用2%的琼脂糖凝胶电泳分析RNA降解程度以及是否有污染, 利用Nanodrop检测RNA的纯度、Qubit定量RNA浓度以及用Agilent 2100精确检测RNA的完整性, 将检测合格的总RNA样品送北京诺禾致源生物信息科技有限公司, 富集mRNA、合成双链cDNA和构建文库, 并采用Illumina HiSeq 2000平台测序。

利用RSEM软件的bowtie2将Trinity拼接得到的转录组作为参考序列(ref), 将每个样品的clean reads往ref上绘制mapping。采用TMM对read count数据进行标准化处理, 再利用DEGseq技术分别对父本花与母本花、父本花与子代不育花、母本花与子代不育花进行差异分析。按照FPKM法(Reads per kb per million reads)计算差异基因的表达量^[16]。筛选阈值< 0.005以及|log₂Fold_change| > 1; 对于差异基因, 如果其log₂ Fold_change > 0, 则认为是上调表达, 反之, 若log₂Fold_change < 0, 认为是下调表达。

根据研究目的, 不育基因重点筛选在父本花与子代不育花、母本花与子代不育花之间差异表达, 且在父本花与母本花之间没有的差异表达基因。在KOBAS2.0 (http://kobas.cbi.pku.edu.cn/home.do)网站进行DEGs的KEGG代谢途径富集分析。GO和KEGG统计检验中均采用Benjamini Hochberg方法进行多重比对检验。

从测序获得的差异表达基因中, 随机挑选5个进行实时荧光定量PCR测定, 采用CTAB-LiC1法^{[14, 15]}提取总RNA。按照Fermentas公司的RevertAid First Strand cDNA Synthesis Kit说明书合成cDNA。依据测序所得序列库对应的代表序列, 使用Primer 5.0软件设计实时荧光定量PCR引物, 以茶树的18S (AB120309.1)基因作为内参基因(表1)。实时荧光定量PCR的反应体系及程序参见Fermentas公司的SYBR Green master mix说明书, 设实验重复3个, 采用2^{-Δ Δ CT}法计算基因相对表达量。

表1

Table 1

表1(Table 1)

表1 实时定量PCR选用基因及引物 Table 1 Genes and primers list for qRT-PCR 编号 Code 基因 Gene 正向引物 Forward primer (5° -3° ) 反向引物 Reverse primer (5° -3° ) 1 c112358_g1 GCCAAGGTCTCCCTCATCAT AGTGTGTGCTCCATAGGTCG 2 c112358_g5 GGAGGAGGAGAACAAGCACC TGGAAATCTCGCACTCGCTG 3 c115764_g1 CCGTCATCACCTTCTCCAAC CACGGAGGCTCTGGAGGATG 4 c116710_g4 AAGTTCTGGCAAAGGGAAGC ATTGCCTGTGGCTTTCTTGC 5 c123710_g4 ATCGCTCTCATCGTCTTCTC ACACTACTGCCAGCAACTTC 13 18S ACATAGTAAGGATTGACAGA TAACGGAATTAACCAGACA

表1 实时定量PCR选用基因及引物 Table 1 Genes and primers list for qRT-PCR

对父本花(FB_H)、母本花(MB_H)、子代不育花(ZD_H)转录组测序, 经去除接头、N和低质量reads, 在父本花中获得58 121 920条clean reads, 在母本中获得73 150 088条clean reads, 在不育子代中获得68 216 824条clean reads。获得序列的Q30值均大于91%, GC含量均大于43% (表2), 符合下一步分析要求。再利用Trinity对clean reads进行拼接共获得272 804条转录本(Transcripts)序列, 取每条基因中最长的转录本作为Unigene, 筛选出204 199条Unigene, 总长为120 173 858 bp, 平均长度为588 bp。对转录本及Unigene的长度分别统计, 其中Unigene小于301 bp的有81 272条, 分布在301~501 bp的有56 681条, 大于2000 bp的有8719条(图1)。

表2

Table 2

表2(Table 2)

表2 数据产出质量情况一览表 Table 2 List of sequencing data quality 样本 Sample 原始序列 Raw reads 过滤序列 Clean reads 过滤碱基 Clean bases 碱基错误率 Error rate (%) Q20 (%) Q30 (%) GC含量 GC content (%) FB_H 59 360 060 58 121 920 8.72G 0.02 95.74 89.72 43.88 MB_H 74 623 240 73 150 088 10.97G 0.02 96.33 91.14 43.74 ZD_H 69 569 630 68 216 824 10.23G 0.02 96.27 90.97 43.32 FB_H is male parent flowers; MB_H is female parent flowers; ZD_H is sterility flowers. FB_H为父本花; MB_H母本花; ZD_H为不育花。

表2 数据产出质量情况一览表 Table 2 List of sequencing data quality

图1

Fig. 1

	Figure Option ViewDownloadNew Window
	图1 Unigene与转录本长度分布Fig. 1 Length distribution of Unigene and transcripts

对于没有参照基因的茶树来说, 根据Trinity拼接得到的转组本作为参考序列, 将每个样本的clean reads往ref上绘制mapping。父本花获得42 291 172条mapping, 占clean reads 72.76%。母本花获得53 537 284条mapping, 占clean reads 73.19%。子代不育花获得48 115 182条mapping, 占clean reads的70.53%。再利用RSEM的bowtie2对bowtie软件进行比对, 并对比对的结果进行FPKM转换及FPKM密度分布(图2)。从图中看出父本花、母本花、子代不育花之间的基因表达存在差异。

图2

Fig. 2

	Figure Option ViewDownloadNew Window
	图2 不同样的基因表达水平比对 FB_H为父本花; MB_H母本花; ZD_H为不育花。Fig. 2 Comparison of genes expression of different samples FB_H is male parent flowers; MB_H is female parent flowers; ZD_H is sterility flowers.

对于无生物学重复的样品, 先采用TMM对read count数据进行标准化处理, 再用DEGseq进行差异分析, 为了控制假阳性率, 需Q-value结合Fold_ change来筛选, 差异基因筛选条件为阈值 < 0.005 & |log₂Fold_change| > 1。统计各组比较得到的差异基因个数并作维恩图(图3)。根据研究目的, 不育相关基因选择在父本花与子代花、母本花与子代之间共有的且在父本与母本之间没有的, 拟定为不育相关基因, 共筛选出1219个。

图3

Fig. 3

	Figure Option ViewDownloadNew Window
	图3 差异基因维恩图 FB_H为父本花; MB_H母本花; ZD_H为不育花。Fig. 3 Venn image of differential genes FB_H is male parent flowers; MB_H is female parent flowers; ZD_H is sterility flowers.

对获得的不育相关基因进行GO功能显著性富集分析, 在1219个差异基因中有952个被注释, 占78.1%, 并被分为生物过程(biological process, BP)、分子功能(molecular function, MF)和细胞组件(cellular component, CC) 3个功能注释, 进一步划分为24个类别。其中生物过程的类别最多, 有16个, 分子功能有7个类别, 细胞组件仅有1细胞壁类别。因为同一基因可能含有1个或多个功能注释, 从图4中可以看出参与代谢过程的最多, 为625个, 其次是催化活性, 为569个, 推测茶树不育花中的代谢是紊乱的。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 不育基因的GO分类 BP: 生物过程; MF: 分子功能; CC: 细胞组件。Fig. 4 GO classification of sterility genes BP: biological process; MF: molecular function; CC: cellular component.

在KEGG数据库中对不育相关基因进行pathway注释和分类表明, 不育相关基因注释到116个pathway, 涉及到核酸、氨基酸、糖以及次生代谢等多个方面。对20条富集最显著的pathway绘制散点图(图5)。从图可以看出富集到碳代谢途径的基因条数最多, 其次是氨基酸的生物合成。从富集显著性看, 半胱氨酸和蛋氨酸代谢、戊糖和葡萄糖醛酸转换、酒精中毒、三羧酸循环等最为显著。

图5

Fig. 5

	Figure Option ViewDownloadNew Window
	图5 不育基因的KEGG pathway富集Fig. 5 KEGG pathway of sterility genes

生长素作为植物中一种重要的内源激素, 在植物花的生长发育中发挥重要作用。不育基因的KEGG代谢通路分析表明, 在植物激素信号转导途径(ko04075)中有16个与生长素转导有关(表3), 其中在不育花中下降的基因包括生长素载体AUX1 LAX家族、生长素响应的GH3家族、生长素反应蛋白IAA和SAUR家族, 上调的为生长素响应因子ARF。

表3

Table 3

表3(Table 3)

表3 生长素转导相关的差异表达基因 Table 3 Differentially expressed genes of auxin signal transduction pathway No. 编号 Code 基因名称 Gene name FB_HFPKM MB_HFPKM ZD_HFPKM 表达 Expression 1 c119992_g1 Auxin influx carrier (AUX1 LAX family) 23.91 24.95 10.01 Down 2 c119992_g2 Auxin influx carrier (AUX1 LAX family) 28.79 43.88 15.12 Down 3 c109936_g1 Auxin influx carrier (AUX1 LAX family) 37.49 16.49 3.03 Down 4 c71040_g1 Auxin influx carrier (AUX1 LAX family) 39.68 21.19 4.28 Down 5 c121918_g1 Auxin responsive GH3 gene family 38.08 15.96 12.69 Down 6 c115411_g2 Auxin responsive GH3 gene family 13.53 13.98 1.00 Down 7 c122907_g1 Auxin response factor 20.11 26.66 64.89 Up 8 c124184_g3 Auxin response factor 35.64 43.92 104.96 Up 9 c112874_g2 Auxin response factor 14.57 11.46 93.39 Up 10 c88952_g1 Auxin response factor 1.80 2.12 31.62 Up 11 c96398_g1 Auxin response factor 31.93 24.36 118.90 Up 12 c89059_g1 Auxin-responsive protein IAA 36.04 31.57 0.93 Down 13 c109320_g1 Auxin-responsive protein IAA 134.19 92.74 26.87 Down 14 c101070_g1 SAUR family protein 228.56 113.47 1.97 Down 15 c103712_g1 SAUR family protein 87.82 13.56 0.06 Down 16 c95331_g1 SAUR family protein 129.99 44.03 7.54 Down FB_H is male parent flowers; MB_H is female parent flowers; ZD_H is sterility flowers. FB_H为父本花; MB_H母本花; ZD_H为不育花。

表3 生长素转导相关的差异表达基因 Table 3 Differentially expressed genes of auxin signal transduction pathway

随机筛选出5个差异表达基因c112358_g1、c112358_g5、c115764_g1、c116710_g4和c123710_g4进行qRT-PCR验证, 分别提取父本花、母本花、子代不育花总RNA, 反转录后进行定量PCR验证。并根据所得数据以log₂(FPKM)进行标准化处理。从图6可以看出, 5个差异表达基因在父本花、母本花以及子代不育花中的相对表达量与表达谱分析的变化趋势一致, 说明基因表达谱的分析结果可靠。

图6

Fig. 6

	Figure Option ViewDownloadNew Window
	图6 Illumina HiSeq测序和qRT-PCR的表达量比对 FB_H为父本花; MB_H母本花; ZD_H为不育花。Fig. 6 Validation of Illumina HiSeq sequencing data by qRT-PCR FB_H is male parent flowers; MB_H is female parent flowers; ZD_H is sterility flowers.