基于轉錄組和代謝組解析川芎對鎘脅迫的響應機制

關鍵詞：川芎；鎘脅迫；轉錄組；代謝組

川芎是傘形科藁本屬植物川芎（Ligusticumchuanxiong Hort.）的干燥根莖。作為川產道地藥材，早在公元1153年（宋朝），四川都江堰一帶就有川芎的種植和藥用記錄[1]。目前，川芎的核心道地產區是四川盆地的彭州、什邡、都江堰、眉山等地，國內90% 以上的商品藥材都來源于此[2]。但近年來川芎在出口時曾多次因重金屬超標問題而被銷毀[3]，重金屬超標問題是川芎國際貿易壁壘的主要成因。川芎重金屬污染中鎘（Cd）污染尤為突出[4]，而栽培土壤生物及理化特性是導致川芎Cd污染的主要原因[5]。在川芎中，Cd主要富集在其根部以及根莖，Cd脅迫能抑制川芎葉片的光合速率，并顯著降低根莖中阿魏酸和川芎嗪含量[6]。因此，需要利用更多技術手段來控制中藥材的品質安全，使中藥材質量達到穩定、有效、安全、可控。

代謝組學利用氣相色譜- 質譜（gaschromatography-mass spectrometry，GC-MS）、液相色譜- 質譜（liquid chromatographytandem-massspectrometry，LC-MS/MS）等先進技術研究生物體內代謝產物的變化及其隨時間的演變[7]，能更好地理解植物受重金屬等逆境脅迫下的生物代謝過程。羅慶[8]利用代謝組學技術探究Cd、鉛（Pb）污染下東南景天根系抵御和富集重金屬的作用機制，發現車前可以通過調節糖類和氨基酸類等差異代謝物來響應土壤Cd 脅迫[9]。此外，轉錄組技術也被用來探究植物體Cd脅迫的應答機制。汪京超[10]鑒定出油菜體內與Cd吸收和代謝相關的轉運蛋白和功能基因；王書鳳[11]通過轉錄組學分析，篩選出Cd 累積高效油菜品種P78 和低效油菜品種P72，并闡明了Cd累積高、低效品種響應Cd脅迫的主要機制。

目前，對中藥材川芎Cd脅迫的研究主要集中在Cd對植物生長發育和生理生化的影響方面，在轉錄組測序分析方面的研究尚處于起步階段，而有關吸收或轉運Cd的基因也未發現，極大限制了耐Cd新品種川芎的選育。因此，本研究以川芎為試驗材料，采用超高效液相色譜質譜聯用（ultraperformance liquid chromatography tandem massspectrometry，UPLC-MS/MS）技術測定代謝組數據并解析川芎根莖代謝物的差異表達，結合轉錄組學技術對川芎中Cd轉運相關同源基因進行研究，從基因表達的角度分析川芎受Cd脅迫后的基因變化規律，揭示川芎響應Cd脅迫的防御機制，以期為川芎Cd污染耕地種植提供理論依據與科學支撐，為川芎資源的合理利用提供方向。

1 材料與方法

1.1 試驗材料

試驗材料來自四川省農業科學院經濟作物研究所篩選的川芎，于2021年9月9日栽種于成都市青白江區光明村溫室大棚內。該研究采用盆栽試驗，花盆規格為42 cm×28 cm×28 cm（上徑×下徑×高），盆栽基質按農田土∶草炭土∶珍珠巖∶河沙=25∶4∶1∶8的體積比混勻，每盆裝基質19 kg，每盆栽種6株川芎。以不添加Cd溶液的清潔土壤為對照（CK），設置1（Cd1）、3（Cd3）、6（Cd6）、10 mg·kg-1（Cd 10）共4 個Cd 脅迫處理水平。CdCl2·2.5H2O（析純試劑）以溶液形式施入土壤中充分混勻，平衡4周后栽入川芎。在2022年6月21日（即川芎收獲期），快速將川芎根莖挖出并沖去表面泥土后分成2份，分別用于轉錄組和代謝組紙標的測定。每個處理取3株長勢相同的川芎（分別用-1、-2、-3表示，如Cd 1的3株材料分別為Cd1-1、Cd1-2和Cd1-3），將其根莖切分、裝袋后放入液氮中備用。

1.2 試驗方法

