本发明涉及林木育种技术领域，具体涉及一种旱柳耐盐枢纽基因的发掘方法。
背景技术：
耐盐基因通过调控蛋白作用，从而调控林木的耐盐表型。
传统的耐盐基因发掘是通过同源基因序列对比，或者用转录组测序比较等方法，这些方法无法精确找到林木耐盐的枢纽基因。同时考虑到林木自身的生长周期长，遗传背景复杂，杂合度高，研究基础薄弱，组学数据噪声和稀疏性较高等特性，传统方法局限性较大，不能适应林木耐盐基因发掘的需要。
旱柳(salixmatsudana)是杨柳科柳属重要树种，原产我国东北部，又被称为中国柳，在土壤重金属修复、盐碱地绿化、生物燃料等方面发挥着重要的作用，已广泛用于我国东部沿海造林。异速理论是一种用来研究生物有机体结构和过程与大小之间关系的科学，基于异速模型可以精确发掘旱柳耐盐枢纽基因，培育出更耐盐的旱柳新品种，为沿海滩涂生态建设做出贡献。而采用传统的耐盐柳树育种无法精准的找到耐盐基因。
技术实现要素：
为克服现有技术所存在的缺陷，现提供一种旱柳耐盐枢纽基因的发掘方法，以解决传统的耐盐柳树育种无法精准的找到耐盐基因的问题。
为实现上述目的，提供一种旱柳耐盐枢纽基因的发掘方法，包括：
将耐盐旱柳与敏盐旱柳杂交得到子代f1群体；
将所述子代f1群体划分为耐盐型子代和敏盐型子代；
在不同时长的盐胁迫下，分别对所述耐盐型子代和所述敏盐型子代的根系进行转录组测序得到差异表达基因序列；
基于所述耐盐型子代构建高密度遗传图谱，进行耐盐性qtls定位并计算得到相应的qtl区段；
将所述qtl区段、所述差异表达基因序列与柳树全基因组比对，分析获得多个耐盐候选基因；
在所述耐盐型子代受到不同时长的盐胁迫后，测定所述耐盐候选基因的多维数据，所述多维数据包括所述耐盐候选基因的表达量、所述耐盐型子代的耐盐表型指标和耐盐生理指标；
基于所述多维数据，建立异速生长模型；
基于所述异速生长模型，构建常微分方程组；
基于所述常微分方程组，构建所述耐盐基因的基因调控网络；
通过对所述常微分方程组结合变量选择，并基于所述基因调控网络的中心度特征，计算多个所述耐盐候选基因之间的互作程度以确定所述基因调控网络中的耐盐枢纽基因。
进一步的，所述将所述子代f1群体划分为耐盐型子代和敏盐型子代的步骤包括：
将所述子代f1群体的第一茎段放置于浓度为50mm的nacl溶液进行生根培养；
在生根培养两周后，测定所述第一茎段的耐盐表型指标；
基于耐盐表型指标，将所述子代f1群体划分为耐盐型子代和敏盐型子代。
进一步的，所述耐盐表型指标包括生根数、总根长、最长根长、最粗根直径、根鲜重、根干重。
进一步的，所述分别对所述耐盐型子代和所述敏盐型子代的根系进行转录组测序得到差异表达基因序列的步骤包括：
将分别所述耐盐型子代及所述敏盐型子代的第二茎段放置于清水中进行生根培养；
在生根培养两周后，将所述第二茎段放置于浓度为200mm的nacl溶液中进行盐胁迫；
在所述第二茎段盐胁迫4小时、8小时和12小时时，分别对所述第二茎段的根系进行转录组测序得到差异表达基因序列。
进一步的，所述异速生长模型为：
式中，n为耐盐候选基因的数量；
m为耐盐型子代的数量；
α为平衡常数项；
β为生长指数；
x为在一时间点，同一耐盐型子代中所有的耐盐候选基因的表达量的总和；
y为在一时间点，同一耐盐候选基因在所有的耐盐型子代中的表达量的总和。
进一步的，所述常微分方程组为：
其中，m为耐盐候选基因的表达量。
进一步的，还包括：
将所述耐盐枢纽基因导入柳树幼苗中，并检测所述柳树幼苗的耐盐生理指标；
基于所述柳树幼苗的耐盐生理指标，于所述耐盐枢纽基因中确定耐盐关键枢纽基因。
进一步的，所述将所述耐盐枢纽基因导入柳树幼苗中，并检测所述柳树幼苗的耐盐生理指标的步骤包括：
构建所述耐盐枢纽基因的植物表达载体；
将所述植物表达载体导入柳树幼苗中；
在所述柳树幼苗受不同时长的盐胁迫时，测定所述柳树幼苗的耐盐生理指标。
本发明的有益效果在于，本发明的旱柳耐盐枢纽基因的发掘方法可以大量的、精准的、快速的发掘出旱柳耐盐枢纽基因以解析柳树的耐盐机理，提高了耐盐柳树的选育效率，同时本发明的旱柳耐盐枢纽基因的发掘方法可以运用在其它耐盐林木上，更好地选育耐盐林木新品种。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。
本发明提供了一种旱柳耐盐枢纽基因的发掘方法，包括：
s1：将耐盐旱柳与敏盐旱柳杂交得到子代f1群体。
在本实施例中，敏盐旱柳为仅能在含盐量2‰以下的环境中正常生长的旱柳。耐盐旱柳为能在含盐量4‰以上的环境中生长的旱柳。
具体的，以敏盐旱柳为母本，以耐盐旱柳为父本。
在本实施例中，敏盐旱柳为沿江柳，耐盐旱柳为9901柳。
s2：基于柳树耐盐表型指标，将子代f1群体判定为耐盐型子代和敏盐型子代。
具体的，步骤s2包括：
步骤s21、将子代f1群体的第一茎段放置于浓度为50mm的nacl溶液进行生根培养。
在本实施例中，随机选择子代f1群体中的400株子代f1，截取400株子代f1的第一茎段，并将400株子代f1的第一茎段放置于浓度为50mm的nacl溶液进行生根培养。
步骤s22、在生根培养两周后，测定第一茎段的耐盐表型指标。
在400株子代f1的第一茎段的生根培养两周后，测定第一茎段的根系的耐盐表型指标。
柳树的耐盐表型指标包括但不限于：生根数、总根长、最长根长、最粗根直径、根鲜重、根干重。
具体的，生根数为茎段经培养后的生成的根的数量。
总根长为茎段的所有生成的根的总长度。
最长根长为茎段培养后所生成的最长根的长度。
最粗根直径为茎段培养后所生成的最粗根的最大直径。
根鲜重为茎段培养后所生成的所有的鲜活根的总重量。
根干重为茎段培养后所生成的所有的根经烘干后的总重量。
步骤s23、基于耐盐标准和敏盐标准，将子代f1群体中的子代f1分别判定为耐盐型子代或敏盐型子代。
具体的，耐盐标准为各耐盐表型指标值均大于所测耐盐表型指标的平均值。
敏盐标准为植株死亡或有耐盐表型指标小于所测耐盐表型指标的平均值。
依据耐盐标准和敏盐标准，分别将400株子代f1判定为耐盐型子代或敏盐型子代。
s3：在不同时长的盐胁迫下，分别对耐盐型子代和敏盐型子代的根系进行转录组测序得到差异表达基因序列。
在本实施例中，分别对判定为耐盐型子代和敏盐型子代的第一茎段的根系提取rna进行转录组测序分析得到差异表达基因序列(degs)。
具体的，步骤s3包括：
步骤s31、将随机选取30株的耐盐型子代、30株的敏盐型子代。截取耐盐型子代、敏盐型子代的第二茎段。将第二茎段放置于清水中进行生根培养。
步骤s32、在第二茎段生根培养两周后，将第二茎段放置于浓度为200mm的nacl溶液中进行盐胁迫。
步骤s33、在不同时长的盐胁迫下，分别对第二茎段的根系提取rna进行转录组测序得到差异表达基因序列。
具体的，在不同时长的盐胁迫下，即第二茎段在盐胁迫4小时、8小时和12小时时，分别对第二茎段的根系提取rna进行转录组测序得到差异表达基因序列(即，将30株耐盐型子代、30株敏盐型子代每种各10株组成共三种生物学重复混合样，在4、8、12小时盐胁迫后，将耐盐型子代、敏盐型子代两两比对相同时长、不同耐盐材料的degs，在三个时长均相同的degs为最终差异表达基因序列(degs))。
使用blast软件将转录组测序分析得到的差异表达基因序列degs序列与go、kegg、nr、swiss-prot、cog以及eggnog基因数据库进行比对，获得差异表达基因序列degs在上述的6个基因数据库的注释信息。
s4：基于耐盐型子代构建高密度遗传图谱，进行耐盐性qtls定位并计算得到相应的qtl区段。
具体的，利用400株子代f1群体基于snp标记构建高密度遗传图谱，定位与旱柳的根系测量指标相关的的耐盐性qtls，通过计算qtls的lod值与贡献率分析旱柳根系的耐盐性重要的qtl区段。
s5：将qtl区段、差异表达基因序列与柳树全基因组比对，分析获得多个耐盐候选基因。
将精细定位后的耐盐型子代的根系重要的qtl区段、耐盐型子代的根系转录组分析出的degs比对到旱柳全基因组上，分析获得旱柳根系耐盐候选基因群，即多个耐盐候选基因。
具体的，通过qtl区段上得到一批基因，转录组的degs注释得到一批基因，再将这些基因与全基因组进行对比，三者重合的基因确定为耐盐候选基因。
s6：在所述耐盐型子代受到不同时长的盐胁迫后，测定耐盐候选基因的多维数据。多维数据包括耐盐候选基因的相对表达量、耐盐型子代的根系耐盐表型指标和耐盐生理指标。
具体的，步骤s6包括：
步骤s61、截取随机选取的30株耐盐子代的第三茎段。将第三茎段放置在50mmnacl溶液中生根，动态观测培养0、3、6、9、12、15天的根系的耐盐表型指标，同时测定耐盐候选基因表达量。
步骤s62、截取随机选取的30株耐盐子代的第四茎段。将第四茎段放置清水中生根。在清水中生根两周后，将生根后的第四茎段放置于浓度为200mm的nacl溶液中进行0小时、6小时、12、24小时和48小时的盐胁迫，其中将6小时、12小时、24小时和48小时的盐胁迫设定为待测组，将0小时的盐胁迫设定为对照组。
在第三茎段受到盐胁迫0小时、6小时、12小时、12小时、24小时和48小时时，分别测定待测组和对照组的第四茎段的所述耐盐候选基因的表达量，同时测定耐盐型子代的第四茎段的耐盐生理指标。
耐盐生理指标包括电导率、丙二醛、脯氨酸、sod(超氧化物歧化酶superoxidedismutase)、pod(过氧化物酶peroxidase)、cat(过氧化氢酶catalase)、h2o2(过氧化氢hydrogenperoxide)及可溶性糖含量等。
具体的，耐盐候选基因的表达量数据，按照f＝2-δδct公式计算。
在本实施例中，耐盐候选基因的相对表达量(f)的计算公式为：
f＝2-[(待测组目的基因平均ct值-待测组内参平均ct值)-对照组目的基因平均ct值-对照组内参基因平均ct值)]。
在耐盐候选基因的相对表达量的计算公式中，2的上标的具体内容为：-[(待测组目的基因平均ct值-待测组内参基因平均ct值)-(对照组目的基因平均ct值-对照组内参基因平均ct值)]。
待测组为6、12、24、48小时的盐胁迫时的基因表达量，对照组为0小时的基因表达量。
s7：基于多维数据，即耐盐候选基因的相对表达量、耐盐生理指标和耐盐表型指标，将多维数据在某一基因或者表型指标进行标准化，基于所有指标的总和，建立异速生长模型。
以基因表达量列表为例，在时间点t的基因表达量可表示为：
其中xm表示在第t个时间点耐盐型子代旱柳个体m的所有的耐盐候选基因的基因表达量的总和，因此在一个时间点，基于一个耐盐候选基因在所有耐盐型子代旱柳个体的表达量的总和ynm与总的表达水平xm，即部分与整体的关系，构建如下的异速生长模型：
①式中，n为耐盐候选基因的数量；
m为耐盐型子代的数量；
α为平衡常数项；
β为生长指数；
α和β为未知参数，可通过非线性最小二乘法估计获得。其中β大于0一般为基因表达趋势上升，反之则为下降。
x为在一个时间点，同一耐盐型子代中所有的耐盐候选基因的表达量的总和；
y为在一个时间点，同一耐盐候选基因在所有的耐盐型子代中的表达量的总和。
s8：基于异速生长模型，构建常微分方程组。
首先，完成差异基因的筛选，然后基于异速生长模型，针对每个耐盐候选差异基因的基因表达表达量，通过非线性最小二乘法估计异速生长模型的两个关键参数α和β(见①式)。利用估计的参数，最终获得耐盐候选基因的连续的表达量。
差异基因的分析筛选过程基于rpackagemasigpro完成。
具体的，以p-value<0.05以及拟合优度r2>0.8作为筛选基准。
其次，使用k均值聚类算法(k-meansclusteringalgorithm)完成差异基因的表达量的聚类分析，最佳的聚类数量通过间隔统计量(gapstatistic)的方法确定。
最后，基于多维变量构建常微分方程组，然后进行变量选择。考虑变量m(m＝1,...,p)与其它变量的互作，变量可以表示为基因的表达量、耐盐生理指标等因子。考虑后续建立的耐盐候选基因的调控网络的稀疏性，采取分组最小角回归算法(grouplasso)进行变量选择。
构建的常微分方程组为：
其中，m为耐盐候选基因的表达量。
s9：基于常微分方程组构建耐盐基因的基因调控网络。
基于每个耐盐候选基因i的基因表达量的相对变化是一个包含所有耐盐候选基因的连接函数，分别表示为f1,...,fi。利用连接函数刻画所有耐盐候选基因的互作效应以形成耐盐基因的基因调控网络。
s10：通过常微分方程组结合变量选择，并基于基因调控网络的中心度特征，计算多个耐盐候选基因之间的互作程度以确定基因调控网络中的耐盐枢纽基因。
通过变量选择可获得耐盐候选基因之间是否存在互作，然后根据稀疏的基因互作网络构建高维微分方程，利用最大似然估计或者非线性最小二乘法结合平滑方法完成高维微分方程的参数估计，进而量化耐盐候选基因的互作程度。最后通过连通性、紧密性、介数、偏心率、特征向量等核心-边缘特性指标完成核心枢纽的筛选，具体的临界值可根据基因的数量和网络的复杂程度进行设定。
核心-边缘特性包括连通性(connectivity)、紧密性(closeness)、介数(betweenness)、偏心率(eccentricity)、特征向量(eigenvector)以及pagerank算法。
具体的，连通性是指一个网络内一个基因直接连接其它基因的数量。
紧密性是指网络内一个基因与其它基因的紧密程度。
紧密性以下式计算量化：
②式中，v表示网络内的耐盐候选基因的数量，dist(i，j)表示网络内耐盐候选基因i与耐盐候选基因j之间的最小距离。
介数反映耐盐候选基因作为“桥梁”的连系作用的重要程度。介数以下式计算量化：
③式中，gjk是连接耐盐候选基因j和耐盐候选基因k之间的最短路径数，gjk(i)是耐盐候选基因i位于最短路径gjk上的数量。
偏心率用于描述从其它的耐盐候选基因到达一个耐盐候选基因的最长距离。偏心率以下式计算量化：
特征向量用于衡量调控网络内的耐盐候选基因的重要性。当前的耐盐候选基因的重要性取决于邻居耐盐候选基因的重要性。对于耐盐候选基因i，其特征向量中心分数与其直接相连的耐盐候选基因总的积分成比例，以下式计算量化：
⑤式中λ是一个常数满足方程ax＝λx，aij表示调控网络内的耐盐候选基因之间是否存在连接，n是总的耐盐候选基因的数量。
pagerank算法对于耐盐候选基因i，pagerank中心分数定义为：
⑥式中d为阻尼因数，通常设置为0.85。
s11：通过瞬时转化验证耐盐枢纽基因的功能，并基于耐盐枢纽基因的功能及耐盐枢纽基因在基因调控网络中的上下关系，确定耐盐关键枢纽基因。
具体的，步骤s11包括：
s111、将耐盐枢纽基因导入柳树幼苗中，并检测柳树幼苗的耐盐生理指标。
具体的，包括：构建耐盐枢纽基因的植物表达载体。
通过植物表达载体将耐盐枢纽基于导入柳树幼苗中。
在柳树幼苗受不同时长的盐胁迫时，测定柳树幼苗的耐盐生理指标。
将导入耐盐枢纽基于的柳树幼苗放置于浓度为200mm的nacl溶液进行盐胁迫处理48小时、72小时和96小时，通过生化染色和酶活分析等方法分析柳树幼苗耐盐表型，测定确定耐盐生理指标，分析外源基因对其耐盐性的影响。
s112、基于柳树幼苗的耐盐生理指标，于耐盐枢纽基因中确定耐盐关键枢纽基因。
具体实施流程：
构建微分方程：
利用转录组测序技术(rna-seq)测定耐盐型子代，在盐胁迫处理条件下t个时间点或者样本的耐盐候选基因表达量的数据，共获得n个耐盐候选基因或者功能模块(耐盐生理指标)。让g1(t),....,gi(t)表示基因i在时间点或者样本t的基因表达量，t＝t0,...,t,i＝1,...,n。n个基因及其互作的动态变化可以使用系统常微分方程(ordinarydifferentialequations,ode)进行描述，具体表示如下：
该方程组为前述的常微分方程组的细化，参数ω代表未知参数，通过估计未知参数可以量化基因之间的互作程度。
其中，每个基因i的相对变化是一个包含所有基因的连接函数，分别表示为f1,...,fi；每个连接函数依赖于参数ω＝{ωj1,...,ωji}，利用这些参数能够刻画(表达)所有耐盐候选基因的互作效应以形成耐盐基因的基因调控网络。基于四阶龙格库塔法的最大似然估计以及非线性最小二乘法估计完成参数估计。
确定显著互作的耐盐枢纽基因
基于基因调控网络的稀疏性特性以及生物系统的生物学内涵，将高维变量选择技术grouplasso和adaptivegrouplasso应用于上述常微分方程组中每个基因所对应的ode初步确定存在显著互作的耐盐枢纽基因。
基于基因调控网络特征筛选耐盐关键枢纽基因
使用基因调控网络的中心度特征对基因调控网络中的节点进行评估。其中，中心度特征包括连通性、紧密性、介数、偏心率、特征向量以及pagerank。利用以上的中心度特征获取基因调控网络中每个耐盐候选基因与其它的耐盐候选基因的连接数、互作程度等指标，依据网络拓扑的复杂程度，通过排列检测确定各个指标的临界值，最终筛选确定出耐盐关键枢纽基因。
本发明的旱柳耐盐枢纽基因的发掘方法可以大量的、精准的、快速的发掘出旱柳耐盐枢纽基因以解析旱柳的耐盐机理，提高了耐盐柳树的选育的效率，同时本发明的旱柳耐盐枢纽基因的发掘方法可以运用在其它耐盐林木上，更好地选育耐盐林木新品种。
最后应说明的是：显然，上述实施例仅仅是为清楚地说明本申请所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引申出的显而易见的变化或变动仍处于本申请型的保护范围之中。