海岛棉GH9基因家族成员鉴定及分析

刘正文, 王省芬, 孟成生, 张艳, 孙正文, 吴立强, 马峙英, 张桂寅

(河北农业大学农学院, 华北作物改良与调控国家重点实验室, 河北省作物种质资源重点实验室, 河北 保定 071000)

棉花是重要的经济作物，为纺织工业提供天然纤维原料。陆地棉(Gossypiumhirsutum)和海岛棉(Gossypiumbarbadense)是目前栽培范围比较广泛的两个棉种，陆地棉产量高但纤维品质中等，海岛棉纤维品质优异但产量低。棉纤维品质主要由细胞壁特性决定，因此挖掘分析海岛棉纤维细胞壁发育相关基因有利于明确优质棉纤维形成的分子机理，同时为分子育种提供候选基因，有利于实现高产和优质性状的有效聚合。

植物GH9 (glycoside hydrolase 9) 家族由内切-β-1，4-葡聚糖酶(endo-β-1, 4-glucanase)构成，能够水解具有β-1，4-葡聚糖主链的多聚糖，并且涉及纤维素的生物合成，在细胞壁的生物合成和重塑中起着重要作用。区别于细菌，植物GH9基因的潜在底物是无定型区的纤维素以及非结晶多聚糖(如木葡聚糖)。依据结构不同，可将植物GH9基因分为A、B、C三个亚族[1]。A亚族N端有一个跨膜区域，是膜结合蛋白。在拟南芥中，GH9A (KOR) 是纤维素合酶复合体的重要成分[2-3]。在杨树[4-5]和陆地棉[6]中沉默KOR显著减少结晶纤维素含量。可见，GH9A参与植物细胞壁纤维素生物合成。B亚族通常包含信号肽，是分泌蛋白，在植物GH9基因家族中数量最多，具有复杂多变的功能。沉默拟南芥GH9B1(CEL1) 造成细胞壁褶皱，同时伴随减少的纤维素和木质素含量[7]。在水稻中，OsGH9B8、OsGH9B9、OsGH9B10、OsGH9B11与OsCESA4、OsCESA7、OsCESA9共表达，表明其在次生壁形成过程中起作用[8]。此外，在水稻中过表达OsGH9B1或OsGH9B3只会轻微改变纤维素、半纤维素、木质素组分，但会显著减少纤维素聚合度和结晶指数，说明OsGH9B1/B3能够修饰纤维素微纤丝[9]。C亚族结构类似B亚族，具有信号肽，此外其C端包含植物特有的CBM49 (carbohydrate binding module 49) 结构。细菌纤维素酶通常具备CBM结构，用于附着结晶纤维素。体外试验发现，来自番茄SlCel9C1的CBM结构能结合结晶纤维素[10]。然而，也有证据显示CBM结构会在翻译后被切除[11]。目前，尚无植物GH9C在体内水解结晶纤维素的强力证据。

随着植物基因组信息的完善，研究者已从拟南芥[1]、水稻[8]、杨树[12]、陆地棉[13]中鉴定出GH9基因家族成员。相较于陆地棉，海岛棉拥有更优质的棉纤维，在海岛棉中鉴定并分析GH9基因能够加深对棉纤维品质性状形成的理解，补充分子育种资源，然而目前尚无海岛棉GH9基因家族相关报道。本研究基于海岛棉品种Hai7124基因组鉴定GH9基因家族成员，并分析其理化性质、基因结构、染色体分布、进化历程、表达模式、转录调控和潜在功能，为棉花GH9基因的深入研究奠定基础。

1 材料与方法

1.1 基因组序列