1.2.1 代謝組測定及分析 ①樣品制備。將川芎根莖置于40 ℃烘箱中烘干，研磨成粉末后過100目篩。稱取50.00 mg的川芎粉末樣品，加入1.00 mL的提取液（提取液的組成為甲醇：乙腈∶水=2∶2∶1，體積比），維持其內標水平在20 mg·L?1。通過30 s的渦旋，保證樣品與提取液充分混勻。在45 Hz下添加鋼珠進行研磨處理，10 min后將樣品放入冰水浴中進行10 min的超聲處理。然后，?20 °C靜置1 h，4 °C、12 000 r·min?1離心15 min，取其上清液500 μL。通過真空濃縮器將分離的代謝產物干燥后加入160 μL提取液（乙腈與水的體積比為1∶1），復溶后再次超聲離心，得到120 μL的上清液。每個樣取10 μL混勻，制成質控（quality control，QC）樣本用于后續檢測[12]。

② 色譜條件。使用Waters Acquity I-ClassPLUS UPLC/Xevo G2-XS Q-TOF 超高效液相色譜質譜聯用儀（沃特世，美國）進行測定，每個樣品平行3次測量。色譜柱為Waters ACQUITITY UPLC高速鋼T3色譜柱（2.1 mm×100 mm，1.8 μm）；在負離子模式下進行測定，流動相A為0.1 %甲酸水溶液，流動相B為0.1 %甲酸乙腈；進樣體積設置為1 μL。

③ 質譜條件。采用Waters Xevo G2-XS Q-TOF高分辨率質譜儀（沃特世，美國）以MSe模式采集一次和二次質譜數據。低碰撞能量：2 V；離子源溫度：150 °C；高碰撞能量范圍：10～40 V；質譜掃描頻率：0.2 s；負離子模式下毛細管電壓：?1 500 V；反沖洗氣體流量：50 L·h?1；脫溶氣體流量：800 L·h?1。

④ 方法學考察。通過觀察QC樣品間內標峰高與空白樣品中內標峰條件的差異，判斷儀器是否穩定，檢測過程中是否存在殘留，結果顯示，內標在QC樣品中的保留時間和響應強度都比較穩定，儀器的數據采集穩定性很好，空白樣品中除內部標準品外，未檢測到明顯的峰。證明該方法可用于本次試驗。

⑤ 數據分析。采用MassLynx V4.2軟件獲取原始數據，利用Progenesis QI軟件對其進行峰提取、峰對等數據分析。利用已建立的Biomark自建庫以及METLIN數據庫進行鑒別，將理論片段鑒定和質量偏差控制在100 ppm以內。在此基礎上，將原始峰面積信息與總峰面積進行歸一化處理后，利用主成分分析和Spearman相關分析，對各組內樣品及質量控制樣品進行可重復性評價。在京都基因與基因組百科全書數據庫（kyotoencyclopedia of genes and genomes，KEGG）、人類代謝組數據庫（human metabolome database，HMDB）和脂質圖數據庫中搜索鑒定出的化合物以獲取分類和途徑信息。經過統計分析，確定不同組分之間的差異倍數（fold change，FC），并計算其差異顯著性P 值。使用正交偏最小二乘判別分析（orthogonal projections to latent structuresdiscriminant analysis，OPLS-DA）模型，篩選出具有顯著性的差異累積代謝物（differential accumulatedmetabolite， DAMs），以FCgt;1、Plt;0.05、變量投影重要性（variable importance in projection，VIP）gt;1作為篩選指標[13]。采用超幾何分布分析方法，分析代謝產物在KEGG通路中的富集情況。

1.2.2 轉錄組測定及分析 ① RNA提取。使用RNA提取試劑盒（酷來博，北京）對川芎根莖RNA進行提取；利用分光光度計NanoDrop 2000（賽默飛，美國）檢測提取的RNA純度（OD260/280）和含量；采用Agient2100/LabChip GX儀器（安捷倫，美國）進行高精度檢測以確保RNA完整性。

② cDNA文庫構建。樣本檢測合格后，利用含有Oligo（dT）的磁珠對真核mRNA進行富集，再添加Fragmentation Buffer 將mRNA 片段隨機打斷。以獲得的mRNA為模板，分別合成2條cDNA鏈，并對其進行純化。在此基礎上，利用AMPureXP beads進行片段篩選，得到240 bp左右的片段，經PCR擴增獲得cDNA文庫。

