2025, Vol. 55
2. 中国海洋大学海洋生物遗传与育种教育部重点实验室，山东 青岛 266003;
3. 山东省海水种业重点实验室，山东 青岛 266003;
4. 青岛蓝色种业研究院，山东 青岛 266237
条斑紫菜(Neopyropia yezoensis)富含蛋白质、维生素和矿物质，经济价值高，是中国主要栽培的经济海藻之一。然而，随着栽培规模的不断扩大，条斑紫菜出现了种质退化的现象，严重影响了条斑紫菜栽培产业的稳定发展，因此加快良种的培育和推广具有重要的意义。目前，条斑紫菜育种还是以常规育种为主，主要采用选择育种、诱变育种或者二者结合的技术方法[1]，这些方法培育出的条斑紫菜良种为产业的发展提供了重要支撑作用。然而，育种周期较长、目标不明确和遗传资源挖掘不充分等因素限制了条斑紫菜的育种效率。与传统的育种方法相比，基于分子标记技术的育种技术避免了盲目选择和长时间的试验种植的问题，节省了时间成本，该技术已成为了一种重要的育种手段[2]。分子育种方法通常利用高通量测序技术对大规模种质资源或育种群体进行单核苷酸多态性(Single nucleotide polymorphisms, SNPs)鉴定，根据SNP与目标性状的相关性，选择具有有利基因组组合的个体进行培育。目前，该方法在动、植物育种工作中得到了广泛的应用[3-4]。
全基因组关联分析(Genome wide association study, GWAS)以连锁不平衡(Linkage disequilib-rium, LD)理论为基础，通过检测群体全基因组范围内的变异位点，选用合适的模型统计每个分子标记与目标性状之间的关联情况，筛选出与目标性状显著相关的标记位点[5]。1996年，Risch等[6]对基因组关联研究的理论进行了探讨，并预见到群体分层可能导致假阳性结果的问题。2002年，Ozaki等[7]利用94名心肌梗死患者及其基因组上的65 000个单核苷酸多态性，报道了心肌梗死风险的遗传关联。2005年，Doebley等[8]利用95份拟南芥(Arabidopsis thaliana)的全基因组多态性数据对开花时间和病原体抗性进行了关联分析，证明了GWAS在拟南芥中的可行性，这也是GWAS在人类疾病领域之外的第一次应用[9]。随着高通量技术的发展和统计算法的不断完善，关联分析开始在畜禽品种[10]和农作物品种[11]中有了广泛应用。目前，GWAS已在多个重要农作物品种中展开，如玉米、小麦、水稻、大豆和棉花等[12-13]。这些研究明确了多项农艺性状(如株高、籽粒数和产量)相关的基因组区域，为作物的遗传改良提供了理论基础。近年来，GWAS逐渐应用于水产生物重要经济性状的遗传变异研究，为水产动、植物的育种工作提供了可靠的分子标记。在对鲶鱼耐热性和柱状病抗性的遗传学研究中，通过GWAS发现了3个与耐热性明显关联的SNP位点，以及多个与柱状病抗性相关的数量性状基因座(Quantitative trait locus, QTL)[14-15]。Feng等[16]借助GWAS方法，揭示了68个与龙须菜(Gracilariopsis lemaneiformis)生长速率相关的SNPs，并进一步筛选出79个可能影响龙须菜生长速率的潜在候选基因。随着测序技术的不断发展，测序成本显著降低，这为GWAS在水产生物中的应用提供了有利条件。
在以往的条斑紫菜研究中，尚未找到能够有效指示生长性状的分子标记，这制约了对条斑紫菜遗传机制的深入探究。因此，通过GWAS开发与生长性状显著相关的分子标记，将为阐释条斑紫菜生长性状的遗传基础和开展育种工作提供更精准的指导。
1 实验方法 1.1 材料的采集和处理野生条斑紫菜叶状体采集于山东省烟台市蓬莱区(120.8°E, 37.8°N)和威海市环翠区(122.7°E, 37.3°N)两个地点。将采集的叶状体使用毛笔刷去表面附着物，放置于0.7%碘化钾溶液中浸泡5 min后，用海水冲洗干净。从336个叶状体的不同部位切下面积约0.5 cm2的小块，并剁碎成约50个细胞大小的碎片后，转移到预先标记编号的7个六孔板中(每个孔移入来自8个叶状体的碎片)。加入含有PES培养基的灭菌海水，置于光照培养箱中培养，以进行单孢子放散。培养条件如下：光周期L∶D=12 h∶12 h，光照强度20 μmol·m-2·s-1，温度15 ℃，添加终浓度为1.5 μg/mL的GeO2以抑制硅藻。培养7 d后，首次更换配置海水，之后每3 d更换一次。培养30 d后，转移至2 L的充气瓶中培养(每1个六孔板对应1个充气瓶)，培养条件如下：光周期L∶D=12 h∶12 h，光照强度40 μmol·m-2·s-1，温度10 ℃；培养40 d后，每个充气瓶中保留约25个样本；培养至55 d时，收集到169个条斑紫菜样本(见附表 1)。
1.2 表型测定和数据处理表型的测量包括叶长(Thallus length, L)、叶宽(Thallus width, W)、鲜质量(Fresh mass, FM)和干质量(Dry mass, DM)。由于条斑紫菜叶状体形状不规则，本研究定义叶长是指从叶片基部到顶端最长部分的长度，叶宽是指叶片最宽部分的长度，测量时，将叶片完全铺展(放置刻度尺进行标注)，使用拍照装置固定在叶片正上方60 cm处进行拍照，使用ImageJ软件进行计算，精确至0.01 cm；鲜质量是指叶片表面除水后的质量，测量时，用吸水纸将叶片表面的水分擦干，迅速放入电子天平中进行测量，精确至0.01 mg；干质量是指细胞除去全部自由水后的质量，测量时，将叶片放置于65 ℃烘箱中烘8 h，直至质量保持不变，再放入电子天平中测量，精确至0.01 mg。使用IBM SPSS Statistics 27.0在线软件(https://www.ibm.com)对叶长、叶宽、鲜质量和干质量4个表型数据进行描述性统计分析和相关性分析。为观测表型数据的分布情况，使用R软件包的PerformanceAnalytics模块绘制表型数据的频率分布图。
1.3 DNA提取和测序文库的构建随机挑选出130个纯化培养的单一基因型样本(见附表 1)，使用DNA提取试剂盒(购自北京诺贝莱生物科技有限公司)提取DNA，并进行质量检测：(1)使用NanoPhotometerⓇ分光光度计(IMPLEN, CA, USA)检测DNA纯度；(2)使用QubitⓇ 3.0 Flurometer(Life Technologies, CA, USA)检测DNA浓度；(3)使用琼脂糖凝胶电泳(胶浓度1%、电压120 V、时间45 min)检测DNA的纯度和完整性，DNA质检合格的样品进行测序文库的构建。测序文库构建完成后，使用Qubit进行初步定量，再利用qPCR测定文库的有效浓度，以确保测序文库的质量，最后使用HiSeq测序平台进行双末端测序。
1.4 测序质量检测二代测序会得到包含有带接头的、低质量的reads，为确保数据分析的质量，需要对原始数据(Raw data)进行过滤。使用fastp软件(Version: 0.20.0)[17]对原始数据过滤的步骤如下: (1)删除接头序列(Adapter)；(2)删除缺失超过10个碱基的reads；(3)删除碱基测序质量较低的reads。得到clean reads后，使用FastQC软件(Version: 0.12.1)[18]统计单个样本的测序质量信息，再使用MultiQC软件(Version: 1.19)[19]对全部样本的测序质量信息进行汇总。在NCBI网站(https://www.ncbi.nlm.nih.gov/)下载条斑紫菜的参考基因组序列文件(GenBank: GCA_009829735.1)。使用bwa软件(Version: 0.7.17-r1188, mem算法)构建参考基因组的索引文件，并将质控后的clean reads映射到参考基因组上。使用Samtools软件(Version: 1.6)[20]将比对后的clean reads进行排序，使用Bamdst软件(Version: 1.0.9)统计比对率、平均测序深度和覆盖度等信息。
1.5 SNP位点的检测、筛选和注释使用GATK软件(Version: 4.4.0)[21]构建条斑紫菜参考基因组的索引文件，标记并去除测序文库中的重复reads。单独为每个样本生成个体的变异信息文件(gvcf文件)，然后对群体进行联合基因型分型。提取出群体的SNP位点，参考GATK给出的过滤条件，通过硬过滤(Hard-filter)的方式过滤掉低质量的SNP位点。为保证后续分析的准确性，使用VCFtools软件(Version: 0.1.16)[22]再次对变异位点进行过滤，删除测序质量分数小于30、测序深度小于3、缺失率大于10%、次要等位基因频率小于5%的位点。在NCBI网站下载条斑紫菜的参考基因组注释文件(gtf格式)，使用gtfToGenePred软件将基因组注释文件转换成GenePred格式。结合条斑紫菜的参考基因组序列文件，使用ANNOVAR软件(Version: Sun, 7 Jun 2020)[23]构建条斑紫菜的基因组注释数据库，并对高质量的SNP位点进行注释，统计变异位点在基因组区域的分布、外显子区域的变异功能和变异类型(转换和颠换)。
1.6 群体结构和亲缘关系分析通过IQ-TREE软件[24]使用极大似然法(Maximum likelihood estimation, MLE)[25]对基因分型数据构建系统进化树，随后使用ITOL在线网站进行处理。使用GCTA软件计算基因分型数据的遗传关系矩阵，并对GRM矩阵进行特征值分解，得到特征值和对应的特征向量。选择前两大特征值对应的特征向量作为主成分，使用R软件作图。使用Admixture软件[26]，通过最大似然法计算材料来源于不同亚群的概率。通过设置亚群数量K的范围为1~10，并计算每个K值的交叉验证标准误差(Cross-validation standard error, CVSE)，找到CVSE最小时的K值作为最佳K值。使用GCTA软件根据基因分型数据计算GRM矩阵(一种基于基因组数据的kinship矩阵)，最后使用R语言绘制亲缘关系热图。使用PopLDdecay软件(Version: 3.42)[27]计算烟台采集群体(YT)、威海采集群体(WH)和总群体SNP标记之间的LD系数，并根据SNP之间的距离和r2值，来计算LD衰减率(LD decay)。
1.7 全基因组关联分析使用GEMMA软件[28]计算亲缘关系矩阵，将PCA分析中前30个主成分和亲缘关系矩阵作为协变量同时引入混合线性模型，对条斑紫菜的叶长、叶宽、干质量和鲜质量四个性状进行全基因组关联分析。模型如下：
| $ y=X \beta+S \alpha+u+e \text { 。} $ | (1) |
式中：y为表型观测值向量，作为模型的响应变量；X为涵盖截距项及群体结构主成分的协变量矩阵，其对应的固定效应系数向量为β，共同构成控制群体分层等系统性偏差的固定效应部分；S为当前待检验SNP的基因型向量，α为该SNP的加性效应值，二者共同量化目标SNP对表型的潜在关联效应；u为用于校正个体间遗传相关的随机效应向量；e为模型未能解释的残差向量，涵盖环境因素、测量误差等随机误差项。
根据关联分析结果中P值，使用R软件中CMplot包绘制曼哈顿图(Manhattan plot)和分位数-分位数(Quantile-quantile, Q-Q)图，曼哈顿图用于展示基因组上与性状显著相关的标记位点，Q-Q图用于判断关联分析中的假阳性结果是否得到较好的控制。根据Bonferroni校正的原理[29]，阈值设为1/SNP的标记数。通过Aspichueta等[30]的计算公式得到每个SNP位点的表型方差解释率(Phenotypic variance explained, PVE)。公式如下：
| $ \begin{gathered} R_{\mathrm{PVE}}=\left(2 \beta^2 F_{\mathrm{MA}}\left(1-F_{\mathrm{MA}}\right)\right) /\left(2 \beta^2 F_{\mathrm{MA}}\left(1-F_{\mathrm{MA}}\right)+\right. \\ \left.(E(\beta))^2 2 N F_{\mathrm{MA}}\left(1-F_{\mathrm{MA}}\right)\right) 。\end{gathered} $ | (2) |
式中：RPVE为表型变异解释率；β为GWAS中的effect值；FMA为SNP的最小等位基因频率；E(β)为β的标准误(Standard error，SE)；N为GWAS中的样本数。
最后，根据设定的阈值筛选与叶长、叶宽、鲜质量和干质量性状显著相关的SNP，参考条斑紫菜基因组的注释信息，选择显著SNP位点上、下游20 kB区域内的基因作为候选功能基因。
1.8 富集分析使用在线软件eggNOG-mapper(http://eggnog-mapper.embl.de/)对候选基因进行功能注释，并利用GO数据库、KEGG数据库和Omicshare在线软件对注释到的基因进行GO富集分析和KEGG富集分析，选择P值最小的前10个GO富集条目和前15个KEGG富集通路进行绘图。
2 实验结果 2.1 表型数据分析在本研究中，用于表型测量的材料为野外材料(基因型嵌合)放散单孢子后长成的叶片，因此每个用于测量的叶片均为单一基因型的叶片。条斑紫菜4个重要生长性状的统计结果显示(见表 1)，叶长和叶宽的平均值分别为5.73和0.95 cm，变异系数分别为36.13%和32.63%。鲜质量和干质量的平均值分别为15.84和2.61 mg，变异系数分别为50.95%和49.04%。各个性状的变异系数较大，这表明不同个体之间存在显著的表型差异，为后续的GWAS分析提供了条件。
|
表 1 生长性状统计分析 Table 1 Statistical analysis of the four growth traits |
对条斑紫菜生长性状之间的相关性进行评估，发现条斑紫菜叶长、叶宽、鲜质量和干质量之间存在极显著的相关性(见图 1)。其中鲜质量与干质量的相关性最强(0.95)，叶长与叶宽的相关性最弱(0.24)。拟合曲线显示，叶宽、鲜质量和干质量之间存在较强的线性关系。叶长、叶宽、鲜质量和干质量性状的分布直方图(见图 1A—D)显示，各个性状的分布近似符合正态分布，这表明它们具有数量性状的特征。
|
( A. 叶长分布特征(柱状图、核密度曲线、正态分布曲线)；B. 叶宽分布特征(柱状图、核密度曲线、正态分布曲线)；C. 鲜质量分布特征(柱状图、核密度曲线、正态分布曲线)；D. 干质量分布特征(柱状图、核密度曲线、正态分布曲线)；E. 叶宽与鲜质量的相关性(散点图+拟合线，r2=0.502)；F. 叶宽与干质量的相关性(散点图+拟合线，r2=0.552)；G. 干质量与鲜质量的相关性(散点图+拟合线，r2=0.899)；H. 各指标(鲜质量、干质量、叶长、叶宽)间的相关性热图，含相关性系数及显著性标记(: P＜0.01; : P＜0.001)。A. Distribution characteristics of thallus length, presented as bar charts, kernel density curves, and normal distribution curves; B. Distribution characteristics of thallus width, presented as bar charts, kernel density curves, and normal distribution curves; C. Distribution characteristics of fresh mass, presented as bar charts, kernel density curves, and normal distribution curves; D. Distribution characteristics of dry mass, presented as bar charts, kernel density curves, and normal distribution curves; E. Correlation between thallus width and fresh mass, shown as a scatter plot with a fitting line (r2=0.502); F. Correlation between thallus width and dry mass, shown as a scatter plot with a fitting line (r2=0.552); G. Correlation between dry mass and fresh mass, shown as a scatter plot with a fitting line (r2=0.899); H. Correlation heatmap of fresh mass, dry mass, thallus length, and thallus width, including correlation coefficients and significance markers (: P < 0.01; : P < 0.001). ) 图 1 生长性状相关性统计分析 Fig. 1 Correlation statistical analysis between the growth traits |
经过全基因组重测序共获得345.4 GB的数据量，测序数据中平均重复reads占比为25.92%。经过滤后，测序产出数据的Q30均大于90%，测序数据中每个碱基位置的平均质量分数均大于30，这表明测序的错误率小于0.001，说明测序质量高。统计结果显示，测序数据的reads比对率为39.33%，平均测序深度为5.27，平均覆盖度为70.83%。
2.3 SNP位点的检测、筛选和功能注释对重测序数据进行变异位点的筛选和过滤后，共得到了405 999个高质量的SNP位点。每个样本的SNPs位点缺失率小于10%，这表明这些个体能够检测出90%以上的变异位点。使用ANNOVAR软件对WGS基因分型数据中SNP位点注释的结果显示(见表 2)，263 516个SNP位点(64.9%)位于基因间区域，56 318个SNP位点(13.9%)位于转录起始位点上、下游1 kb的区域内，19 959个SNP位点(4.9%)位于5′非翻译区或(和)3′非翻译区，12 945个SNP位点(3.2%)位于内含子区域，139个SNP位点位于剪接位点的2 bp之内，53 077个SNP位点(13.1%)位于编码区。编码区SNP位点中，24 415个SNP位点为非同义突变，将改变氨基酸的编码类型；578个SNP位点将导致翻译过程提前终止或肽链异常延伸。SNP位点突变数量最多的类型为G∶C>A∶T转换，突变数量最少的类型是A∶T>T∶A颠换，转换与颠换数量之比为1.23(见图 2)。
|
|
表 2 SNPs注释结果 Table 2 SNPs annotation results |
|
( 纵坐标代表SNP位点的数目；横坐标代表SNP位点的突变类型(A∶T>G∶C表示A∶T碱基对转换为G∶C碱基对，G∶C>A∶T则相反；A∶T>T∶A表示A∶T碱基对颠换为T∶A碱基对，其他同为颠换)。The vertical axis represents the number of SNP loci, while the horizontal axis denotes the mutation type (A∶T>G∶C indicates a transition from A∶T to G∶C, and G∶C>A∶T represents the reverse; A∶T>T∶A indicates a transversion from A∶T to T∶A, with other cases following the same pattern) of SNP loci. ) 图 2 SNP位点突变类型 Fig. 2 The mutation types of SNP loci |
130个条斑紫菜样本基因型主成分分析的结果显示(见图 3A)，前两个主成分的方差贡献率分别为17.36%和8.47%，累计的方差贡献率为25.83%。测序群体可分为4个亚群，与地理区域划分不同，这可能是由于不同地理群体之间产生基因交流导致的。聚类分析显示(见图 3B)，不同亚群间存在明显的遗传差异，与主成分分析的结果一致，而且不同地理群体的遗传差异也非常明显。使用Admixture软件进行群体结构分析，绘制了不同K值下的CVSE变化折线图(见图 3C)和遗传结构图(见图 3D)。当K=3时，CVSE最小，表明130个条斑紫菜样本可分为3个亚群。当K=4时，CVSE与当K=3时接近，且群体分类的结果更接近PCA分析的结果。从遗传结构图看出，亚群之间具有不同的遗传背景，呈现出较为明显的差异。
|
( A. 选择前两个主成分因子绘制的主成分分析图，不同颜色的椭圆代表不同亚群，红色方框和绿色圆点代表不同的采集群体；B. 聚类分析图，圆圈上不同颜色代表不同的亚群，红色和绿色分支代表不同的采集群体；C. CVSE变化折线图，当K=3时，CVSE最小；D. 当K=3和K=4时的遗传结构图，不同颜色线段代表不同的亚群和采集群体。A. A PCA plot with the first two principal components, where differently colored ellipses denote different subgroups, red rectangles and green dots signify different collection populations; B. A cluster analysis plot, with different colors representing distinct subgroups, red and green branches representing different collection populations; C. A line graph of CVSE, showing that it is minimal when K=3; D. Genetic structure plots for K=3 and K=4, with different-colored segments indicating subgroups and collection populations. ) 图 3 群体结构分析 Fig. 3 Analysis of population structure |
使用GCTA软件分析了测序个体之间的亲缘关系，亲缘关系矩阵显示多数个体间的亲缘关系系数小于0.5，表明它们之间的亲缘关系较远，有利于后续的GWAS分析(见图 4A)。连锁不平衡分析的结果显示(见图 4B)，威海采集群体的LD衰减速率低于烟台采集群体，总群体的LD衰减速度快，r2的最大值为0.51，LD衰减距离较短(80 bp)。
|
( A. 亲缘关系热图，横坐标和纵坐标代表不同的样品，红色越深代表样品间的亲缘关系越近，反之则代表亲缘关系越远；B. LD衰减图，代表了不同采集群体的连锁不平衡衰减趋势。A. A heatmap of kinship, where horizontal and vertical coordinates represent different samples, with darker red indicating closer relationships and lighter red indicating more distant ones; B. An LD decay plot, illustrating the pattern of linkage disequilibrium (LD) declining trends for different collection populations. ) 图 4 亲缘关系分析与连锁不平衡分析 Fig. 4 Genetic relationship analysis and linkage disequilibrium analysis |
全基因组关联分析的曼哈顿图和Q-Q图(见图 5)显示，根据设定的基因组水平显著关联阈值-lgP>5.60，分别筛选到了67个与叶长相关的位点，17个与叶宽相关的位点，4个与鲜质量相关的位点，以及2个与干质量相关的位点。这些位点的表型解释率分别为16.1%~24.8%(叶长)、16.1%~20.7%(叶宽)、16.2%~17.5%(鲜质量)和18.1%(干质量)。值得注意的是，其中2个SNP(Chr1:18395778和Chr1:20680490)在鲜质量和干质量两个性状中同时被检测到，存在一因多效的现象。去除重复的SNP后，总共得到了88个位点。曼哈顿图显示，与叶长、叶宽显著相关的SNP分布在多条染色体上，而与鲜质量和干质量相关的SNP主要分布在CM020618.1染色体上。
|
( A1. 叶长曼哈顿图；A2. 叶长Q-Q图；B1. 叶宽曼哈顿图；B2. 叶宽Q-Q图；C1. 鲜质量曼哈顿图；C2. 鲜质量Q-Q图；D1干质量曼哈顿图；D2. 干质量Q-Q图。红色虚线代表筛选阈值(-lgP = 5.60)。A1. Manhattan plot of thallus length; A2. Q-Q plot of thallus length; B1. Manhattan of thallus width; B2. Q-Q plot of thallus width; C1. Manhattan plot of fresh mass; C2. Q-Q plot of fresh mass; D1. Manhattan plot of dry mass; D2. Q-Q plot of dry mass. The red dotted line represents the screening threshold (-lgP = 5.60). ) 图 5 条斑紫菜GWAS的结果 Fig. 5 Results of GWAS for Neopyropia yezoensis |
根据这些SNP位点，对其上、下游20 kb范围内的基因进行了筛选，得到175个候选基因，其中具有注释信息的候选基因有39个(见附表 2)。对于叶长、叶宽性状，分别有20和16个基因得到了注释；对于鲜质量性状，有2个基因得到了注释，其中基因SKIV2L2还与干质量性状相关。此外，一些常见的基因MPK17、SQD1和PCNA被认为与细胞代谢、信号传递和生长发育有关。
对具有注释信息的候选基因进行了富集分析。GO富集的结果显示，生物过程富集136项，分子功能富集32项，细胞组分富集98项。在生物过程分类上(见图 6A)，主要富集在小分子化合物生物合成(乙醇、蝶啶、硫胺素二磷酸等)、脱氧核糖核酸酶活性的正向调节和细胞凋亡负调控等过程；在分子功能分类上(见图 6B)，主要富集在RNA解旋酶活性、ATP酶活性和DNA N-糖基化酶活性等方面；在细胞成分分类上(见图 6C)，主要富集在高尔基体亚室、PCNA复合体和DNA聚合酶延长性因子复合体等组分中。KEGG富集结果(见图 6D)显示，候选基因富集于核苷酸切除修复、倍半萜类和三萜类的生物合成、错配修复、内吞作用和核黄素代谢等通路。
|
( A. 在生物学过程中富集到的P值前10小的GO条目；B. 在分子功能上富集到的P值前10小的GO条目；C. 在细胞成分中富集到的P值前10小的GO条目；D. 候选基因富集到的P值最小的前15个KEGG通路。A. The top 10 GO terms with the lowest P-values enriched in the biological process; B. The top 10 GO terms with the lowest P-values enriched in molecular function; C. The top 10 GO terms with the lowest P-values enriched in cellular component; D. The top 15 KEGG pathways with the lowest P-values enriched for candidate genes. ) 图 6 富集分析 Fig. 6 Enrichment analysis |
在自然条件下，条斑紫菜个体之间的基因交流频繁，而壳孢子的萌发过程会经过减数分裂，由此生成的叶状体可能是由不同基因型细胞组成的嵌合体，这对条斑紫菜的基因分型产生了影响。本研究利用条斑紫菜配子体无性生殖释放单孢子的特点进行纯化培养，构建了由单一基因型个体组成的群体，避免了嵌合体对GWAS研究的影响。
对条斑紫菜表型数据进行统计的结果显示，不同个体之间的表型差异明显，由于样本均在相同环境下培养，这些表型差异可能是由遗传因素导致的。不同基因的组合和变异可能会影响其形态、生长速度和色素合成等方面的特征，从而导致表型差异的出现。表型差异可以增加GWAS研究中关联信号的数量和质量，当表型差异较大时，基因与表型间的关联更容易被GWAS方法检测到，从而提高了GWAS研究的可靠性。
在育种中，表型相关性分析对于种质的选择和改良具有重要作用。通过分析不同表型之间的相关性，可以确定哪些性状与育种目标密切相关。例如，在农作物育种中，通过研究产量、品质、抗逆性和抗病性之间的相关性，可以确定相互影响的性状，从而选择合适的目标性状进行改良[31]。这种方法对那些难以测量的性状很有帮助，可以帮助育种者更准确地评估其潜在的表型特征。本研究发现，条斑紫菜叶长、叶宽、鲜质量和干质量这4种性状两两之间均存在显著的相关性，其中鲜质量和干质量的相关性系数为0.95，这表明两者的相关程度很高。此外，这些表型数据呈现连续性特征，且大致遵循正态分布，这与数量性状的特征一致。这些结果表明，本研究构建的GWAS群体具有相对丰富的遗传变异，为全面了解其性状奠定了基础。
3.2 GWAS统计模型的选择群体分层和亲缘关系是造成GWAS出现假阳性结果的重要原因[32]。群体分层是指群体中存在不同等位基因频率的亚群，这可能导致等位基因频率差异较大的位点与目标性状之间产生虚假关联，从而产生假阳性的结果[33]。因此，在进行GWAS之前，有必要对构建群体的群体结构进行分析。主成分分析将构建群体分为4个亚群，分群情况并没有严格按照地理区域划分，这可能是因为采集地点比较接近，不同群体之间产生了基因交流。主成分分析结果同系统发育树和种群结构分析结果一致，均表明在构建群体中存在群体分层。
亲缘关系较近的个体之间存在较高的遗传相似性，包括共享的遗传变异和基因频率。当进行GWAS时，可能会将基因的相似性误认为是基因与目标性状之间的关联，从而产生假阳性结果[5]。亲缘关系分析的结果显示，多数个体间的亲缘关系系数小于0.5，说明亲缘关系较远，这保证了样本之间有足够的遗传多样性。此外，个别样本间的亲缘关系较近，这时需要排除亲缘关系对GWAS的干扰。为减少群体分层和亲缘关系对GWAS结果的影响，本研究将混合线性模型(Mixed linear model association, MLMA)应用于关联分析，该模型将群体结构作为固定效应，亲缘关系作为随机效应，因此能够显著提高模型分析的统计效力[34]。同时，GEMMA软件以其强大的计算能力在大样本应用场景中独具优势，能显著加速模型的运行，从而大幅度提升工作效率[32]。
GWAS分析流程中，曼哈顿图和Q-Q图作为常用的可视化工具，在数据解读过程中发挥关键作用。曼哈顿图显示，关联位点的分布比较分散，这可能是因为生长性状容易受到多个微效基因的影响，导致多个关联位点被检测到。Q-Q图(见图 5A1—D1)显示，P值分布初始趋向于均匀分布的直线，这表明数据没有明显的偏离，而尾部区域的升高，则反映出存在显著的关联信号。从Q-Q图(见图 5A2—D2)可以看出，混合线性模型有效地减少了群体分层和亲缘关系造成的假阳性关联，提高了GWAS结果的准确性和可靠性。
3.3 SNP位点及候选基因调控机制的初步探讨在GWAS分析中，通常利用Bonferroni校正法来确定关联筛选的显著性阈值(即0.05/SNP总数)，然而该方法的筛选条件过于严苛，显著减少了关联位点的数量，所以本研究将显著性阈值设为1/SNP总数[35]。根据该阈值，筛选得到了88个与生长相关的显著性位点，其中与干质量相关的SNP位点同时与鲜质量显著相关，这表明干质量和鲜质量之间具有紧密的联系。对关联位点上、下游20 kb范围内的基因进行筛选，得到了175个候选基因，其中39个基因得到了注释，包括3个重要的功能基因MPK17、SQD1和PCNA[36]。
MPK17参与丝裂原活化蛋白激酶(Mitogen-activated protein kinase, MAPK)级联反应，在调控植物生长发育、生物和非生物应激反应以及生物信号传递中发挥着重要作用[37]。在水稻中，通过RNAi技术降低MPK17基因表达水平，导致节间长度、穗粒数、穗长、结实率和千粒质量均显著降低，这表明MPK17基因在水稻生长发育过程中起着重要作用[38]。SQD1是短链脱氢酶/还原酶(Short-chain dehydrogenases/reductases，SDR)超家族的成员，对于控制细胞代谢途径、转录和信号传导具有重要意义[39]。Sun等[40]对水稻的SQD1基因进行了RNA干扰和过表达，OsSQD1的沉默抑制了根、花粉和花粉囊的发育，降低了穗长、结实率及单株产量，而过表达材料的单株产量、千粒质量和种子宽度明显增加则表明SQD1能够影响水稻生长和发育。研究表明，细胞周期基因可能在调控植物生长和形态结构中发挥重要作用[41]。增殖细胞核抗原(Proliferating cell nuclear antigen, PCNA)是一种细胞周期蛋白，参与DNA复制和修复过程，具有激发DNA聚合酶δ的活性和与细胞周期调节剂p21相互作用的能力，在动、植物细胞增殖过程中起着关键作用[42]。GO富集分析中，在细胞成分分类上富集到PCNA复合体、PCNA-p21复合体和DNA聚合酶延长性因子复合体等。
在细胞分裂和增殖过程中，DNA损伤修复机制的精细调控对于维持细胞的正常生命周期至关重要[43]。至少有5种主要的DNA修复途径能够使细胞修复DNA损伤：碱基切除修复、核苷酸切除修复、错配修复、同源重组和非同源末端连接[44]。GO富集分析中，在生物过程上富集到脱氧核糖核酸酶活性的正向调节，在分子功能上富集到DNA N-糖基化酶活性的正向调节。在碱基切除修复中，主要通过DNA糖基化酶和AP核酸内切酶的联合作用去除损伤[45]。KEGG富集分析中，候选基因被富集到核苷酸切除修复和错配修复通路中。其中核苷酸切除修复是通用的DNA修复途径，它可以去除破坏DNA双链稳定性的大段DNA损伤[46]。而错配修复涉及多种细胞过程，包括微卫星稳定性、减数分裂重组、DNA损伤信号传导和细胞凋亡等[47]。
植物能够合成各种小分子化合物，包括植物生长繁殖所必需的初级代谢产物和生长发育调节中的次级代谢产物[48-49]。在GO富集分析中，生物过程分类主要富集在乙醇、蝶啶和硫胺素二磷酸等小分子化合物的生物合成过程中，这些小分子化合物在糖代谢和氨基酸代谢途径中发挥重要作用。此外，KEGG富集显示，候选基因富集到倍半萜类和三萜类的生物合成和核黄素代谢通路中。萜类化合物(如Thalianin)在拟南芥(Arabidopsis thaliana)中已被证明在调节根系生长和塑造健康的根际细菌群落方面发挥了重要作用[50-51]。核黄素确保了许多黄素蛋白的功能，包括脱氢酶、氧化酶、单加氧酶和还原酶，它们在线粒体电子传递链、脂肪酸β氧化、氧化还原稳态、柠檬酸循环、支链氨基酸分解代谢、染色质重塑、DNA修复、蛋白质折叠和细胞凋亡中起着关键作用[52]。
本研究的结果基于数据的计算和统计，因此为进一步验证结果的可靠性，本研究团队将考虑在未来构建不同的测试群体或者通过分子生物学实验进行验证，以提高结果的可重复性和可信度。
4 结语本研究对包含130株单一基因型条斑紫菜个体的群体进行全基因组重测序，从中筛选得到405 999个高质量的SNP位点，并通过GWAS分析获得了与重要生长性状相关的90个SNP位点和39个候选功能基因。这一结果为深入阐释条斑紫菜重要生长性状的遗传基础提供了有力支撑，同时也为条斑紫菜遗传育种工作提供了极具价值的基因资源。
|
|
附表 1 样本信息 Supplementary table 1 Sample information |
|
|
附表 2 与叶长、叶宽、鲜质量和干质量性状相关的候选基因 Supplementary table 2 Candidate genes significantly associated with thallus length, thallus width, fresh mass and dry mass traits |
| [1] |
周伟, 胡传明, 陆勤勤, 等. 条斑紫菜的种质创新与应用[J]. 广西科学院学报, 2021, 37(1): 46-52. Zhou W, Hu C M, Lu Q Q, et al. Germplasm innovation and application of Pyropia yezoensis[J]. Journal of Guangxi Academy of Sciences, 2021, 37(1): 46-52. (  0) |
| [2] |
Collard B C Y, Mackill D J. Marker-assisted selection: An approach for precision plant breeding in the twenty-first century[J]. Philosophical Transactions of the Royal Society B: Biological Sciences, 2007, 363(1491): 557-572. ( 0) |
| [3] |
Raina V S, Kour A, Chakravarty A K, et al. Marker-assisted selection vis-à-vis bull fertility: Coming full circle-a review[J]. Molecular Biology Reports, 2020, 47(11): 9123-9133. DOI:10.1007/s11033-020-05919-0 ( 0) |
| [4] |
Fujino K, Hirayama Y, Kaji R. Marker-assisted selection in rice breeding programs in Hokkaido[J]. Breeding Science, 2019, 69(3): 383-392. DOI:10.1270/jsbbs.19062 ( 0) |
| [5] |
Uffelmann E, Huang Q Q, Munung N S, et al. Genome-wide association studies[J]. Nature Reviews Methods Primers, 2021, 1(1): 59. DOI:10.1038/s43586-021-00056-9 ( 0) |
| [6] |
Risch N, Merikangas K. The future of genetic studies of complex human diseases[J]. Science, 1996, 273(5281): 1516-1517. DOI:10.1126/science.273.5281.1516 ( 0) |
| [7] |
Ozaki K, Ohnishi Y, Iida A, et al. Functional SNPs in the lymphotoxin-α gene that are associated with susceptibility to myocardial infarction[J]. Nature Genetics, 2002, 32(4): 650-654. DOI:10.1038/ng1047 ( 0) |
| [8] |
Doebley J, Aranzana M J, Kim S, et al. Genome-wide association mapping in Arabidopsis identifies previously known flowering time and pathogen resistance genes[J]. PLoS Genetics, 2005, 1(5): e60. DOI:10.1371/journal.pgen.0010060 ( 0) |
| [9] |
Tibbs Cortes L, Zhang Z, Yu J. Status and prospects of genome-wide association studies in plants[J]. The Plant Genome, 2021, 14(1): e20077. DOI:10.1002/tpg2.20077 ( 0) |
| [10] |
Aulchenko Y S, De Koning D J, Haley C. Genomewide rapid association using mixed model and regression: A fast and simple method for genomewide pedigree-based quantitative trait loci association analysis[J]. Genetics, 2007, 177(1): 577-585. DOI:10.1534/genetics.107.075614 ( 0) |
| [11] |
Beló A, Zheng P, Luck S, et al. Whole genome scan detects an allelic variant of fad2 associated with increased oleic acid levels in maize[J]. Molecular Genetics and Genomics, 2007, 279(1): 1-10. ( 0) |
| [12] |
Ersoz E S, Yu J, Buckler E S. Applications of linkage disequilibrium and association mapping in crop plants[M]// Tuberosa R Genomics-Assisted Crop Improvement: Vol 1: Genomics Approaches and Platforms. [s. l.]: Springer, 2007: 97-119.
( 0) |
| [13] |
Liu H J, Yan J. Crop genome-wide association study: A harvest of biological relevance[J]. The Plant Journal, 2018, 97(1): 8-18. ( 0) |
| [14] |
Jin Y, Zhou T, Geng X, et al. A genome-wide association study of heat stress-associated SNPs in catfish[J]. Animal Genetics, 2016, 48(2): 233-236. ( 0) |
| [15] |
Geng X, Sha J, Liu S, et al. A genome-wide association study in catfish reveals the presence of functional hubs of related genes within QTLs for columnaris disease resistance[J]. BMC Genomics, 2015, 16(1): 196. DOI:10.1186/s12864-015-1409-4 ( 0) |
| [16] |
Feng X, Xiao B, Jiang M, et al. Identification of candidate genes related to two economic traits using GWAS in Gracilariopsis lemaneiformis (Rhodophyta)[J]. Algal Research, 2023, 76: 103309. DOI:10.1016/j.algal.2023.103309 ( 0) |
| [17] |
Chen S F, Zhou Y Q, Chen Y R, et al. fastp: An ultra-fast all-in-one FASTQ preprocessor[J]. Bioinformatics, 2018, 34(17): 884-890. DOI:10.1093/bioinformatics/bty560 ( 0) |
| [18] |
de Sena Brandine G, Smith A D. Falco: High-speed FastQC emulation for quality control of sequencing data[J]. F1000 Research, 2019, 8: 1874. DOI:10.12688/f1000research.21142.1 ( 0) |
| [19] |
Ewels P, Magnusson M, Lundin S, et al. MultiQC: Summarize analysis results for multiple tools and samples in a single report[J]. Bioinformatics, 2016, 32(19): 3047-3048. DOI:10.1093/bioinformatics/btw354 ( 0) |
| [20] |
Li H, Handsaker B, Wysoker A, et al. The sequence alignment/map format and SAMtools[J]. Bioinformatics, 2009, 25(16): 2078-2079. DOI:10.1093/bioinformatics/btp352 ( 0) |
| [21] |
Van Der Auwera G A, Carneiro M O, Hartl C, et al. From fastQ data to high-confidence variant calls: The genome analysis toolkit best practices pipeline[J]. Current Protocols in Bioinformatics, 2013, 43(1110): 1-33. ( 0) |
| [22] |
Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools[J]. Bioinformatics, 2011, 27(15): 2156-2158. DOI:10.1093/bioinformatics/btr330 ( 0) |
| [23] |
Wang K, Li M, Hakonarson H. ANNOVAR: Functional annotation of genetic variants from high-throughput sequencing data[J]. Nucleic Acids Research, 2010, 38(16): e164. DOI:10.1093/nar/gkq603 ( 0) |
| [24] |
Nguyen L T, Schmidt H A, Von Haeseler A, et al. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[J]. Molecular Biology and Evolution, 2015, 32(1): 268-274. DOI:10.1093/molbev/msu300 ( 0) |
| [25] |
Nagy P, Szabó á, Váradi T, et al. Maximum likelihood estimation of FRET efficiency and its implications for distortions in pixelwise calculation of FRET in microscopy[J]. Cytometry Part A, 2014, 85(11): 942-952. DOI:10.1002/cyto.a.22518 ( 0) |
| [26] |
Dominguez Mantes A, Mas Montserrat D, Bustamante C D, et al. Neural ADMIXTURE for rapid genomic clustering[J]. Nature Computational Science, 2023, 3(7): 621-629. DOI:10.1038/s43588-023-00482-7 ( 0) |
| [27] |
Zhang C, Dong S S, Xu J Y, et al. PopLDdecay: A fast and effective tool for linkage disequilibrium decay analysis based on variant call format files[J]. Bioinformatics, 2019, 35(10): 1786-1788. DOI:10.1093/bioinformatics/bty875 ( 0) |
| [28] |
Vogt F, Shirsekar G, Weigel D, et al. vcf2gwas: Python API for comprehensive GWAS analysis using GEMMA[J]. Bioinformatics, 2022, 38(3): 839-840. DOI:10.1093/bioinformatics/btab710 ( 0) |
| [29] |
Curtin F, Schulz P. Multiple correlations and bonferroni's correction[J]. Biological Psychiatry, 1998, 44(8): 775-777. DOI:10.1016/S0006-3223(98)00043-2 ( 0) |
| [30] |
Aspichueta P, Shim H, Chasman D I, et al. A multivariate genome-wide association analysis of 10 LDL subfractions, and their response to statin treatment, in 1868 caucasians[J]. PLoS One, 2015, 10(4): e0120758. DOI:10.1371/journal.pone.0120758 ( 0) |
| [31] |
肖扬, 陈康, 丛倩倩, 等. 平菇品种的黄斑病抗病性及其与主要农艺性状相关性分析[J]. 食药用菌, 2021, 29(6): 509-512. Xiao Y, Chen K, Cong Q Q, et al. Analysis of resistance to yellow blotch disease of Pleurotus ostreatus and its correlation with the main agronomic characters[J]. Edible and Medicinal Mushrooms, 2021, 29(6): 509-512. ( 0) |
| [32] |
Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies[J]. Nature Genetics, 2012, 44(7): 821-824. DOI:10.1038/ng.2310 ( 0) |
| [33] |
Barsh G S, Sul J H, Martin L S, et al. Population structure in genetic studies: Confounding factors and mixed models[J]. PLoS Genetics, 2018, 14(12): e1007309. DOI:10.1371/journal.pgen.1007309 ( 0) |
| [34] |
Yu J, Pressoir G, Briggs W H, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness[J]. Nature Genetics, 2005, 38(2): 203-208. ( 0) |
| [35] |
Zhu F, Sun H, Jiang L, et al. Genome-wide association study for growth-related traits in golden pompano (Trachinotus ovatus)[J]. Aquaculture, 2023, 572: 739549. DOI:10.1016/j.aquaculture.2023.739549 ( 0) |
| [36] |
Wang D, Yu X, Xu K, et al. Pyropia yezoensis genome reveals diverse mechanisms of carbon acquisition in the intertidal environment[J]. Nature Communications, 2020, 11(1): 4028. DOI:10.1038/s41467-020-17689-1 ( 0) |
| [37] |
Zhang M, Su J, Zhang Y, et al. Conveying endogenous and exogenous signals: MAPK cascades in plant growth and defense[J]. Current Opinion in Plant Biology, 2018, 45: 1-10. ( 0) |
| [38] |
Zhu Z, Wang T, Lan J, et al. Rice MPK17 plays a negative role in the Xa21-mediated resistance against Xanthomonas oryzae pv. oryzae[J]. Rice, 2022, 15(1): 41. DOI:10.1186/s12284-022-00590-4 ( 0) |
| [39] |
Mulichak A M, Theisen M J, Essigmann B, et al. Crystal structure of SQD1, an enzyme involved in the biosynthesis of the plant sulfolipid headgroup donor UDP-sulfoquinovose[J]. Proceedings of the National Academy of Sciences, 1999, 96(23): 13097-13102. DOI:10.1073/pnas.96.23.13097 ( 0) |
| [40] |
Sun Y, Jain A, Xue Y, et al. OsSQD1 at the crossroads of phosphate and sulfur metabolism affects plant morphology and lipid composition in response to phosphate deprivation[J]. Plant, Cell & Environment, 2020, 43(7): 1669-1690. ( 0) |
| [41] |
Den Boer B G W, Murray J A H. Control of plant growth and development through manipulation of cell-cycle genes[J]. Current Opinion in Biotechnology, 2000, 11(2): 138-145. DOI:10.1016/S0958-1669(00)00072-0 ( 0) |
| [42] |
Strzalka W, Ziemienowicz A. Proliferating cell nuclear antigen (PCNA): A key factor in DNA replication and cell cycle regulation[J]. Annals of Botany, 2011, 107(7): 1127-1140. DOI:10.1093/aob/mcq243 ( 0) |
| [43] |
Jackson S P, Bartek J. The DNA-damage response in human biology and disease[J]. Nature, 2009, 461(7267): 1071-1078. DOI:10.1038/nature08467 ( 0) |
| [44] |
Chatterjee N, Walker G C. Mechanisms of DNA damage, repair, and mutagenesis[J]. Environmental and Molecular Mutagenesis, 2017, 58(5): 235-263. DOI:10.1002/em.22087 ( 0) |
| [45] |
Kow Y W. Repair of deaminated bases in DNA[J]. Free Radical Biology and Medicine, 2002, 33(7): 886-893. DOI:10.1016/S0891-5849(02)00902-4 ( 0) |
| [46] |
Krasikova Y, Rechkunova N, Lavrik O. Nucleotide excision repair: From molecular defects to neurological abnormalities[J]. International Journal of Molecular Sciences, 2021, 22(12): 6220. DOI:10.3390/ijms22126220 ( 0) |
| [47] |
Jiricny J. The multifaceted mismatch-repair system[J]. Nature Reviews Molecular Cell Biology, 2006, 7(5): 335-346. DOI:10.1038/nrm1907 ( 0) |
| [48] |
Fang C, Fernie A R, Luo J. Exploring the diversity of plant metabolism[J]. Trends in Plant Science, 2019, 24(1): 83-98. DOI:10.1016/j.tplants.2018.09.006 ( 0) |
| [49] |
Zhou X, Liu Z. Unlocking plant metabolic diversity: A (pan)-genomic view[J]. Plant Communications, 2022, 3(2): 100300. DOI:10.1016/j.xplc.2022.100300 ( 0) |
| [50] |
Huang A C, Jiang T, Liu Y X, et al. A specialized metabolic network selectively modulates Arabidopsis root microbiota[J]. Science, 2019, 364(6440): eaau6389. DOI:10.1126/science.aau6389 ( 0) |
| [51] |
Bai Y, Fernández-Calvo P, Ritter A, et al. Modulation of Arabidopsis root growth by specialized triterpenes[J]. New Phytologist, 2021, 230(1): 228-243. DOI:10.1111/nph.17144 ( 0) |
| [52] |
Shanti B, Joy Y L. Riboflavin metabolism: Role in mitochondrial function[J]. Journal of Translational Genetics and Genomics, 2020, 4(4): 285-306. ( 0) |
2. Key Laboratory of Marine Genetics and Breeding, Ministry of Education, Ocean University of China, Qingdao 266003, China;
3. Shandong Provincial Key Laboratory of Marine Seed Industry, Qingdao 266003, China;
4. Qingdao Institute of Blue Seed Industry, Qingdao 266237, China