全 文 ：基金项目：国家自然科学基金(No. 31572066)、转基因生物研究与商业化重大专项(No. 2016ZX08011-001)和教育部高等学
校博士点基金优先发展领域课题(No. 20130097130006)
收稿日期：2016-02-26 接受日期：2016-04-05
利用拟南芥种群间转录组和表皮毛密度差异信息筛选表皮毛相关功能
基因
薛思鸣 成丹 解洪杰 何宝叶 强胜*
南京农业大学杂草研究室，南京 210095
*通讯作者，wrl@njau.edu.cn
摘 要 如何利用基因组、转录组以及代谢组学等海量信息确定基因功能，并阐述其调控机制是生命科
学研究领域的核心内容之一。本研究利用种群间数量性状遗传变异与相应基因表达量差异之间的关联
性，尝试一种功能基因筛选的快速方法，首先利用已经发表的拟南芥(Arabidopsis thaliana)生物信息数据，
运用线性相关模型，根据拟南芥种群遗传变异导致的叶片表皮毛数量性状差异和对应种群基因表达量差
异的相关性，筛选与表皮毛性状功能相关的基因，从33 554个基因中共筛选到1 747个相关基因，其中负
相关基因195个，正相关基因1 552个。对筛选基因集合进行基因本体(Gene Ontology, GO)富集分析，结
果表明，筛选到的基因主要富集于植物防御、先天免疫、胁迫反应、茉莉酸刺激应答、机械损伤反应和芥子
油苷代谢等生物功能，特别筛选到和表皮毛发育及其防御代谢物分泌有关的基因，如诱导叶片表皮毛分
化起始的MYB结构域蛋白基因MYB23，对叶片表皮毛发育具有重要影响的缺陷表皮毛因子1 (distorted
trichomes 1, DIS1)，参与吲哚芥子油苷合成的细胞色素 P450家族蛋白基因(cytochrome P450 family 81
subfamily F polypeptide 2, CYP81F2)和MYB51。将本研究筛选基因功能富集结果与前人对表皮毛的功能
研究成果做对比，证明了本研究中所用方法的可行性。因此，在与数量性状有关目的基因筛选的研究中，
此方法可以作为参考，对于其他物种的数量性状研究也具有借鉴意义。
关键词 转录组，数量性状，相关分析，表皮毛，功能基因
中图分类号 Q332 文献标识码 A
Screening Trichome- related Functional Genes Using Difference of
Transcriptome and Trichome Density Among Arabidopsis thaliana
Populations
XUE Si-Ming CHENG Dan XIE Hong-Jie HE Bao-Ye QIANG Sheng*
Weed Research Laboratory, Nanjing Agricultural University, Nanjing 210095, China
* Corresponding author, wrl@njau.edu.cn
Abstract How to confirm the functions of the genes and their regulation mechanisms by the massive
information such as genome, transcriptome and the metebolomics is one of the core tasks in the field of life
science. In present study, to make use the biological information data published in databases of Arabidopsis
thaliana, the trichome related functional genes were screened based on the correlation between the difference
of leaf trichome density variations and expression level differences of relevant genes among 31 Arabidopsis
populations, which was established by linear correlation analysis model. All 1 747 genes were screened as their
Online system: http://www.jabiotech.org
农 业 生 物 技 术 学 报
Journal of Agricultural Biotechnology
2016, 24(8): 1138~1146
DOI: 10.3969/j.issn.1674-7968.2016.08.004
农作物许多形态特征及其产量、食用品质、饲
用品质和抗逆性等大多数性状都是一些数量性状
(喻树迅,袁有禄, 2002)，数量性状有关基因的筛选
一直以来都是作物遗传育种的研究热点。在性状
有关基因的研究中，全基因组关联(genome-wide as-
sociation, GWAS)是一种主要发掘有效等位基因的
分析方法。Aranzana等(2005)利用 95份全基因组
多态性数据，通过筛选与开花时间和病原菌耐药性
有关基因来测试GWAS分析在拟南芥(Arabidopsis
thaliana)中的可行性，结果确定了与待测性状有关
的主要基因；Tian等(2011)利用GWAS分析，确定了
玉米(Zea mays)叶片结构形成的重要遗传基础，发
现了一些关键基因，而且证明无叶枕基因的变种可
导致更多的直立叶片。但是GWAS存在问题，很多
常见单核苷酸多态性 (single nucleotide polymor-
phism, SNP)很难阐明大多数性状的变异，原因可能
是遗传统计分析方法中多重假设检验和种群结构
导致的假阳性、假阴性问题，此外，GWAS难以检测
到罕见变异以及基因与环境之间、基因与基因之间
互作的统计效能等(段忠取,朱军, 2015)。因此仍然
需探索更多高通量的性状功能基因筛选方法。
物种内表型差异不仅与DNA序列多态性有
关，而且还会受到基因表达水平的影响(Druka et
al., 2010)，SNP对作物性状的影响最终还是体现到
影响基因的表达水平上(涂雨辰等, 2013)；特别是在
植物生长发育阶段，表达量持续且稳定的基因，其
表达量的高低可能决定了其数量性状差异。Zhang
等(2015)对番茄(Lycopersicon esculentum)的研究结
果表明，生长素响应因子(auxin response factor 3,
SIARF3)在表皮毛中含有较高的表达水平，SIARF3
的下调可以降低表皮细胞密度，而且明显降低叶片
表皮毛密度。豚草(Ambrosia artemisiifolia)是一种一
年生杂草，原产北美，入侵地较之于原产地，具有增
长快速和繁殖成功率高的特点，Hodgins等(2013)利
用芯片技术筛选5个原产地和 6个入侵地豚草种群
在3种不同环境(对照、光胁迫、养分胁迫)及两个时
间点下的差异表达基因，试图寻找能够解释豚草入
侵分子机制的备选基因，涉及的180个候选基因中
部分基因涉及次生代谢产物、胁迫反应及解毒反
应。Xie等(2015)对采自不同纬度和海拔高度的入
gene expression level was related to leaf trichome density from totally 33 554 genes. The number of negatively
related genes were 195, and positively related genes were 1 552. The positively correlated genes were
performed for function enrichment analysis. The results showed that the cluster terms mainly concentrated in
the 4 aspects, including stress response, cellular component organization, biological regulation and metabolic
processes. There were the most hierarchical clustering terms in the stress responses respect, including immune
response and innate immune response, stress response, defense response, wound response, response to
biological stimuli, fungal defense reaction, bacterial response, chemical stimulation of the response, organic
material response, jasmonic acid response to stimuli, in response to carbohydrate stimulation and chitin
response cluster terms. There were metabolism of glucosinolates, indole derivatives synthesis and metabolism,
indole glucosinolates metabolism, protein phosphorylation and phosphate metabolism cluster terms about the
metabolic process. There were cell wall modification and cell wall depolymerization cluster terms about the
cellular component organization, and there were signal transduction, stress response regulation and defense
response regulation cluster terms about biological regulation. In particular, there were genes directly related to
the development of trichome and genes directly related to the secretion of the defense metabolites of trichome
in the screened genes. Such as MYB23, a transcription factor which could induce the initiation of leaf trichome,
DIS1 (distorted trichomes 1) which had an important influence on the development of trichome, and CYP81F2
(cytochrome P450 family 81 subfamily F polypeptide 2) and MYB51 which were both involved in the synthesis
of indole glucosinolates. And the expression level of the 4 genes mentioned above all had a significant positive
correlation with the density of the leaf trichome. Comparison of the function enrichment results of screened
genes in the present research with the previous functions studies on leaf trichome proved practicability of the
present method. Therefore, this method may provide a reference to screen target functional genes related to
quantitative traits not only in A. thaliana but also in other species.
Keywords Transcriptome, Quantitative trait, Correlation analyses, Trichome, Functional gene
利用拟南芥种群间转录组和表皮毛密度差异信息筛选表皮毛相关功能基因
Screening Trichome-related Functional Genes Using Difference of Transcriptome and Trichome Density AmongArabidopsis thaliana 1139
农业生物技术学报
Journal of Agricultural Biotechnology
侵植物紫茎泽兰(Ageratina adenophora)，对其进行
地理种群耐冷性研究发现，C-重复基序结合因子(C-
repeat binding factor, CBF)介导的冷反应通路中大部
分成员基因本底表达量和种群本身的耐冷性成正
比，这种耐冷性分化与其在中国西南部从南到北的
入侵广度密切相关。尽管基因表达量的差异与某一
性状差异不一定有必然联系，但是可以依据差异数
量性状和差异基因表达量的相关性，利用统计学原
理先找到两者间的相关性，进一步有针对性地研究
选定基因表达量引起表型差异的规律，从而可以提
供挖掘新基因的方向和缩小探索功能基因的范
围。基于这个想法以及公共数据库中越来越丰富
的转录组测序数据(祁云霞等, 2011)，本研究选择表
皮毛密度差异性状，分析其与基因表达量的相关
性，尝试一种筛选表皮毛相关功能基因的方法，以
期为利用生物信息挖掘功能基因提供更多手段。
1 材料与方法
1.1 实验材料
Schmitz等(2013)对来自于北半球 144个种群
的野生拟南芥(Arabidopsis thaliana)进行 RNA-seq
测序，本研究以该研究项目中的转录组数据(NCBI
GEO 登录号: GSE43858, http://www.ncbi.nlm.nih.
gov/geo/query/acc.cgi?acc=GSE43858) 为 数 据 来
源。拟南芥对应种群的数量性状数据来自Atwell
等(2010)所做的全基因组关联分析中所使用的107
个差异表型性状指标 (https://github.com/Gregor-
Mendel- nstitute/atpolydb/blob/master/miscellane ous_
data/phenotype_published_raw.tsv)。
1.2 数据获取及整理
由于数据来源于不同实验，所涉及到的拟南芥
种群并不完全相同。鉴于此，本研究挑选了 31个
具有RNA-seq测序数据且具有表皮毛密度性状数
据的共有种群，并直接从NCBI GEO数据库下载31
个拟南芥种群对应的、在原实验中已经经过生物信
息学处理后得到的、以“每1百万个匹配读段中，匹
配到外显子上的每 1千个碱基的读段个数(frag-
ments per kilo bases per million reads, FPKM)”存贮
的基因表达量数据文件；共包含33 554个基因对应
的估计转录表达量 FPKM值(注：下文中所提到的
表达值均指FPKM)。
从1.1中提到性状数据网址可查看对应性状指
标数据，共包含107个拟南芥性状数据。本研究采
用了其中的表皮毛密度：
表皮毛密度=表皮毛个数/叶圆片面积(cm²)
数据详情可查阅Atwell等(2010)和 1.1中给出
的网址，表1中列出了本研究的31个拟南芥的种群
号及对应的表皮毛密度值。
1.3 数据分析
在不同拟南芥种群中同一基因的FPKM是有
差异的，本研究根据不同种群、同一基因不同表达
值和数量性状的成对数据组成的两组数据，进行线
性相关分析，以期寻找基因FPKM和性状指标之间
的相关性，找到相应基因集合。
为了分析每一个基因在多个种群中的 FPKM
与对应种群的待研究数量性状的数值之间是否具
有相关性，本研究利用R数据统计软件(R version
3.2.1，Platform: x86_64-w64-mingw32/x64 (64-bit))
编写了一段简单代码，能够提取本研究的表达量数
据集和性状数据集，以线性相关模型为依据，针对
每一个基因，对表达量和性状数据之间进行线性相
关性分析(R数据统计软件中的 cor.test函数)，得到
P值和皮尔森相关系数R值。然后设定P＜0.05为
显著相关，将具有显著相关的基因筛选出来。
为了进一步探究所筛选到的基因生物学功能，
并验证此方法的可行性，本研究结合DAVID (Data-
表1 31个种群对应表皮毛密度
Table 1 Trichome density of 31 populations
种群号
Population ID
An_1
Bor_4
Col_0
Ei_2
Ga_0
Got_7
Gu_0
Gy_0
HR_5
Kin_0
Kondara
Kz_9
Ler_1
Lp2_2
Mz_0
Nok_3
表皮毛密度
Trichome density
36.3
24.5
32.5
14
18.3
9.8
42.5
43.8
56.8
14
34.3
72
14.3
14.8
77.8
13
种群号
Population ID
Pna_17
Pu2_23
Pu2_7
Ra_0
Rmx_A02
Se_0
Sq_8
Tamm_2
Uod_1
Van_0
Wa_1
Ws_2
Wt_5
Yo_0
Zdr_1
表皮毛密度
Trichome density
53.8
10.8
37
17
24.8
30.5
30.7
48
15.5
20.8
17.5
33.3
17.3
25
10
1140
base for Annotation, Visualization and Integrated Dis-
covery)，及 Cytoscape 软件的 BINGO (Biological
Networks Gene Ontology Tool)插件进行基因本体
(Gene Ontology, GO)(Thomas et al., 2007)及生物学
过程富集分析(杨蓉,蔡琳, 2012)。
2 结果与分析
2.1 基因表达量和拟南芥表皮毛密度有相关性的
基因筛选和分析
本研究用表皮毛密度作为待研究的数量性状，
利用 1.3中所叙述的方法，共筛选到 1 747个基因，
其中正相关基因1 552个，占筛选到基因的89%，呈
负相关的基因有195个。
2.2 对筛选到的相关基因进行功能聚类分析
用两种方法Cytoscape软件的BINGO插件和
DAVID数据库的“Functional Annotation Clustering”
对这些基因进行基因功能富集分析。利用BINGO
对筛选到的正相关基因进行生物过程方面的富集
分析，形成可视化图(图 1)，图 1中有颜色的圆圈代
表在P＜0.01水平下极显著的富集结果，颜色越深，
显著性越强。链接拟南芥GO数据库得到筛选的
正相关基因集合中1 023个基因的生物过程注释信
息，靶基因显著富集在应激反应、细胞组分组织、生
物调节、代谢过程 4大方面；其中应激反应包含的
各层次聚类词条最多，包括免疫反应、天然免疫应
答、胁迫反应、防御反应、伤口反应、响应生物刺激、
真菌防御反应、细菌应答、化学刺激应答、有机物质
应答、茉莉酸刺激应答、响应碳水化合物刺激、几丁
质应答；在代谢过程上有硫代葡糖苷代谢、吲哚衍
生物合成与代谢、吲哚族芥子油苷代谢、蛋白质磷
酸化修饰、磷酸盐代谢；在细胞组分组织上有细胞
壁解聚和修饰；在生物调节上有信号转导、胁迫反
应调节、防御反应调节(图1)。
表 2是BINGO软件对所有相关基因在显著水
平P＜0.05标准下的生物过程方面全部富集结果，
据此富集结果可以归纳出这批基因显著参与的生
物过程，包括免疫反应、胁迫反应、茉莉酸刺激反
应、机械损伤刺激、信号传导、细胞壁加工、芥子油
苷代谢、蛋白磷酸化和细胞过程调节；其中参与胁
迫应答的基因有 152个，参与响应刺激的基因有
226个，参与化学刺激反应的基因有130个，参与蛋
白磷酸化的基因有 77个，调节细胞过程的基因有
160个，大分子代谢过程的基因有103个。
表3是DAVID软件对正相关基因富集结果，分
别是甲壳素反应、糖刺激反应、免疫反应、防御反
应、先天免疫应答、蛋白氨基酸磷酸化、磷代谢、细
胞死亡、程序性细胞死亡、细胞凋亡。与BINGO的
分析结果基本一致，显著性富集交集有：甲壳素反
应、糖刺激反应、免疫反应、防御反应、先天免疫应
答、磷酸化、刺激反应等；略有不同的是BINGO分
析结果中还有胁迫反应、响应其他生物过程的类
别，DAVID分析结果中还有细胞死亡等类别。
2.3 筛选基因集合中包含与表皮毛发育及代谢直接
相关的基因
另外，本研究筛选到的基因集合还包含已知的
已鉴定出的和表皮毛发育及防御代谢物有关的基
因，如诱导叶片表皮毛分化起始的MYB结构域蛋
白基因MYB23(Kirik et al., 2001)和缺陷表皮毛因子
1 (distorted trichomes 1, DIS1)，DIS1对表皮毛发育
具有重要影响，其突变体会造成叶片表皮毛生长缺
陷(Jie et al., 2003)；与防御相关的次级代谢物合成
及调节有关的基因，如编码细胞色素 P450家族蛋
白基因 (cytochrome P450 family 81 subfamily F
polypeptide 2, CYP81F2)，其控制吲哚族芥子油苷
合成，有助于增强对桃蚜 (Myzus persicae)的抗性
(Marina et al., 2009)；MYB51基因属于 R2R3-MYB
转录因子家族，参与芥子油苷合成(Tamara et al.,
2007)。
为了更形象直观地揭示表达量和叶片表皮毛
密度的关系，将不同种群的这4个基因FPKM值和
相应种群叶片表皮毛密度做散点图线性回归分析，
结果均呈现出显著的正相关关系(图 2)；图中的点
大致分布在左下角和右上角区域，且除表皮毛密度
最大值对应的两个点距离趋势线较远外，其他点较
为均匀的分布在趋势线两侧，并表现出表皮毛密度
随基因表达量增加而增加的趋势，统计结果表明，
表皮毛密度与基因表达量呈现显著正相关性。
3 讨论
表皮毛作为植物与环境交互作用的直接接触
界面，在植物与生物和非生物相互作用中所起的积
极作用被很多生态学实验证实。比如，表皮毛数量
的增加可以增强对草食昆虫的抵抗性(Mauricio,
利用拟南芥种群间转录组和表皮毛密度差异信息筛选表皮毛相关功能基因
Screening Trichome-related Functional Genes Using Difference of Transcriptome and Trichome Density AmongArabidopsis thaliana 1141
农业生物技术学报
Journal of Agricultural Biotechnology
Rausher, 1997; Clauss et al., 2006)。 Reymond 等
(2004)发现拟南芥专食性昆虫菜粉蝶(Pieris rapae)
在无表皮毛突变体gl1中较野生型(col-0)啃食表现
更好。而 Jakoby等(2008)的GWAS表达分析表明，
花青素、黄酮类和芥子油苷化合物合成通路中的一
些基因在表皮毛中的表达量较叶片中具有相对更
高的水平。这些和防御有关的次级代谢的化合物
存在，表明表皮毛在植物防御和保护中具有重要作
用。Traw和Bergelson(2003)首次实验证实人工损
伤和茉莉酸处理可以增加拟南芥叶片表皮毛密度
和数量。在本研究中，筛选到的目的基因的生物过
程显著富集词条与前人对表皮毛的功能研究结果
十分相符，强有力的证实了本研究方法的有效性。
传统的基因挖掘方法有突变体法、图位克隆法
及比较基因组学法，目前高通量的方法主要是全基
因组关联分析法和基因表达差异法。其中GWAS
分析是根据连锁不平衡(linkage disequilibrium, LD)
原理，将群体的目标性状与遗传标记联系起来寻找
与目标性状有关基因的方法(Myles et al., 2009)。
此方法在人类遗传疾病及动植物复杂性状研究中
应用广泛，也暴露出很多缺点，由于GWAS识别关
联位点的能力与等位基因频率有关，其对于在表型
变异上虽有较大贡献但基因频率较低的变异很难
识别(Pritchard et al., 2000)；复杂性状的表型变异除
图1 表皮毛差异筛选正相关基因GO生物过程分类的层次网络
Figure 1 The hierarchical network of the screening for the positive correlation genes by GO analysis in the process of bi⁃
ological process
圆圈大小代表基因数量多少，颜色深浅代表P值大小
Circle size represents the number of genes, and color depth represents the P value
1142
了受位点的遗传变异影响，还受其他非碱基变异分
子机制影响，如拷贝数变化(copy number variation,
CNV)(Freeman et al., 2006)和基因甲基化的影响，
Xiangqian等(2015)研究发现，RAV6/ABI3 (related to
abscisic acid insensitive 3)基因的甲基化变异导致
的基因表达差异会影响水稻(Oryza sativa)叶角和种
表2 BINGO软件对筛选到的表皮毛密度相关基因的生物过程富集
Table 2 Biological process enrichment by BINGO for screened genes which are correlated with trichome density
GO号
GO ID
10200
9743
6952
6955
2376
6950
45087
10033
51707
9607
50896
51704
6468
42221
16310
9611
43687
42430
80134
6796
9753
31347
9617
44277
9620
6464
42343
23060
50794
23052
31349
6541
48583
2218
52543
10556
19757
9723
50778
9816
43412
矫正P值
Corrected P value
9.01E-21
4.23E-15
5.00E-08
1.55E-07
1.55E-07
7.99E-07
7.99E-07
4.80E-06
6.30E-06
1.00E-05
3.02E-05
5.88E-05
1.30E-04
3.47E-04
1.02E-03
1.19E-03
1.47E-03
2.12E-03
2.44E-03
2.89E-03
3.97E-03
5.39E-03
6.00E-03
6.00E-03
6.58E-03
6.87E-03
8.60E-03
2.23E-02
2.78E-02
2.78E-02
2.78E-02
3.53E-02
3.71E-02
4.01E-02
4.81E-02
4.81E-02
4.87E-02
4.88E-02
4.92E-02
4.92E-02
4.92E-02
x
39
43
73
41
42
152
38
95
58
59
226
67
76
130
77
20
88
10
14
79
20
12
27
3
20
93
4
46
160
76
7
4
18
6
5
103
8
15
6
6
95
功能描述
Description
甲壳素应答Response to chitin
响应有机物刺激Response to carbohydrate stimulus
防御反应Defense response
免疫反应 Immune response
免疫系统过程 Immune system process
胁迫应答Response to stress
天然免疫反应 Innate immune response
有机物质应答Response to organic substance
对其他生物的反应Response to other organism
响应生物刺激Response to biotic stimulus
响应刺激Response to stimulus
多生物过程Multi-organism process
蛋白氨基酸磷酸化Protein amino acid phosphorylation
化学刺激反应Response to chemical stimulus
磷酸化Phosphorylation
机械损伤反应Response to wounding
蛋白质翻译后修饰Post-translational protein modification
吲哚及其衍生物的代谢过程 Indole and derivative metabolic process
胁迫反应调节Regulation of response to stress
磷酸盐代谢过程Phosphate metabolic process
响应水杨酸刺激Response to jasmonic acid stimulus
调节防御反应Regulation of defense response
细菌应答Response to bacterium
细胞壁解聚Cell wall disassembly
真菌反应Response to fungus
蛋白修饰过程Protein modification process
吲哚芥子油苷的代谢过程 Indole glucosinolate metabolic process
信号传导Signal transmission
调节细胞过程Regulation of cellular process
信号传递Signaling
防御反应正调节Positive regulation of defense response
谷氨酰胺代谢Glutamine metabolic process
调节刺激反应Regulation of response to stimulus
激活天然免疫Activation of innate immune response
细胞壁胼胝质沉积Callose deposition in cell wall
大分子生物合成过程的调控Regulation of macromolecule biosynthetic process
芥子油苷代谢Glycosinolate metabolic process
响应乙烯刺激Response to ethylene stimulus
对免疫反应的正调控Positive regulation of immune response
细菌的防御反应，不相容的相互作用
Defense response to bacterium, incompatible interaction
大分子修饰Macromolecule modification
n
104
177
637
272
283
1853
257
1037
528
550
3207
694
842
1710
910
133
1093
41
78
979
148
65
241
3
156
1241
7
541
2448
1023
31
10
157
25
18
1504
45
125
27
27
1378
P值的矫正采用Benjamini法；x：待分析基因在当前功能类的数量；n：当前功能类基因总数。筛选基因中总共有1148个基因
具有对应注释信息
Benjamini methods was used for the P value correction; x: The total number of stay analysis of gene in the function class number;
n: The total number of functional classes of genes. A total of 1148 genes have corresponding annotation information
利用拟南芥种群间转录组和表皮毛密度差异信息筛选表皮毛相关功能基因
Screening Trichome-related Functional Genes Using Difference of Transcriptome and Trichome Density AmongArabidopsis thaliana 1143
农业生物技术学报
Journal of Agricultural Biotechnology
子大小，对这种非基因位点变异导致的性状差异，
GWAS就显得无能为力；由于遗传位点还可以与环
境之间产生互作，因此性状的表型差异会受到环境
影响，不考虑此因素，会导致GWAS研究中一部分
未能解释的遗传率(Thomas, 2010)。
基因表达差异分析技术不考虑基因位点变异
和涉及庞大的群体，比较对照组和实验组两种状态
下个体的基因表达差异，就可以确定目标基因功能
(Costa et al., 2010)。表达差异发现功能基因更多
是基于对照组与处理组之间的表达差异，来发现与
环境应激有关的基因。物种内表型差异不仅与
DNA序列多态性有关，而且还会受到基因表达水
平的影响(Druka et al., 2010)，SNP对作物性状的影
响最终还是体现到影响基因的表达水平上(涂雨辰
等,2013)。本研究在表达差异技术的应用上进行了
拓展，将多种群数量性状与多种群基因表达量联系
起来进行相关性研究，不仅能够利用前人GWAS研
究积累的差异表型数据资料，还可以利用越来越丰
富的多种群转录组数据，具有高通量、低成本，不受
限于已知基因组信息及庞大的SNP位点信息，且方
法只涉及简单的线性相关，能够避免GWAS分析中
低基因频率导致的低关联度，非位点变异导致的表
达量差异，基因与环境之间的互作等问题。缺点是
基因表达具有组织、时间、条件特异性，所以在选材
上要比较苛刻。当然随着测序及基因芯片技术的
发展，相信多种群、多条件、多组织的表达数据会越
来越丰富，从而在与植物数量性状有关功能基因的
筛选中此方法可能会得到应用。
表3 DAVID对筛选到的表皮毛密度正相关基因的生物过程富集
Table 3 Biological process enrichment by DAVID for screened genes which are positively correlated with trichome density
GO号
GO ID
0010200
0009743
0010033
0006955
0006952
0045087
0006468
0016310
0006796
0006793
0016265
0008219
0012501
0006915
P值
P value
3.30E-24
1.75E-18
4.01E-08
2.30E-11
1.64E-10
5.13E-10
1.18E-05
1.47E-04
1.87E-04
1.97E-04
0.00661061
0.00661061
0.02513425
0.04203342
矫正后P值
Corrected P value
3.44E-21
9.10E-16
6.95E-06
7.98E-09
4.28E-08
1.07E-07
0.00174634
0.015147814
0.014866191
0.014493629
0.279975129
0.279975129
0.551673844
0.691262937
基因数
Gene number
43
46
100
45
98
41
84
86
91
91
23
23
19
15
基因编号
Gene ID
AT5G59820, AT3G50060,
AT3G28210, AT4G17500 et al.
AT1G62300, AT4G33050,
AT5G67450, AT4G18880 et al.
AT3G28210, AT1G74020,
AT4G17500, AT5G27380 et al.
AT4G34460, AT3G04220,
AT1G60320, AT5G44070 et al.
AT1G64107, AT3G04220,
AT4G17500, AT3G52400 et al.
AT4G34460, AT3G04220,
AT1G60320, AT5G37770 et al.
AT3G25250, AT1G30640,
AT4G34440, AT1G68690 et al.
AT1G23540, AT1G49160,
AT5G48940, AT1G61610 et al.
AT3G59740, AT2G18530,
AT2G02800, AT5G61560 et al.
AT2G30040, AT5G60310,
AT3G09010, AT4G23280 et al.
AT3G04210, AT3G04220,
AT1G17610, AT5G44070 et al.
AT1G17610, AT5G44070,
AT1G12210, AT1G58410 et al.
AT1G64070, AT5G38350,
AT1G65850, AT1G15890 et al.
AT1G58410, AT1G59620,
AT4G10780, AT3G25510 et al.
生物过程注释
Biological process annotation
甲壳素响应
Response to chitin
糖刺激的反应
Response to carbohydrate stimulus
有机物质的反应
Response to organic substance
免疫反应
Immune response
防御反应
Defense response
先天免疫应答
Innate immune response
蛋白质氨基酸的磷酸化
Protein amino acid phosphorylation
磷酸化
Phosphorylation
磷酸盐的代谢过程
Phosphate metabolic process
磷的代谢过程
Phosphorus metabolic process
死亡
Death
细胞死亡
Cell death
细胞程序性死亡
Programmed cell death
细胞凋亡
Apoptosis
1544个具有注释信息，根据P值大小只给出富集分析结果中的前14个生物过程类别
All 1544 of the screened genes were annotated, and the first 14 biological processes in the enrichment analysis were listed based
on the P value
1144
4 结论
本研究基于拟南芥种群间叶片表皮毛密度和
基因表达量的相关性，从 33 554个基因中共筛选
到 1 747个相关基因，其中负相关基因 195个，正
相关基因 1 552个。对筛选基因集合进行GO富
集分析，结果表明，筛选到的基因主要富集于植
物防御、先天免疫、胁迫反应、茉莉酸刺激应答、机
械损伤反应和芥子油苷代谢等生物功能，特别筛
选到和表皮毛发育及其防御代谢物分泌有关的基
因。筛选基因的功能富集结果符合前人对表皮毛
的功能研究成果，证实了方法的可行性，在与数量
性状有关目的基因筛选的研究中，此方法具有借
鉴意义。
参考文献
段忠取,朱军 . 2015.全基因组关联分析研究进展[J].浙江大
学学报, 41(4): 385-393. (Duan Z H, Zhu J. 2015. Re-
search progress of genome- wide association study[J].
Journal of Zhejiang University, 41(4): 385-393.)
祁云霞, 刘永斌, 荣威桓 . 2011. 转录组研究新技术: RNA-
Seq及其应用 [J]. 遗传, 33(11): 1191- 1202. (Qi Y X,
Liu Y B, Rong W H. 2011. RNA- Seq and its applica-
tions: A new technology for transcriptomics[J]. Heredi-
tas, 33(11): 1191-1202.)
涂雨辰,田云,卢向阳 . 2013.全基因组关联分析在植物中的
应用[J]. 化学与生物工程, 30(6): 7-10. (Tu Y C, Tian
Y, Lu X Y. 2013. The application of whole genome asso-
ciation analysis in plants[J]. Chemistry & Bioengineer-
ing, 30(6): 7-10.)
图2 直接参与表皮毛发育基因的相对表达量和相应种群叶片表皮毛密度的散点图
Figure 2 The scatter plot by the expression value of the genes and the leaf trichome density
A~D：MYB23、缺陷表皮毛因子1(DIS1)、细胞色素P450家族蛋白基因(CYP81F2)和MYB51基因的FPKM值与叶片表皮毛密
度相互关系；FPKM：每1百万个匹配读段中，匹配到外显子上的每1千个碱基的读段个数
A~D: Relative expression (FPKM values) and leaf surface trichome density relationship of MYB23, distorted trichomes 1 (DIS1), cy-
tochrome P450 family 81 subfamily F polypeptide 2 (CYP81F2) and MYB51; FPKM: Fragments per kilo bases per million reads
表
皮
毛
密
度
/(
个·
cm
- ²)
L
ea
f
tr
ic
ho
m
e
de
ns
it
y
80
70
60
50
40
30
20
10
5 15 25 ×105
MYB51的FPKM
FPKM of MYB51 D
A
表
皮
毛
密
度
/(
个·
cm
- ²)
L
ea
f
tr
ic
ho
m
e
de
ns
it
y
80
70
60
50
40
30
20
10
0 2 4 6 ×105
MYB23的FPKM
FPKM of MYB23
表
皮
毛
密
度
/(
个·
cm
- ²)
L
ea
f
tr
ic
ho
m
e
de
ns
it
y
80
70
60
50
40
30
20
10
15 25 35 ×104
DIS1的FPKM
FPKM of DIS1 B
表
皮
毛
密
度
/(
个·
cm
- ²)
L
ea
f
tr
ic
ho
m
e
de
ns
it
y
80
70
60
50
40
30
20
10
2 6 10 14 ×105
CYP91F2的FPKM
FPKM of CYP91F2 C
利用拟南芥种群间转录组和表皮毛密度差异信息筛选表皮毛相关功能基因
Screening Trichome-related Functional Genes Using Difference of Transcriptome and Trichome Density AmongArabidopsis thaliana 1145
农业生物技术学报
Journal of Agricultural Biotechnology
杨蓉,蔡琳, 2012. BINGO及DAVID在miR-155靶基因富集
分析中的作用[J].福建医科大学学报, 46(6): 408-414.
(Yang R, Cai L. 2012. The application of BINGO and
DAVID in the enrichment analysis of miR- 155 target
genes[J]. Academic Journal of Fujian Medical Universi-
ty, 46(6): 408-414.)
喻树迅,袁有禄 . 2002.数量性状遗传研究的新进展[J].棉花
学报, 14(3): 180-184. (Yu S X, Yuan Y L. 2002. New
progress of genetic study on quantitative traits[J]. Cot-
ton Science, 14(3): 180-184.)
Aranzana M J, Kim S, Zhao K, et al. 2005. Genome-wide as-
sociation mapping in Arabidopsis identifies previously
known flowering time and pathogen resistance genes[J].
PLOS Genet, 1(5): e60.
Atwell S, Huang Y S, Vilhjálmsson B J, et al. 2010. Genome
wide association study of 107 phenotypes in Arabidopsis
thaliana inbred lines[J]. Nature . 465(7298): 627-631.
Clauss M J, Dietel S, Schubert G, et al. 2006. Glucosinolate
and trichome defenses in a natural Arabidopsis lyrata
population[J]. Journal of Chemical Ecology, 32(11):
2351-2373.
Costa V, Angelini C, De F I, et al. 2010. Uncovering the com-
plexity of transcriptomes with RNA-Seq[J]. Biomed Re-
search International, 2010(1): 853916.
Druka A, Potokina E, Luo Z, et al. 2010. Expression quantita-
tive trait loci analysis in plants[J]. Plant Biotechnology
Journal, 8(1): 10-27.
Freeman J L, Perry G H, Feuk L, et al. 2006. Copy number
variation: New insights in genome diversity.[J]. Ge-
nome Research, 16(8): 949-961.
Hodgins K A, Zhao L, Nurkowski K, et al. 2013. The molecu-
lar basis of invasiveness: differences in gene expression
of native and introduced common ragweed (Ambrosia ar⁃
temisiifolia) in stressful and benign environments[J].
Molecular Ecology, 22(9): 2496-2510.
Jakoby M J, Doris F, Mader M T, et al. 2008. Transcriptional
profiling of mature Arabidopsis trichomes reveals that
NOECK encodes the MIXTA-like transcriptional regula-
tor MYB106[J]. Plant Physiology, 148(3): 1583-1602.
Jie L, El- Assal E D, Basu D, et al. 2003. Requirements for
Arabidopsis ATARP2 and ATARP3 during epidermal de-
velopment.[J]. Current Biology Cb, 13(15):1341-1347.
Kirik V, Schnittger A, Radchuk V, et al. 2001. Ectopic expres-
sion of the Arabidopsis AtMYB23 gene induces differen-
tiation of trichome cells[J]. Developmental Biology, 235
(2): 366-377.
Marina P, Vogel H, Kroymann J. 2009. The gene controlling
the indole glucosinolate modifier1 quantitative trait lo-
cus alters indole glucosinolate structures and aphid resis-
tance in Arabidopsis[J]. Plant Cell, 21(3): 985-999.
Mauricio R, Rausher M D. 1997. Experimental manipulation
of putative selective agents provides evidence for the
role of natural enemies in the evolution of plant disease
[J]. Evolution, 51(5): 1435.
Myles S, Peiffer J, Brown P J, et al. 2009. Association map-
ping: Critical considerations shift from genotyping to ex-
perimental design.[J]. Plant Cell, 21(8):2194-2202.
Pritchard J K, Stephens M, Rosenberg N A, et al. 2000. Asso-
ciation mapping in structured populations[J]. American
Journal of Human Genetics, 67(1):170-181.
Reymond P, Rm B N P, Krishnamurthy V, et al. 2004. A con-
served transcript pattern in response to a specialist and a
generalist herbivore[J]. Plant Cell, 16(11): 3132-3147.
Schmitz R J, Schultz M D, Urich M A, et al. 2013. Patterns of
population epigenomic diversity[J]. Nature, 495(7440):
193-198.
Tamara G, Berger B, Mock H P, et al. 2007. The transcription
factor HIG1/MYB51 regulates indolic glucosinolate bio-
synthesis in Arabidopsis thaliana[J]. Plant Journal, 50(5):
886-901.
Thomas D. 2010. Gene- environment- wide association stud-
ies: Emerging approaches[J]. Nature Reviews Genetics,
11(4): 259-272.
Thomas P D, Mi H, Lewis S. 2007. Ontology annotation:
Mapping genomic regions to biological function[J]. Cur-
rent Opinion in Chemical Biology, 11(1): 4-11.
Tian F, Bradbury P J, Brown P J, et al. 2011. Genome-wide as-
sociation study of leaf architecture in the maize nested
association mapping population[J]. Nature Reviews Ge-
netics, 43(2): 159-162.
Traw M B, Bergelson J. 2003. Interactive effects of jasmonic
acid, salicylic acid, and gibberellin on induction of tri-
chomes in Arabidopsis[J]. Plant Physiology, 133(3):1367-
1375.
Xiangqian Z, Jing S, Xiaofeng C, et al. 2015. Epigenetic muta-
tion of RAV6 affects leaf angle and seed size in rice[J].
Plant Physiology, 169(3): 2118-2128.
Xie H J, Li H, Liu D, et al. 2015. ICE1 demethylation drives
the range expansion of a plant invader through cold toler-
ance divergence[J]. Molecular Ecology, 24(4): 835-850.
Zhang X, Yan F, Tang Y, et al.2015. Auxin response gene
SlARF3 plays multiple roles in tomato development and
is involved in the formation of epidermal cells and tri-
chomes.[J]. Plant & Cell Physiology, 56(11).
(责任编辑 靳晓霞)
1146