③ 文庫質控。使用實時熒光定量PCR（qRTPCR）定量文庫有效水平（文庫有效水平大于2 nmol·L-1），通過質量檢測后，利用IlluminaNovaSeq 6000（因美納，美國）進行高質量測序，獲得Raw data。經過篩選最終得到高質量的Cleandata[14]。

④ 序列拼接及功能注釋。將所得到的CleanData采用Trinity軟件進行序列組裝拼接，利用DeBruijn 圖識別各個片段集合中的轉錄本序列（transcript）。Unigene序列通過與NR、Swiss-Prot、COG、KOG、eggNOG4.5、KEGG 數據庫進行比對，整合數據庫分析得到基因的注釋結果。

⑤ 差異表達分析。用FPKM 分析基因表達量，與Unigene庫進行對比。使用DESeq2軟件進行差異分析，對得到的差異表達基因進行KEGG富集分析。

2 結果與分析

2.1 代謝組學分析

2.1.1 差異累積代謝物統計 由圖1A可知，在不同Cd 處理水平的川芎根莖中，共標注和定量了1 238種差異累積代謝物（DAMs），根據其分子結構可分為氨基酸、有機酸、脂肪酸、酮類、糖類、酰胺類和其他共7大類。與CK相比，Cd1、Cd3、Cd6、Cd10處理的DAMs總數分別為283、252、172、531個，其中上調代謝物數目均高于下調代謝物數目。Cd10處理的DAMs 最多，包括346 個上調代謝物和185個下調代謝產物； Cd6處理的DAMs最少，包括111個上調代謝物和61個下調代謝物。

對DAMs進行分析，由圖1B韋恩圖可知，CKvs Cd1 與CK vs Cd10 共有DAMs 數量最多，為134個，而4組共有差異代謝物僅為14個。這14種差異代謝物包括泛酸鹽、1-O-阿魏-D葡萄糖、多邊形內酯、多球殼菌素、阿非迪霉素、亞油酸等。

2.1.2 不同處理間差異代謝物分析 利用正交偏最小二乘判別分析（OPLS-DA）發現（圖2），CK的根莖樣品可以從Cd處理組中顯著分離出來。其主要貢獻來自 Cd10處理的531種代謝物（其中185種上調，346種下調），造成差異的主要原因是其他氨基酸的代謝、乙醛酸鹽和二羧酸鹽代謝。與CK相比，Cd10處理中所有氨基酸，包括丙氨酸、天冬氨酸、苯丙氨酸、酪氨酸、色氨酸、精氨酸、甘氨酸、絲氨酸、蘇氨酸均上調，而泛醌和其他萜類醌下調。

2.1.3 KEGG代謝途徑分析 將鑒別代謝物注釋入KEGG數據庫，共獲得76條代謝通路。由圖3可知，在Cd10處理組的531種代謝物中，乙醛酸鹽和二羧酸鹽代謝、氨基酸代謝、氨酰基-tRNA的生物合成等通路顯著富集（Plt;0.05）。Cd1處理組的283種代謝物中，上調的主要代謝通路包括檸檬酸循環、甘油磷脂代謝等。在Cd3處理組的252種代謝物中主要上調的通路有纈氨酸、亮氨酸和異亮氨酸降解通路，賴氨酸生物合成通路；主要下調的通路有纈氨酸、亮氨酸和異亮氨酸生物合成通路，脂肪酸生物合成通路等。Cd處理組172種代謝物中上調的通路有煙酸鹽和煙酰胺代謝通路，丙酮酸代謝通路，丙氨酸、天門冬氨酸和谷氨酸代謝通路；主要下調的通路有丙酸鹽代謝通路，精氨酸生物合成通路等。

2.2 轉錄組學分析

2.2.1 轉錄組序列分析 對不同水平鎘脅迫下的川芎樣品進行測序（表1），共得到97.46 Gb Cleanreads，各個樣品的Clean reads 達到5.77 Gb，其中CK 處理組共有73 976 690 條； Cd1 處理組共有67 688 858條；Cd3處理組共有61 710 731條；Cd6處理組共有63 874 313 條；Cd10 處理組共有60 114 700條。總堿基數為97 461 585 558 nt，各樣品的GC 含量占比均不低于43.61%，其中Q30堿基所占比例不低于93.05%。