海岛棉(GossypiumbarbadenseZJU_Hai7124_V1.1)、陆地棉(GossypiumhirsutumZJU_TM-1_V2.1)序列下载自CottonFGD (https://cottonfgd.org/) 。大豆(GlycinemaxWm82.a2.v1)序列下载自Phytozome (https://data.jgi.doe.gov/refine-download/phytozome) 。葡萄(Vitisvinifera12X)、黄麻(CorchoruscapsularisCCACVL1_1.0)、可可(TheobromacacaoCriollo_cocoa_genome_V2)序列下载自Ensembl Plants (http://plants.ensembl.org/index.html) 。榴莲(DuriozibethinusGCF_002303985.1_Duzib1.0)序列下载自NCBI FTP (ftp://ftp.ncbi.nih.gov/genomes/refseq/plant/) 。

1.2 GH9基因家族成员鉴定及理化性质分析

从Pfam (http://pfam.xfam.org/) 数据库下载GH9家族(PF00759)结构域序列文件，利用HMMER 3.0[14]软件建立隐马尔可夫模型并鉴定候选基因，利用SMART (http://smart.embl-heidelberg.de) 和CDD (http://www.ncbi.nlm.nih.gov/cdd) 分析候选蛋白序列，剔除不含完整保守结构域的序列。利用SignalP-5.0 (http://www.cbs.dtu.dk/services/SignalP/)预测信号肽，TMHMM-2.0 (http://www.cbs.dtu.dk/services/TMHMM/) 预测跨膜域。蛋白长度、分子量、等电点利用本地Perl脚本批量计算。

1.3 保守基序、基因结构、染色体定位及系统进化分析

使用MEME (http://meme-suite.org/)工具分析家族保守基序，motif数量设置为15。基因结构信息和染色体位置信息提取自基因组结构注释文件，并分别用Tbtools[15]和MapChart[16]软件展示。MEGA 7.0[17]软件构建系统发育树，基于邻接法(neighbor-joining method)，Bootstrap设置为1 000次。在线工具iTOL (https://itol.embl.de/) 用于展示系统发育树。

1.4 基因复制分析及计算Ka、Ks、Ka/Ks

使用MCScanX[18]软件鉴定海岛棉基因复制事件，Tbtools软件提取并展示GH9基因复制。使用ParaAT[19]实现复制基因对的密码子序列比对，使用KaKs_Calculator[20]计算非同义替换率(Ka)、同义替换率(Ks)以及二者比值(Ka/Ks)。

1.5 GH9基因在棉纤维中的表达分析

海岛棉品种Pima90-53、海7124(Hai7124)，陆地棉品种HY405、农大棉13号(ND13)、中棉所8号(CCRI8)、农大601(ND601)用于转录组测序，所用材料均由本课题组保存。测序文库分别构建自开花当天的胚珠(0 days post-anthesis, DPA)以及5、10、15、20、25、30 DPA的纤维。RPKM值(reads per kilobase of exon model per million mapped reads，每百万reads中来自某基因每千碱基长度的reads数)用于评估基因表达量。计算log2(1+RPKM),并用在线工具iTOL绘制热图。设计A、D亚基因组同源基因通用引物(表1)，利用实时荧光定量PCR方法确认GbGH9s在棉纤维发育过程中的表达模式，内参为组蛋白基因H3。

1.6 启动子顺式作用元件分析

JASPAR (http://jaspar.genereg.net/) 网站下载转录因子结合位点的位置频率矩阵，被选择的转录因子涉及植物细胞壁发育以及脱落酸、细胞分裂素、生长素、赤霉素、乙烯、油菜素内酯信号。提取海岛棉GH9基因起始密码子ATG上游2 000 bp基因组序列，基于JASPAR矩阵，使用FIMO[21]软件鉴定转录因子结合位点，设定阈值p<1E-5。

1.7 群体构建、功能标记开发、基因分型以及表型数据分析

利用陆地棉品种CCRI8作为受体亲本，海岛棉品种Pima90-53作为供体亲本，构建包含167个株系的BC3F5群体；基于A、D亚基因组序列差异，设计SNP引物GH9B6_R用以区分GH9B6与D亚基因组同源基因GH9B25，基于海岛棉和陆地棉基因组序列差异，设计SNP引物GbGH9B6_F和GhGH9B6_F用以区分直系同源基因GbGH9B6和GhGH9B6，在引物倒数第3位制造错配用以提高扩增条带区分度(表1)；使用两组引物(组合1：GbGH9B6_F和GH9B6_R；组合2：GhGH9B6_F和GH9B6_R)分别对167个株系进行基因分型，只在组合1有扩增条带的株系此位点为海岛棉纯合型GbGH9B6/GbGH9B6，只在组合2有扩增条带的株系此位点为陆地棉纯合型GhGH9B6/GhGH9B6，在组合1和2都有扩增条带的株系此位点为杂合型GbGH9B6/GhGH9B6；根据基因分型结果把群体分为3个亚群，比较亚群间的纤维品质差异，表型数据来自两个不同的生态环境(河北保定和新疆轮台)，数据分析使用GraphPad Prism软件，显著性P由双尾t检验计算。

表1 本研究中所用引物序列Table 1 Primer sequences in this study

2 结果与分析

2.1 海岛棉GH9基因家族成员鉴定及系统进化分析

基于海岛棉品种Hai7124的参考基因组，共鉴定53个具有糖苷水解酶GH9完整结构域的基因(表2)。基于陆地棉遗传标准系TM-1参考基因组共鉴定52个GH9基因。为了明晰家族成员进化关系，构建包含海岛棉、陆地棉、拟南芥[1]、水稻[8]、杨树[12]GH9基因的系统发育树(图1)。根据聚类结果和已有分类信息，53个海岛棉GH9基因被分为3类，包括8个A亚族基因(命名为GbGH9A1～GbGH9A8)、37个B亚族基因(命名为GbGH9B1～GbGH9B37)以及8个C亚族基因(命名为GbGH9C1～GbGH9C8)。此外，A亚族还包括8个陆地棉基因、3个拟南芥基因、3个水稻基因、4个杨树基因；B亚族包括37个陆地棉基因、19个拟南芥基因、18个水稻基因、18个杨树基因；C亚族包括7个陆地棉基因、3个拟南芥基因、4个水稻基因、3个杨树基因。可见，A、B、C三个亚族形成于单子叶植物和双子叶植物分化之前，而且GH9基因家族在植物中具有较高的保守性。

表2 53个海岛棉GH9蛋白的理化性质Table 2 Characteristics of 53 GH9 proteins in Gossypium barbadense

图1 海岛棉、陆地棉、拟南芥、水稻、杨树GH9基因的系统发育树Fig.1 Phylogenetic tree of GH9 genes from Gossypium barbadense, Gossypium hirsutum, Arabidopsis thaliana, Oryza sativa and Populus trichocarpa

2.2 海岛棉GH9家族成员基因结构和保守基序分析

GbGH9Cs的基因结构相对保守(图2A)，糖苷水解酶结构域分布在6或7个外显子上；B亚族成员基因结构变异幅度较大，可能对应其复杂多变的功能，糖苷水解酶结构域分布在4至8个外显子上；A亚族糖苷水解酶结构域主要分布在6个外显子上，然而GbGH9A4和GbGH9A8仅含有2个外显子，推测其最初来自反转录转座。通过分析海岛棉GH9家族成员保守基序(图2B)，发现15个被鉴定的保守基序全部集中在糖苷水解酶结构域范围内，表明此结构域在海岛棉中相当保守。糖苷水解酶结构域之外，A亚族成员N端通常具有跨膜域(图2C)，为膜结合蛋白；B亚族中22个成员被预测具有信号肽，可能为分泌蛋白(图2D和表2)；C亚族成员相对保守，全部具备信号肽和CBM49结构(图2E)。

A：基因结构；B：保守基序；C：A亚族成员结构域；D：B亚族成员结构域；E：C亚族成员结构域。A: Gene structure; B: Conserved motif; C: Domain of GH9A; D: Domain of GH9B; E: Domain of GH9C.图2 海岛棉GH9基因家族成员结构特征Fig.2 Structure characteristics of GbGH9s

2.3 海岛棉GH9基因染色体定位分析

根据基因位置信息，使用MapChart软件绘制GbGH9s的染色体分布图(图3)。发现53个基因零散地分布在除A04、A09、D04、D09之外的22条染色体上。A、B、C亚族基因均未形成基因簇，可见串联复制在海岛棉GH9基因家族扩张中作用较小。此外，GbGH9s在A、D亚基因组间无明显偏好性。

注：GbGH9As、 GbGH9Bs和GbGH9Cs基因分别用绿色、黑色和蓝色标注。Note: GbGH9As, GbGH9Bs and GbGH9Cs are colored in green, black and blue, respectively.图3 海岛棉GH9基因染色体分布Fig.3 Chromosomal distribution of GbGH9s

2.4 海岛棉GH9基因家族基因复制分析

为研究海岛棉GH9基因家族扩张，使用MCScanX软件鉴定基因复制事件。A亚基因组中发现9对片段复制/全基因组复制基因，涉及13个GbGH9s；D亚基因组中发现12对片段复制/全基因组复制基因，涉及12个GbGH9s(图4)。此外，GbGH9s中未鉴定到串联复制基因，可见片段复制/全基因组复制是海岛棉GH9基因家族扩张的主要驱动力。为评估复制发生时间，计算上述基因对的同义替换率(Ks)，发现2个区分明显的Ks值范围：0.44～0.73之间涉及15个GbGH9Bs和6个GbGH9Cs，对应棉属特有多倍化事件(约16.6 MYA)[22]；1.33～2.20之间涉及11个GbGH9Bs，对应双子叶植物共有的古六倍体化事件(约130.8 MYA)[22]。此外，计算复制基因对的Ka/Ks从而评估其受到的选择压力，发现所有复制基因对的Ka/Ks均远小于1，其中GbGH9Bs所属基因对的Ka/Ks范围为0.09～0.30，GbGH9Cs所属基因对的Ka/Ks范围为0.09～0.14(数据未列出)。综上，古六倍体化事件推动了B亚族的扩张，随后棉属特有多倍化事件推动了B亚族和C亚族的扩张；海岛棉GH9基因家族经历了严格的选择约束，其中C亚族可能在进化上更保守。

注：绿线代表A亚基因组片段复制基因对，蓝线代表D亚基因组片段复制基因对，橙线代表A亚基因组与D亚基因组同源基因对。数值代表复制基因对的同义替换率Ks。Note: Green line indicates segmental duplication within At-subgenome, blue line indicates segmental duplication within Dt-subgenome, orange line links homologous genes between At- and Dt-subgenome. The numbers represent the Ks values of duplicated gene pairs.图4 海岛棉GH9基因家族片段复制事件Fig.4 Segmental duplication events among GbGH9s

2.5 海岛棉GH9基因家族进化历程

为明确海岛棉GH9基因家族的进化历程，基于单子叶植物水稻及双子叶植物葡萄、大豆、毛果杨、拟南芥、黄麻、可可、榴莲、海岛棉GH9家族成员蛋白序列，构建系统发育树分析其进化关系。如图5所示，181个B亚族成员聚类成13组(Ⅰ～ⅩⅢ)，其中进化枝Ⅰ、Ⅲ、Ⅴ、Ⅵ、Ⅶ、Ⅷ、Ⅸ、Ⅹ、Ⅺ产生于单子叶和双子叶植物分化之前，而进化枝Ⅱ、Ⅳ产生于双子叶植物共有的古六倍体化事件。由进化枝扩张收缩情况推断，古六倍体化事件之后，可可由于未发生新的全基因组加倍事件[23]，可能维持了较为初始的GH9基因数量。海岛棉拥有全部13个进化枝，A、D亚基因组未分化前棉属经历多倍化事件[22]，造成进化枝Ⅱ、Ⅳ、Ⅹ、Ⅺ扩张。如图6所示，32个C亚族成员聚类成2组，C-Ⅱ在进化上可能更活跃，在本研究所涉及的大部分物种中扩张，其中海岛棉C-Ⅱ的扩张由棉属特有多倍化事件推动；35个A亚族成员聚类成3组，进化枝A-Ⅲ稳定存在于全部用于研究的物种并且趋于扩张，突显其重要性，而进化枝A-Ⅰ和A-Ⅱ仅存在于部分物种，海岛棉在进化中丢失A-Ⅱ。此外，近乎全部海岛棉GH9基因家族进化枝都拥有相同数量的A和D亚基因组成员。通过同源比对，发现无论是A、D亚基因组成员间，还是海岛棉和陆地棉成员间都存在高度匹配的系统进化关系，可见海岛棉和陆地棉GH9基因家族形成于棉种分化之前，并且在进化中高度保守。

A：水稻、葡萄、大豆、毛果杨、拟南芥、黄麻、可可、榴莲、海岛棉GH9基因家族B亚族成员的系统发育树，蓝线和红线分别指示古六倍体化事件和棉属特有多倍化事件引起的基因复制；B：植物GH9基因家族B亚族成员分类。A: Phylogenetic relationships of GH9Bs from Oryza sativa, Vitis vinifera, Glycine max, Populus trichocarpa, Arabidopsis thaliana, Corchorus capsularis, Theobroma cacao, Durio zibethinus and Gossypium barbadense. Blue and red lines indicate segmental duplications from the paleo-hexaploidization event and Gossypium-specific polyploidization, respectively; B: Clades of GH9Bs.图5 植物GH9基因家族B亚族成员进化历程Fig.5 Evolutionary history of GH9Bs

2.6 GbGH9s在棉纤维发育过程中的表达分析

考虑到植物GH9基因与细胞壁发育密切相关，利用海岛棉品种Pima90-53和Hai7124的RNA-seq数据分析GbGH9s在开花当天的胚珠以及5、10、15、20、25、30 DPA纤维中的表达谱。如图7所示，A亚族成员GbGH9A2、GbGH9A3、GbGH9A6和GbGH9A7在棉纤维发育各时期尤其是次生壁加厚期(20～30 DPA)高表达，可能与纤维素的生物合成相关；B亚族成员GbGH9B6和GbGH9B25拥有与上述A亚族基因相似的表达模式，可能同样涉及纤维素合成，或者涉及微纤丝修饰；C亚族成员GbGH9C1、GbGH9C3、GbGH9C4、GbGH9C6、GbGH9C7、GbGH9C8在棉纤维起始期(0 DPA)和伸长期(5～15 DPA)高表达，可能涉及初生壁的形成与重塑；B亚族成员GbGH9B13、GbGH9B23、GbGH9B32主要在棉纤维起始期高表达，可以作为改良棉纤维产量性状的候选基因。此外，海岛棉GH9家族在A和D亚基因组上的对应成员具有基本一致的表达模式，进一步表明其在进化中高度保守。为确认GbGH9s在棉纤维发育过程中的表达模式，使用实时荧光定量PCR方法检测GbGH9A2/GbGH9A6、GbGH9A3/GbGH9A7、GbGH9C1/GbGH9C6、GbGH9C3/GbGH9C7、GbGH9C4/GbGH9C8、GbGH9B6/GbGH9B25、GbGH9B13/GbGH9B32在Hai7124开花当天的胚珠以及10和25 DPA的纤维中的表达水平，其结果(图8)与RNA-seq数据相符合。

注：转录组数据来自开花当天的胚珠以及5、10、15、20、25、30 DPA的纤维。Note: RNA-seq data come from ovules of 0 DPA, and fibers of 5, 10, 15, 20, 25 and 30 DPA.图7 海岛棉GH9基因在棉纤维中的表达模式Fig.7 Expression patterns of GbGH9s in cotton fibers

图8 qRT-PCR验证GbGH9s在棉纤维发育中的表达模式Fig.8 Expression pattern validation of GbGH9s during cotton fiber development by qRT-PCR

此外，通过分析4个陆地棉品种HY405、ND13、CCRI8和ND601的RNA-seq数据，发现GhGH9s在棉纤维发育过程中的表达模式类似于GbGH9s。然而GbGH9A3、GbGH9A7、GbGH9B6在次生壁加厚期的表达高于其直系同源基因GhGH9A3、GhGH9A7、GhGH9B6；GbGH9B14在起始期的表达高于其直系同源基因GhGH9B14；GbGH9B27、GbGH9C1、GbGH9C4在伸长期的表达低于其直系同源基因GhGH9B27、GhGH9C1、GhGH9C4，上述直系同源基因的差异表达可能在一定程度上造就了海岛棉和陆地棉纤维品质的差异。进一步利用CottonFGD[24]中陆地棉遗传标准系TM-1转录组数据[25]分析GhGH9s在不同组织中的表达，发现GhGH9As和GhGH9Cs在用于分析的组织(根、茎、叶、副萼、花托、花瓣、雄蕊、雌蕊、胚珠、纤维)中具有更广泛的表达，推测其功能保守，而GhGH9Bs具有更强的组织特异性，对应其复杂多变的功能。

2.7 转录调控分析

众所周知，转录因子调控基因表达，为了解析GbGH9s在棉纤维中的表达机制，基于207个JASPAR位置频率矩阵，挖掘GbGH9s启动子上涉及植物激素信号或细胞壁发育的转录因子结合位点。在海岛棉纤维中，B亚族成员GbGH9B6和GbGH9B25高表达，而相邻进化枝上的GbGH9B8、GbGH9B20和GbGH9B27基本不表达。如图9A所示，GbGH9B6和GbGH9B25启动子上鉴定到植物细胞壁发育(MA1045.1)、生长素(MA1278.1、MA0565.2)和油菜素内酯(MA0964.1)相关转录因子结合位点，也鉴定到乙烯相关转录因子结合位点，涉及多个AP2/ERF (APETALA2/ethylene-responsive factor) 类转录因子DREB亚家族A-4/A-5亚组成员。GbGH9B8、GbGH9B20和GbGH9B27启动子上鉴定到细胞分裂素(MA1391.1、MA0945.1、MA1210.1)和脱落酸(MA1210.1、MA0570.2、MA1326.1)相关转录因子结合位点。通常，生长素、油菜素内酯、乙烯可以促进棉纤维生长发育，而细胞分裂素、脱落酸抑制棉纤维生长发育，因此GbGH9B6/GbGH9B25与GbGH9B8/GbGH9B20/GbGH9B27在棉纤维中的差异表达可能来源于上述顺式作用元件的差异。C亚族大部分成员在棉纤维起始期和伸长期高表达，而GbGH9C2和GbGH9C5基本不表达，分析发现GbGH9C2和GbGH9C5启动子上缺乏植物细胞壁发育(MA0990.1、MA1375.1)和生长素(MA1278.1)相关转录因子结合位点，也缺乏乙烯相关转录因子结合位点，涉及BPC(Basic Pentacysteine)和AP2/ERF家族的多个成员(图9B)。此外，DREB亚家族A-4亚组多个成员的靶序列在棉纤维高表达基因GbGH9A2、GbGH9A3、GbGH9A6和GbGH9A7的启动子上被鉴定出来，而未见于低表达(或不表达)基因GbGH9A1、GbGH9A4、GbGH9A5和GbGH9A8的启动子上。

2.8 GbGH9B6功能分析

通过分析海岛棉和陆地棉纤维的RNA-seq数据，发现GbGH9B6在次生壁加厚期的表达明显高于其直系同源基因GhGH9B6(图10A)。为研究GbGH9B6功能，设计用于区分GbGH9B6与GhGH9B6的SNP引物(图10B)，对拥有167个株系的BC3F5群体(供体亲本Pima90-53；受体亲本CCRI8)进行基因分型(图10C)，结果鉴定到17个GbGH9B6/GbGH9B6型株系、6个GbGH9B6/GhGH9B6型株系以及144个GhGH9B6/GhGH9B6型株系。分析群体表型数据，发现无论是在河北保定(E1)还是在新疆轮台(E2)，GbGH9B6/GbGH9B6型亚群的棉纤维强度都显著高于GhGH9B6/GhGH9B6型亚群(图10D)，而两个亚群的棉纤维长度和马克隆值并无显著变化，可见GbGH9B6可以作为改良棉纤维品质性状的候选基因。

A：GbGH9Bs启动子上转录因子结合位点预测；B：GbGH9Cs启动子上转录因子结合位点预测。基于JASPAR矩阵，使用FIMO软件鉴定潜在位点，设定阈值p<1E-5，括号中为产生JASPAR矩阵的转录因子。A: Identification of TFBS in the promoter regions of GbGH9Bs; B: Identification of TFBS in the promoter regions of GbGH9Cs. FIMO is used to identify potential sites based on JASPAR matrices with a threshold of p-value<1E-5, and the transcription factors producing JASPAR matrices are shown in brackets.图9 GbGH9Bs和GbGH9Cs的转录调控分析Fig.9 Transcriptional regulation of GbGH9Bs and GbGH9Cs

A：GH9B6在海岛棉和陆地棉纤维中的表达差异；B：SNP引物设计；C：基因分型；D：不同基因型株系棉纤维强度比较。A: Differential expression of GH9B6 in cotton fibers between Gossypium barbadense and Gossypium hirsutum; B: Design of SNP primers; C: Genotyping of the BC3F5 population; D: Box plot for fiber strength, based on different genotypes of GH9B6.图10 GbGH9B6功能分析Fig.10 Functional characterization of GbGH9B6

3 讨论

植物GH9基因家族古老且保守，分为A、B、C三个亚家族，且B亚族成员最多[1, 8, 12, 26]。本研究共鉴定到53个海岛棉GH9基因，包括8个A亚族基因，37个B亚族基因以及8个C亚族基因。分析基因复制事件发现，25个GbGH9s涉及片段复制/全基因组复制，未鉴定到串联重复基因，其结果类似于杨树[12]及玉米[26]GH9基因家族。计算复制基因对的同义替换率，发现2个明显区分的Ks范围1.33～2.20和0.44～0.73，分别对应古六倍体化事件(Ks范围1.5～1.9)[22]及棉属特有多倍化事件(Ks范围0.4～0.6)[22]，可见多倍化事件是海岛棉GH9基因家族扩张的主要推动力。不同于海岛棉，在拟南芥GH9基因家族[1]B亚族中观察到串联重复形成的基因簇AtGH9B9/10/11/12和AtGH9B16/17/18，使其进化枝Ⅲ和Ⅹ扩张(图5)。古六倍体化事件(γ)之后，拟南芥又经历了两次全基因组复制(α和β)[27]，然而其B亚族成员相较于其他物种并没有明显增多，原因可能是多倍化后染色体的大量断裂、融合，造成进化枝Ⅰ、Ⅷ、Ⅻ、ⅩⅢ丢失。在B亚族中，进化枝Ⅻ和ⅩⅢ最初可能由串联复制产生，分别来自Ⅻ和ⅩⅢ的成员可以构成基因簇(Tc03v2_t013060与Tc03v2_t013070、Glyma.06G277100与Glyma.06G277200、Glyma.12G004700与Glyma.12G004800、PtrGH9B8与PtrGH9B9)。然而，海岛棉Ⅻ和ⅩⅢ进化枝上的成员位于不同染色体，可能原因是基因簇先后经历片段复制/全基因组复制和基因丢失。在A亚族和C亚族中，未观察到明显的基因簇，推断其扩张来源于片段复制/全基因组复制或者转座，古六倍体化事件之后，大豆[28]、杨树[29]、拟南芥[27]、黄麻[30]、榴莲[31-32]、棉花[22]都经历了一到两轮全基因组加倍事件，对应A-Ⅲ和C-Ⅱ进化枝的扩张，而葡萄[33]和可可[23]未经历新的多倍化事件，因此保持了初始基因数量(图6)。由此可见，多倍化事件在植物尤其是棉花GH9基因家族扩张过程中起到了重要作用。

棉纤维的发育按照时间顺序可以分为四个不同但重叠的阶段：起始期、伸长期、次生壁加厚期和成熟期[34]。伸长阶段，棉纤维细胞的初生壁主要包括纤维素、木葡聚糖、木聚糖、果胶多糖和蛋白质等成分，其纤维素含量大约占初生壁干物质的20%～25%；次生壁增厚期，纤维细胞进入活跃的纤维素合成期，成熟棉纤维中的纤维素含量高达90%以上。GH9A/KORRIGAN涉及纤维素合成，可能的机制包括切除谷甾醇糖苷引物[35]、辅助微纤丝聚合[36]、构成纤维素合酶复合体[3]等。在拟南芥中，GH9A1/KOR1突变影响细胞扩展[2]及次生细胞壁形成[37]。在杨树中抑制PtrKOR1表达影响次生细胞壁增厚[5]。在棉花中下调GhKOR1表达影响棉纤维长度和结晶纤维素含量[6]。本研究中,GbGH9A2、GbGH9A3、GbGH9A6和GbGH9A7在棉纤维伸长和次生壁加厚期高表达(图7)，可能涉及纤维素生物合成，并且与AtGH9A1、PtrKOR1、GhKOR1(即本研究中的GhGH9A2)位于同一进化枝(图1)，可以作为改良棉纤维长度和强度的候选基因。细胞壁纤维素结晶度和细胞扩展速率密切相关[38]。在拟南芥中，过表达PtGH9C2提高细胞壁纤维素结晶度，进而降低植株高度和莲座叶大小；相反，沉默AtGH9C2降低细胞壁纤维素结晶度，增加植株高度和莲座叶大小[39]。棉纤维起始和伸长期优势表达基因GbGH9C1、GbGH9C3、GbGH9C4、GbGH9C6、GbGH9C7、GbGH9C8与AtGH9C2高度同源，可能具有相似功能，并且GbGH9C1和GbGH9C4在伸长期的表达低于其直系同源基因GhGH9C1和GhGH9C4，推测这些基因影响纤维素结晶度，对棉纤维发育具有重要作用，同时高结晶度也会一定程度上抑制棉纤维伸长。海岛棉C-Ⅱ进化枝的扩张完全来自棉属特有多倍化(图6)，并且GbGH9C7位于棉纤维长度相关QTL簇Clust_LEN_24_2范围内[40]，进一步突显C-Ⅱ成员的重要性及研究价值。GbGH9B6和GbGH9B25在棉纤维中的表达模式类似于GbGH9A2、GbGH9A3、GbGH9A6和GbGH9A7(图7)，可能涉及纤维素合成，且GbGH9B6和GbGH9B25高度同源于AtGH9B5(图1)，可能影响纤维素结晶度[39]，推测这两个基因可以用于改良棉纤维长度和强度。GbGH9B25位于长度相关QTL簇Clust_LEN_19_3范围内[40]，GbGH9B6定位于多个棉纤维品质QTL，包括长度(FB-FUqQtlc05_1b)和强度(BB-FMtQtlc05_1h)[41]，此外，通过回交将GbGH9B6导入陆地棉可提高棉纤维强度(图10)，表明其具有较好的应用价值。综上所述，GbGH9s是棉纤维品质改良的重要候选基因，也是研究棉纤维发育的重要基因资源，值得科研工作者关注。

从海岛棉全基因组中鉴定出53个GH9基因，可以分为A、B和C三个亚族。分析基因复制发现，多倍化是该家族扩张的主要驱动力，古六倍体化事件推动B亚族扩张，而棉属特有多倍化事件推动B亚族和C亚族的扩张。分析GbGH9s在棉纤维发育过程中的表达模式，筛选出多个品质改良候选基因，并通过BC3F5群体验证了GbGH9B6具有改良棉纤维强度的潜力。本研究为后续深入研究该家族基因奠定了基础，对于揭示棉纤维品质形成的遗传机制及分子育种具有一定的理论意义和应用价值。