2.2.2 文庫構建及Unigene功能注釋 經過組裝后，共獲得42 069條Unigenes，其中長度在2 000 bp以上的Unigene數量最多，占32.19%。Unigene總長度為69 203 497 bp，平均長度為1 645 bp。UnigeneN50長度一般用于評估組裝結果質量，當UnigeneN50≥800 bp時，認為組裝序列的完整性較好。本研究中Unigene N50長度為2 876 bp，表明川芎樣品的組裝完整性較好，可用于后續的試驗分析。通過數據庫比對，共注釋到36 336條Unigenes（表2），其中注釋率最高的是NR 數據庫，注釋Unigenes 有34 643條，占總注釋Unigenes的95.34%。

2.2.3 差異表達基因統計 如表3所示，本研究共篩選到8 569 個差異表達基因（differentiallyexpressed genes， DEGs），其中上調基因有6 859個，下調基因有1 710 個。在CK vs Cd3中篩選到的DEGs數量最多，為5 703個，包括上調基因4 205個和下調基因1 498 個；DEGs 數量最少的是CK vsCd10組，共得到195個DEGs，包括150個上調基因和45個下調基因。

2.2.4 不同處理間差異表達基因分析 對DEGs進行分析，由圖4可知，CK vs Cd1與CK vs Cd3共有DEGs 為875 個，CK vs Cd1 與CK vs Cd6 共有DEGs 為137 個，CK vs Cd1 與CK vs Cd10 共有DEGs 為76個，而3組共有DEGs僅為22個；CK vsCd3 與CK vs Cd6 共有DEGs 為 42 個，CK vs Cd3與CK vs Cd10 共有DEGs 為 51 個，CK vs Cd6 與CK vs Cd10共有DEGs為 49個，而4組共有DEGs僅有1個，為CML19。

2.2.5 KEGG通路注釋分析 通過KEGG對川芎根莖4個不同處理組DEGs進行富集分析（圖5），有4 025個DEGs被注釋到167條KEGG代謝通路中。顯著性富集分析表明，CK vs Cd3有20 條達到顯著水平（Plt;0.05），包括碳代謝、磷酸戊糖途徑、糖酵解/糖異生、氨基酸生物合成、丙酮酸代謝、酪氨酸代謝等通路；CK vs Cd1中基因數目相對較多的代謝通路是核糖體代謝；CK vs Cd6 中相關基因數目相對較多的代謝通路有甘油脂代謝、亞麻酸代謝、植物?病原體互作等；CK vsCd10中相關基因數目相對較多的代謝通路有戊糖和葡萄糖醛酸互轉化、植物激素信號傳導、內質網中的蛋白質加工、植物?病原體互作等。

2.3 轉錄組與代謝組聯合分析

為更好地探究川芎抵抗Cd脅迫時的作用機制，將差異代謝物與差異基因進行關聯分析發現，4個不同處理組共注釋得到24條通路。CK vsCd1 中共有10 條通路，其中碳代謝通路共連接73個DEGs和DAMs，數目最多；CK vs Cd3有10條通路，碳代謝連接的最多，為144個，其次是氨基酸生物合成，有116個；CK vs Cd6中連接最多的是酪氨酸代謝，共5個；CK vs Cd10中苯丙烷類生物合成與ABC 通道轉運均連接7 個DEGs 和DAMs。碳代謝與ABC通道轉運通路僅在Cd6處理組未見。4組處理共有通路有1條，為氨基酸生物合成，主要積累丙氨酸、組氨酸、苯丙氨酸、絲氨酸等氨基酸。

3 討 論

3.1 氨基酸是川芎在Cd 脅迫下的標志性物質

代謝調節是植物應對各種脅迫的重要機制之一。本研究發現，氨基酸類在川芎根系抗Cd脅迫中發揮重要作用，它們所占的比例較大，并參與了不飽和脂肪酸、異喹啉生物堿、氨酰基-tRNA和酪氨酸的生物合成，這與在Cd脅迫時水稻幼苗中的研究結果一致[15]。紫花苜蓿葉片在響應漬水脅迫時顯著上調丙氨酸、天冬氨酸和谷氨酸代謝等[16]，而川芎在受到Cd脅迫時主要上調丙氨酸、天門冬氨酸、甘氨酸、絲氨酸、蘇氨酸和谷氨酸代謝，精氨酸、苯丙氨酸、纈氨酸、酪氨酸和色氨酸等生物合成，說明植物受到非生物脅迫時，可能產生部分相同氨基酸來應答逆境脅迫。特定氨基酸及其代謝可提高植物受到逆境脅迫時的耐受性[17]，通過轉錄組與代謝組聯合分析發現，其共有通路為氨基酸生物合成，主要積累丙氨酸、組氨酸、苯丙氨酸、絲氨酸等氨基酸。組氨酸可以調節其他氨基酸的生物合成以及金屬離子的螯合和轉運[18]；馬洪娜等[19]認為，苯丙氨酸不僅提供了抗氧化劑還為朱砂根葉片適應鈣脅迫提供能量；絲氨酸是合成蛋白質和其他生物分子所必需的，包括核苷酸和絲氨酸衍生的脂類（如磷脂酰絲氨酸和鞘脂），并參與抵抗各種外界刺激[20-22]。其中，絲氨酸在輕度Cd 脅迫（1、3 mg·kg-1）中上調更大；丙氨酸在中（6 mg·kg-1）、重度（10 mg·kg-1）Cd脅迫時的上調貢獻較大，說明在不同的Cd脅迫水平下，川芎根莖中氨基酸的抗Cd脅迫模式不同。因此，氨基酸類物質可能是用于鑒別川芎根莖是否受到不同程度Cd脅迫的標志性物質。

3.2 Cd 脅迫對川芎糖類代謝和脂質代謝的影響

本研究顯示，川芎在受到Cd脅迫后的差異基因主要參與了糖類、氨基酸、脂類生物合成等代謝途徑，這與玉米、菜心等植物受到Cd脅迫后的轉錄組分析結果相似[23-24]。在植物受到非生物脅迫時，將累積糖類物質以適應環境的變化[25]。糖類的積累可以調節細胞滲透壓，以維持細胞的滲透平衡[26]。在受到干旱脅迫時，植物通過運輸蔗糖至根尖部位來抵御脅迫，維持正常生長[27]；Cd脅迫下，冬小麥根中水蘇糖、蔗糖和甘露糖含量增加[28]；Cd脅迫下的水稻中糖類含量也表現出相上升趨勢[29]。而本研究中發現，當川芎受到3 mg·kg-1 Cd脅迫時，根莖中的蔗糖含量呈現出顯著增加的趨勢，說明川芎通過積累糖類來調節細胞滲透壓，維持細胞的滲透平衡并減少Cd對自身的危害。

脂質代謝作為植物生長發育過程中的基本代謝之一，是重要的儲能物質和信號分子[30]。在受到低溫脅迫時，小麥根尖中的不飽和脂肪酸含量增加，有利于提高細胞膜的流動性，以適應外界的低溫環境[31]。本研究中，Cd脅迫下川芎可能通過增加脂肪酸、亞麻酸、甘油酯等代謝物的含量來增加細胞膜的流動性，以適應外界的Cd 脅迫環境。

3.3 川芎抗Cd 脅迫的關鍵基因CML19

對差異基因進行分析發現，CML19 在每個Cd處理組中表達都上調，推測CML19 可能為川芎抗Cd 的關鍵基因。研究表明，類鈣調蛋白（calmodulin-like protein，CML）是植物細胞中的鈣離子結合蛋白，在植物生長發育和逆境脅迫的信號轉導中具有重要功能[32]。CML9 突變體能提高擬南芥對干旱和鹽害的耐受性，并產生支鏈氨基酸的積累[33-34]；CML13 對逆境脅迫的響應非常靈敏，其在無苞芥響應逆境生理生態環境中可能起著重要作用[35]；在受到鹽脅迫時，AtCML37、AtCML38 和AtCML39 基因均上調表達[36]。

綜上所述，川芎通過調節糖類代謝、脂質代謝來響應Cd脅迫，以減少對自身的傷害。通過代謝組和轉錄組分析發現，氨基酸類可作為川芎受Cd脅迫時的標志性物質，而丙氨酸、絲氨酸等氨基酸則可用于判斷川芎受Cd脅迫程度大小。川芎受Cd脅迫時的差異基因CML19 與已有研究中有關作物調節Cd脅迫的關鍵基因相似，可能是川芎響應Cd脅迫的關鍵基因。在未來，可以對CML19進行進一步的克隆與功能鑒定，將有助于進一步闡明川芎對Cd脅迫的響應以及產生差異的分子機制。