999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

茶樹品種白毫早葉綠體基因組結構特征及其密碼子偏好性分析

2025-08-31 00:00:00曾文娟周偉何詩恬賀寧龔意輝陳致印
江蘇農業學報 2025年7期

摘要:本研究首次解析了茶樹品種白毫早(Camellia sinensis cv.Baihaozao)的完整葉綠體基因組特征。采用Ilumi-na高通量測序技術對白毫早的葉綠體基因組進行從頭組裝,并借助Geneious平臺進行基因注釋及微衛星標記檢測。同時基于IRscope工具實現基因組結構的可視化分析。在密碼子使用偏好性評估中,綜合采用相對同義密碼子使用度1 RSCU )、有效密碼子數(ENC)及中性繪圖法系統解析密碼子選擇模式及其演化驅動力。結果表明,白毫早茶樹葉綠體基因組是總長度為157025 bp 的典型四分體環狀結構,含有大小為86 586 bp 的大單拷貝區(LSC)、18 277 bp 的小單拷貝區(SSC)及1對大小各為 26081bp 的反向重復序列(IR)區域;全基因組的 含量為 37.30% ,其中IR區域的C ∴c 含量最高 42.95% ),SSC區域的 含量最低( 30.55% )。白毫早的葉綠體環狀基因組共注釋了133個功能基因,包括87個蛋白質編碼基因、37個轉運核糖核酸(tRNA)基因、8個核糖體核糖核酸(rRNA)基因及1個假基因,基因中內含子與外顯子的分布呈現顯著差異,如 trnK-UUU 的內含子最大,為2488bp,ycf2的外顯子最長,為6 897bp 。此外,本研究檢測到247個簡單重復序列(SSR)位點,其中單核苷酸重復占比較高(占比為 63.56% ),且表現出明顯的 A/T 偏好性 (97.45% ;本研究鑒定出40個長重復序列,包含20個正向重復序列、20個回文重復序列。密碼子偏好性分析結果顯示,密碼子第3位堿基中 的平均含量( GC3 值 27.59% )顯著低于密碼子第1位堿基中 G+C 的平均含量( GC1 值) 46.85% )、密碼子第2位堿基中 的平均含量 (GC2)(39.50%), 。ENC的均值為44.57,表明密碼子的偏好性較弱。中性回歸分析結果顯示,校正決定系數 (Radj2)=0.0160 ,自然選擇是起主導作用的進化動力(貢獻率為 91.47% )。 90% 最優密碼子以 A/T 結尾,這一特征與通過奇偶規則2(PR2)分析得出的密碼子第3位的G含量 (G3)gt; 密碼子第3位堿基中的C含量 (C3 )、密碼子第3位的T含量 (T3) gt;密碼子第3位的 A 含量 (A3) 的偏好性結果一致。基于31個保守蛋白質編碼基因構建系統發育樹,結果顯示,白毫早茶樹與福鼎白毫茶樹形成高支持的分支[自舉值 (bootstrap)=100% 1,屬于山茶屬核心類群。本研究結果為茶樹遺傳資源評價及葉綠體基因組進化機制的研究提供了參考。

中圖分類號:S571.1 文獻標識碼:A 文章編號: 1000-4440(2025)07-1398-14

Structural characteristics and codon usage bias analysis of the chloroplast genome in the tea cultivar Camellia sinensis cv. Baihaozao

ZENG Wenjuan1,2.3.4, ZHOU Wei1,2,3.4, HE Shitian1,2,3,4

HE Ning1,2.3.4, GONG Yihui1,2,.3.4, ,CHEN Zhiyin1,2,.3,4 (1.College ofAgriculture and Biotechnology,Hunan Universityof Humanities,Scienceand Technology,Loudi4170oo,China; 2.Innovationand Entrepreneurship Education CenterforHorticultural ProductionandProcessinginHunanProvince,HunanUniversityofHumanities,Science and Technology,Loudi 417ooo,China;3.Hunan University ofHumanities,Scienceand Technology,HunanProvincial Innovationand EntrepreneurshipDemonstrationBase,Loudi417OoO;4.KeyLaboratoryof Characteristic Agricultural ResourceDevelopmentand Quality Safety Con

Abstract:This study reports the first complete chloroplast genome of the tea cultivar Camelia sinensis cv.Baihaozao. Using Ilumina high-throughputsequencing,wedenovoassembledthechloroplastgenomeand performedgeneannotationandmicrosatellitedetectioninGeneious.Aditionally,genomestructurevisualizationanalysisasachievedwiththeIRscopetool.Inthe assessment of codon usage bias,we comprehensively employed the relative synonymous codon usage ( RSCU ),the effective numberofcodons(ENC),andneutralityplotanalysistosystematicallanalyzecodonselectionpaterns andtheirevolutionarydriving forces.Theresultsrevealedthatthechloroplast genomeofCameliasinensis cv.Bahaozao exhibitedatypicalquadripartitecircularstructure withatotallengthof157O25bp,comprisingalargesingle-copy(LSC)regionof 86586bp,asmallsingle-copy (SSC)region of 18 277 bp,and a pair of inverted repeat(IR)regions each 26081 bp in length. The overall G+C content of the chloroplast genome was 37.30% ,with the IR regions exhibiting the highest value (42.95 % )and the SSC region showing the lowest ( 30.55% ).The circular chloroplast genome wasannotated with133 functional genes,comprising 87 protein-coding genes,37 transfer RNA(tRNA)genes,eight ribosomal RNA(rRNA)genes,and onepseudogene. Gene structure analysis revealed marked variationin intronand exon distribution. The largest intron 2488bp )was identified in trnK- UUU ,while the longest exon(6897bp)wasfoundinyef2.Furthermore,247simplesequencerepeat (SSR)loci were detected,with mononucleotiderepeats constituting the majority( 63.56% )and showing a pronounced A/T bias( 97.45% ).We identified 40 long repeat sequences,comprising 2O forwardrepeatsand 20palindromic repeats.Codon usagebiasanalysis showedthat the average G+Ccontent at the third codon position ( GC3 ,27. 59 % )was significantly lower than that at the first ( GC1 ,46.85 % ) and second ( GC2 , 39.50% )positions.The meanvalueof ENC was44.57,indicating that thecodonbiaswasrelativelyweak.Theresultsof theneutrality regression analysis showed that the adjusted coefficient of determination( Radj2 )value was 0.0160 ,indicating that natural selection was the dominant evolutionary force (with a contribution rate of 91.47 % ).Ninety percent of optimal codons ended with A/T ,consistent withthe bias patern identified byparityrule2(PR2)analysis:G contentat the third position (G3)gt;C content at the third position ( C3 ),and T content at the third position (T3)gt;A content at the third position( A3 ).Phylogenetic analysisbasedon31conservedprotein-codinggenesdemonstratedthatCameliasinensiscv.Baihaozaoformedahighlysupported clade(bootstrap value =100% )with C .sinensiscv.FudingBaihao,belonging to thecore Camellia clade.This studyprovides references for evaluating tea genetic resources and investigating the evolutionary mechanisms of chloroplast genomes.

Key words: Camellia sinensis;chloroplast genome;codon usage bias;optimal codons

茶科(Theaceae)植物是重要的經濟作物,其中茶樹(Camelliasinensis)因其所產茶葉的飲用價值及茶籽的工業用途,在全球范圍內具有較好的生態和經濟意義[1]。中國作為茶樹的起源中心,擁有較高的茶樹種質資源多樣性,但是目前關于茶樹的分類、起源與進化機制仍存在諸多爭議[2]。近年來,葉綠體基因組因其結構穩定、進化速率適中且單親遺傳的特性,成為解析植物系統發育、基因表達調控及適應性進化的重要媒介[1-2]

葉綠體基因組的密碼子使用偏好性(Codon usagebias,CUB)是遺傳信息傳遞中的關鍵特征,反映了同義密碼子的非隨機使用模式[3]。研究發現,高等植物葉綠體基因組普遍存在第3密碼子位點的腺嘌呤(A)/胸腺嘧啶(T)使用偏好性,這可能與突變壓力和自然選擇的共同作用有關[1,4]。例如,在茶科植物中,絕大多數(約 96.55% )高頻使用的密碼子(即同義密碼子中明顯更受偏愛的類型)以堿基A或T結尾。這種密碼子的偏好模式越強,基因的表達水平就越高[1]。然而,目前對特定茶樹品種(如茶樹品種白毫早)葉綠體基因組CUB的研究仍較匱乏,其基因組結構特征與密碼子偏好模式的關聯尚未明晰[2.5]

高通量測序分析結果表明,不同茶樹品種的葉綠體基因組具有高度保守的四分體環狀結構,包括大單拷貝區(LSC)、小單拷貝區(SSC)、反向重復區( IRa/IRb )及130~146個功能基因[含蛋白質編碼基因、轉運核糖核酸(tRNA)基因及核糖體核糖核酸(rRNA)基因][6。基因組中同時存在高變異區域(如psbA、matK、petN-psbM、trnK-rps16等),其序列差異度 gt;1.5% ,可作為系統發育分子標記用于物種分類與遺傳分化研究[4]。然而,目前關于重要的茶樹栽培品種白毫早葉綠體基因組的組裝、注釋及密碼子偏好性的分析仍未見系統報道。此外,現有研究多聚焦于物種間的比較,關于品種內CUB的驅動因素(如自然選擇與突變壓力的相對貢獻)及其對基因異源表達的指導意義仍需深人探討[7-8]

本研究擬以茶樹品種白毫早為研究對象,結合Illumina測序與生物信息學工具,系統解析其葉綠體基因組的結構特征、高變異區域及密碼子使用模式[7.9]。通過中性繪圖、有效密碼子數(ENC)計算、所有密碼子第3位堿基中 G+C 的平均含量占第3位堿基總含量的比例( GC3s 值)分析及奇偶規則2(PR2)偏倚性檢驗,明確CUB的主要驅動機制。研究結果以期能夠豐富茶樹的葉綠體基因組數據庫,為茶樹的分子育種、系統進化分析及功能基因挖掘奠定基礎

1 材料與方法

1.1 試驗材料

茶樹品種白毫早(Camelliasinensiscv.Baiha-ozao)的幼嫩葉片樣本于2022年6月采自省市新化縣川坳村桃花源農業開發有限公司的標準化種植園( 111.28°E,27.73°N) 。樣本經液氮速凍后用干冰保存,并轉運至南京集思慧遠生物科技有限公司進行葉綠體基因組測序。本研究所獲得的葉綠體基因組序列已提交至美國國家生物技術信息中心(NCBI)數據庫(登錄號:PQ066596),數據開放獲取。

1.2試驗方法

1.2.1DNA提取與測序采用改良十六烷基三甲基溴化銨(CTAB)法提取茶樹葉片總脫氧核糖核酸(DNA)[10-11],經純度及完整性檢測合格后,用超聲波機械打斷法進行DNA的片段化處理。片段化產物經末端修復 ?3 端加A、測序接頭連接及瓊脂糖凝膠電泳篩選后,構建Illumina測序文庫。文庫質檢合格后,用Illumina NovaSeq60OO平臺進行雙端測序。原始數據用 fastp(v0.20.0) 進行質量控制,具體步驟如下:(1)剔除含接頭及引物序列的短序列數據片段(Reads);(2)移除整體堿基平均錯誤率高于 31.6% 的低質量測序數據;(3)去除N端堿基數超過5個的低置信度Reads[10]。獲得的高質量測序數據用于后續組裝與分析。

1.2.2葉綠體基因組組裝為了降低序列的復雜性,首先用Bowtie2 v2.2.4 的高靈敏度局部比對模式(Very-sensitive-local)對原始測序數據與自建葉綠體基因組數據庫進行比對,篩選出特異性葉綠體DNA序列[12]。然后,基于二代測序數據特征,采用無參考基因組的從頭組裝(Denovo)策略,用SPAdesv3.10.1軟件將DNA序列分割成不同 k -mer值(55、87、121)的連續短片段進行迭代組裝,獲得初始輸入(SEED)序列[13]。針對由基因組重復序列造成的組裝間隙,用SSPACE 構建支架(Scaffolds),并用GapFillerv2.1.1 進行間隙填補[14]。對于持續性缺口區域,采用引物設計軟件( Primer5.0 )結合聚合酶鏈式反應(PCR)擴增技術驗證的方法進行序列補全,確保基因組的完整性。最后,用Bowtie2對組裝序列進行二次比對校正,并基于葉綠體典型的環狀四分體結構[大單拷貝區(LSC)、小單拷貝區(SSC)及反向重復區1 ΔIRa/IRb )]對基因組坐標進行標準化重排,以消除因反向重復區邊界滑動導致的基因位置偏移[15]

1.2.3葉綠體基因組的注釋本研究基于Geneious平臺,采用雙策略整合分析流程對葉綠體基因組進行系統注釋。首先,運用prodigal v2.6.3 預測蛋白質的編碼區,通過hmmerv3.1b2鑒定核糖體RNA(rRNA)編碼區,并用aragornv1.2.38預測轉運RNA(tRNA)編碼區[12]。同時,基于NCBI數據庫獲取近緣物種的葉綠體基因序列,用BLAST v2.6 進行同源比對,獲取保守基因的注釋結果[12]。然后,通過人工校驗消除注釋差異,重點修正起始密碼子/終止密碼子定位誤差、反密碼子錯誤及基因冗余問題[6],并用tRNAscan-SE1.23對tRNA進行二次驗證[11,17]。針對反向重復序列(IR)區域,采用雙重同源性閾值法驗證其邊界完整性,確保發生轉剪接的基因座(或產生轉剪接轉錄本的基因)及rRNA基因的定位正確[16]。最后,整合多證據注釋結果,精確定義多外顯子基因的剪接位點,形成標準化基因組注釋文件[12,16],并借助IRscope工具實現基因組結構的可視化分析。

1.2.4葉綠體基因組的特征分析散在重復序列(Dispersedrepeats)用vmatch v2.3.0 軟件結合開源動態腳本編輯語言(Perl)進行系統鑒定,參數設置如下:最小重復單元長度為 30bp ,海明距離閾值為3,涵蓋正向、回文、反向及互補4種重復類型[16-17]。簡單重復序列(SSR)的分析用MISAv1.0工具完成,長重復序列的分析借助REPuter,檢測標準為單核苷酸重復次數 ?8 次,二核苷酸重復次數 ?5 次,三核苷酸至六核苷酸重復次數均 ?3 次,以明確葉綠體簡單重復序列的分布特征[11,16]。基于近緣物種葉綠體基因組數據,用CGVIEW軟件進行全基因組結構的共線性比較,以揭示基因組重排及重復序列的演化規律[17] O

1.2.5 密碼子偏好性分析

1.2.5.1中性分析使用密碼子第3位堿基中 G+C 的平均含量( GC3 值)與密碼子第1、第2位堿基中G+C 的平均含量( GC12 值)之間的線性回歸模型,評估突變壓力和自然選擇在密碼子使用偏好中的相對貢獻。如果回歸斜率趨近1且決定系數 (R2)gt;0.5 表明密碼子偏好性由突變壓力主導;反之則提示自然選擇占主導地位[18]

1.2.5.2 ENC- GC3 分析基于ENC值與 GC3 值的分布關系,通過理論曲線 ENC=2+GC3+29/[GC32+ (1-GC32] 評估核苷酸組成對密碼子使用的限制效應。若數據點位置顯著低于理論曲線,表明存在翻譯優化或自然選擇壓力[19]

1.2.5.3PR2分析以密碼子第3位堿基中G含量占G+C 含量的比例 [G3/(G3+C3)] 和密碼子第3位堿基中A含量占 含量的比例 [A3/(A3+T3) ]的偏離程度為指標,基于二等分規則(PR2)的偏差圖,評估突變壓力與選擇壓力在功能位點間的非對稱程度。中心點(0.5,0.5)用于表征無偏狀態,顯著偏離中心點則反映選擇或突變壓力具有非對稱作用[18,20]

1.2.5.4最優密碼子的篩選依據ENC值降序排列,取前十分位組(ENC值 90% 分位數)作為高表達組,后十分位組(ENC值 ?10% 分位數)作為低表達組,用CodonW計算相對同義密碼子使用度的變化(△RSCU), ΔRSCU=RSCU-RSCU? 。最優密碼子的篩選標準為 ΔRSCUgtrsim0.08 且 RSCUefgt;1.00 ,以確定受翻譯優化驅動的核心最優密碼子[21-22]

1.2.6葉綠體基因組的系統發育分析首先對環形葉綠體基因組進行同源起點校正,用MAFFT軟件(v7.427) 的自動優化模式進行多序列比對,該算法利用動態選擇比對策略提升序列的匹配精度[9]。隨后基于廣義時間可逆Gamma模型(GTRGAMMA核昔酸替代模型),用RAxMLv8.2.10軟件進行最大似然法分析,設置快速自舉檢驗(Rapidbootstrap)重復次數為1000次,以評估節點支持率,構建系統發育樹[6.23]

2 結果與分析

2.1茶樹品種白毫早葉綠體基因組的結構特征

本研究基于IlluminaNovaSeq 6O0O高通量測序平臺對白毫早茶樹(Camellia sinensiscv.Baihaozao)的葉綠體基因組進行深度測序。結果顯示,測序數據質量優異,質量數大于20的堿基所占百分比( )和質量數大于30的堿基所占百分比( Q30) 分別達到 97.00% 92.27% (對應堿基錯誤率分別低于1.0%.0.1% ),累計獲得大小為16140581400bp的高質量序列,為后續基因組組裝提供了基礎[24]。

白毫早茶樹的葉綠體基因組呈現典型的四分體環狀結構,總長度為 157 025bp ,由大小為86586bp的大單拷貝區(LSC)、大小為18277bp的小單拷貝區(SSC)和1對大小各為26081bp的反向重復區( IRa/IRb )組成。其全基因組的 G+C 含量為37.30% ,與菊科植物藜蒿(Artemisia selengensis)物種具有較高的相似度(相似度為 37.5% )[25]。區域 G?C 含量分析結果顯示,IR區富含核糖體RNA基因, 含量最高( 42.95% ),顯著高于LSC區( 35.33% )、SSC區( 30.55% ),該特征與木姜子屬(Litsea)植物(IR區的 G+C 含量為 42.68% )、腺果藤科(Biebersteiniace-ae)植物(IR區的 G+C 含量為 42.68% )的葉綠體基因組規律一致[26-27]。值得注意的是,SSC 區的 G+C 含量( 30.55% )接近清風藤屬( Sabia) 物種的 G+C 含量 32.00% ),表明該區域在進化過程中可能承受了更強的 A/T 堿基選擇壓力[28]

2.2茶樹品種白毫早葉綠體基因組的基因組成

白毫早茶樹葉綠體基因組的基因組成分析結果顯示,其環狀基因組共注釋了133個功能基因,包括87個蛋白質編碼基因、37個轉運核糖核酸(tRNA)基因、8個核糖體核糖核酸(rRNA)基因及1個假基因(表1),基因總數及功能分類(包括蛋白質編碼基因、tRNA基因、rRNA基因和假基因)在系統發育近緣物種間高度保守,體現了葉綠體基因組的保守結構特征[29-31]。由表1可知,光系統I亞基與光系統Ⅱ亞基構成光能捕獲與轉化的核心單元。其中光系統I亞基編碼基因包含psaA、psaB等5個,參與光能吸收與電子傳遞;光系統I亞基編碼基因則包含15個,負責水的光解與氧氣釋放。值得注意的是,光系統I亞基編碼基因數量明顯多于光系統I亞基編碼基因,暗示其在光反應中承擔更復雜的調控功能。在能量代謝相關基因中,NADH脫氫酶亞基編碼基因(ndh基因家族)包含12個成員,其中ndhB有2個拷貝,ndhA、ndhB各有1個內含子,提示其可能通過基因重復與可變剪接實現功能的多樣性,該基因可能在電子傳遞鏈中起到關鍵作用。細胞色素 Δb/f 復合體編碼基因(pet基因家族)與ATP合酶亞基編碼基因( atp 基因家族)各包含6個,共同構成質子梯度驅動的ATP合成系統。

白毫早茶樹葉綠體 clpP,ycf3 基因各含有2個內含子, atpF,ndhA,ndhB 等基因含有單內含子。白毫早茶樹葉綠體基因組的內含子長度差異顯著,其中 trnK–UUU 的內含子最大,大小為2488bp,trmL-UAA的內含子最小,大小為 ;外顯子長度差異亦明顯,其中 trnG-GCC 的外顯子長度僅為 60bp ,而ycf2的外顯子長度達 6897bp ,提示不同基因在進化過程中受到的選擇壓力不同[32-33]

此外,LSC區的基因密度(143個基因中含有35個外顯子、19個內含子)顯著高于SSC區和IR區,與多數被子植物葉綠體基因組的區域功能分工模式相符[34]。假基因的存在及部分基因內含子數量的變異可能反映基因組重排或自然選擇驅動的功能存在冗余[31]。本結果可為后續密碼子偏好性分析及系統發育研究提供結構基礎,亦可為茶樹屬葉綠體基因組演化研究補充數據[6]

表1茶樹品種白毫早葉綠體基因組的組成

Table1 Gene composition of thechloroplast genome of Camellia sinensis cv.Baihaozao

基因名稱后面標有*表示該基因有1個內含子,標有**表示該基因有2個內含子,標有(2)表示該基因有2個拷貝,標有#表示該基因為假基因。NADH表示還原型煙酰胺腺嘌呤二核苷酸, ATP 表示三磷酸腺苷,RNA表示核糖核酸。

2.3簡單重復序列(SSR)和長重復序列

本研究對白毫早茶樹葉綠體基因組中的簡單重復序列(SSR)和長重復序列進行了系統分析。由表2可以看出,從葉綠體基因組中共檢測到247個SSR位點,其中單核苷酸重復( 63.56% )占絕對優勢,顯著高于二核苷酸重復( 1.62% )、三核苷酸重復( 29.15% )、四核苷酸重復( 4.86% )及五核苷酸重復 (0.81%) 。值得注意的是,單核苷酸重復中A/T型占比高達 97.45% ,這種 A/T 偏倚性與葉綠體基因組整體表現出的較高的A、T含量特征[35-36]一致,反映了植物葉綠體SSR的進化保守性[37]

由表3可以看出,通過對長重復序列分析,共鑒定出40個長重復序列,包括20個正向重復序列和20個回文重復序列,未檢測到反向重復序列或互補重復序列。其中,正向重復序列的長度為 30~ 82bp 。這一分布模式與部分物種中同時存在正向、反向及回文重復序列的研究結果[38]形成對比,提示白毫早茶樹可能存在葉綠體基因組長重復序列類型的特異性。此外,回文重復的序列長度跨度較大,范圍為 30~26081bp ,可能與其在基因組結構重排或調控元件中的功能相關[37]

表2茶樹品種白毫早葉綠體基因組中簡單重復序列(SSR)的類型及數量分布

Table2Typeandquantitydistributionofsimplsequenceepeats(SSRs)inthechoroplast genomeofCamellasinensiscvBaoao

表3茶樹品種白毫早葉綠體基因組中長重復序列

able 3Long repeat sequences in the chloroplast genome of Camelia sinensis cv.Baihaozai

續表3 Continued3

P:回文重復;F:正向重復; E 值:評估2個序列比對結果統計學顯著性的指標,當 E 值小于 1.0×10-10Hf ,認為比對結果具有顯著性意義。

與已報道的植物葉綠體SSR特征相比,本研究發現的SSR總數(247個)顯著高于豆科胡枝子屬(Lespedeza)植物(75~78個單核苷酸重復)[1],但低于甘蔗(2199個單核苷酸重復)[39],表明SSR的豐度存在顯著的物種特異性。值得注意的是,四核苷酸重復(12個)和五核昔酸重復(2個)的檢出,補充了部分研究中“葉綠體SSR以單核苷酸/雙核苷酸/三核苷酸為主”的結果[40-41] O

2.4茶樹品種白毫早葉綠體基因組密碼子使用偏好性

2.4.1同義密碼子的偏好性由圖1可以看出,白毫早茶樹葉綠體基因組中52個蛋白質編碼基因的GC3 均值為 27.59% ,顯著低于 GC1 的均值( 46.85% ))和 GC2 的均值( 39.50% ),呈現 GC1gt;GC2gt;GC3 的規律。該模式與多種被子植物[如向日葵(HelianthusannuusJ-01)[42]]及單子葉植物葉綠體基因組的堿基組成特征[43-44]一致,反映了第3位密碼子對 A/T的顯著偏好,可能與葉綠體基因組中普遍存在的 A T偏向性突變壓力相關。進一步分析發現,基因組的ENC值為 34.24~53.16 ,均值為44.57,其中51個基因的ENC值 gt;35.00 ,表明整體密碼子的偏好性較弱,提示同義密碼子的使用主要受到突變-漂變平衡影響而非受到強選擇性壓力影響[42.45]。值得注意的是,僅 基因的ENC值(34.24)顯著低于其他基因,可能由于其功能需求(如翻譯效率或mRNA穩定性)受到局部自然選擇的影響。

2.4.2中性分析( GC3–GC12 分析)如表4所示,白毫早茶樹葉綠體基因組 GC3 值 19.85%~37.02%) 整體低于 GC12 值 32.11%~55.04% ,該差異與茶科植物葉綠體基因組中普遍觀察到的 GC3 偏向性一致,回歸分析結果顯示其差異主要源于弱選擇約束下的 A/T 偏向突變壓力[1,46]。通過構建 GC12(Y) 與GC3(x) 的中性回歸模型,發現二者的線性回歸方程為 Y= 0.408 0+ 0.085 3x ,對應的校正決定系數(Radj2)=0.0160 ,回歸斜率為0.0853,表明自然選擇對密碼子偏好形成的貢獻率為 91.47% ,顯著高于突變壓力的貢獻( 8.53% )[4748]。低 Radj2 ( Radj2lt; 0.0500)提示 GC12 與 GC3 間無顯著相關性,進一步支持葉綠體基因組密碼子的偏好性主要受到自然選擇驅動而非受到突變偏置主導的假設[44,47]。該研究結果反映了葉綠體基因在進化過程中受翻譯效率優化等選擇壓力約束的保守性[1,49]。

基因

a1:acCD;b1:atpA;c1:atpB;d:atpE;e1:atpF;f1:atpI;g1:ccsA;h1:emA;i:clpP;j1 :matK;k1:ndhA;1 :ndhB; m1 : ndhC;n1 : ndhE;01 : ndhF ;P1 :ndhG;q1 : ndhH; r1 : ndhl;s1 :ndl ιJ;t1:ndhK;u1:petA;v1:petB;w1:petD;x1:psaA;y1:psaB;z1:psbA;a2:psaA. L;e2:rpl14;f allostsg2:rpl2;h2:rpl20;i2:γpl22;i2:γpl2;i2:γpoli;i2:γpoli1;r2:γpol2;i2:γpol2;i2:γpol2;i2:γpsl2;i2:γpsl2;γpsl2;γpsl2;γpsl2;γpsl2;γps2:γps32:γps32:γpsl42:γpsl42:γpsl5. (20(20 rps7;v2:rps8;w2:ycfl;x2:ycf2;y2:ycf3;z2:ycf4amp;GC1 :密碼子第1位堿基中 G+C 的平均含量; GC2 :密碼子第2位堿基中 G+C 的平均含量;GC3 :密碼子第3位堿基中 G+C 的平均含量; GCall :密碼子所有堿基位中 σG+C 的平均含量; ENC :有效密碼子數; GC3s :所有密碼子第3位堿基中 含量占第3位堿基總含量的比例。

圖1 G+C 含量與有效密碼子數( ENC 值)在茶樹品種白毫早葉綠體基因間的對比分析結果

ig.1Comparative analysis results of G+C content and effective number of codons(ENC values)among genes in the chloroplast genome ofCamellia sinensiscv.Baihaozao

2.4.3 GC3 -ENC分析由表4還可以看出,52個基因的ENC值多數分布在 39.00~54.00 ,平均值為44.57,其中29個基因的ENC值 gt;45.00 ,表明白毫早茶樹葉綠體基因組的整體密碼子偏好性較弱[50]值得注意的是,所有基因的實測ENC值均極顯著低于基于 GC3 組成的理論預期值[假設當密碼子的使用僅由突變壓力(即 G+C 含量)驅動,不受自然選擇影響時,通過數學模型計算出的ENC的預測值]( Plt;0.01? ,提示密碼子使用模式受到非突變壓力因素的強烈影響[51-52]。 GC3 值分布范圍狹窄( 19.85%~37.02% ),這一現象與茶科植物(Theace-ae)及向日葵(Helianthusannuus)葉綠體基因組的GC3 分布模式相似,均呈現對 A/T 的偏好性[1]

表4茶樹品種白毫早葉綠體基因組密碼子分析結果

Table 4Analysis of codons in the chloroplast genome of Camelia sinensis cv.Baihaozao

NADH:還原型煙酰胺腺嘌呤二核苷酸;ATP:三磷酸腺苷;ENC:有效密碼子數; GC12 :密碼子第1、2位堿基中 的平均含量; GC3 :密碼子第3位堿基中 G+C 的平均含量; G3/(G3+C3) :密碼子第3位堿基中G含量占 G+C 總含量的比例 ;A3/(A3+T3) :密碼子第3位堿基中A含量占 A+T 總含量的比例。

2.4.4PR2分析根據PR2偏倚分析結果,白毫早茶樹葉綠體基因組中密碼子第3位堿基使用偏好性呈現顯著特異性。通過計算 G3/(G3+C3) 和 A3/ ( ∣A3+T3∣ )發現,基因密碼子第3位普遍存在 G3gt;C3 [G3/(G3+C3) 均值 gt;0.50] 和 T3gt;A3[A3/(A3+T3) 均值 lt;0.50] 的特點(表4)。值得注意的是, psbA 基因的 G3/(G3+C3) 值和 A3/(A3+T3) 值均最低(0.33),表明其密碼子第3位對 C,T 的偏好性最強,這與葉綠體 psbA 基因密碼子普遍存在以C結尾的偏好性一致[53]。然而, rpsI4 基因的 G3/(G3+C3) 值最高(0.72),rps3基因的 A3/(A3+T3 )值最高(0.60),提示核糖體蛋白質編碼基因可能受到不同于其他功能基因的選擇性壓力[27,54]

2.4.5最優密碼子的確定根據篩選標準( RSCUgt; 1.00且 ),本研究一共從白毫早茶樹葉綠體基因組中鑒定出20個最優密碼子,其中 90% 的密碼子以A或T結尾(如GCA、GCT、CGA等),比例顯著高于以 G/C 結尾的密碼子(表5)。這一結果表明,白毫早茶樹葉綠體基因組表現出對A/T結尾密碼子的強烈偏好,與蕨類植物鐵線蕨(Adiantumcapillus-veneris)( 75% 以 A/T 結尾)、理查德水蕨(Ceratopterisrichardii)( 89% 以A/T結尾)等物種的葉綠體密碼子使用模式高度一致[51]。進一步分析發現,這種偏好性可能與葉綠體基因組固有的高 A T含量和突變偏置密切相關[44,49],同時自然選擇壓力可能通過優化翻譯效率強化了該趨勢[44,55]。值得注意的是,本研究篩選的最優密碼子中,GCA(編碼丙氨酸)、GCT(編碼丙氨酸)等高頻密碼子與柴胡(Bupleurumfalcatum)葉綠體基因組中高表達密碼子存在重疊[56],提示不同物種間在葉綠體基因表達調控機制方面存在一定的保守性。

表5茶樹品種白毫早葉綠體基因的最優密碼子

Table5Optimal codonsinthechloroplast genomeofCamellia sinensiscv.Baihaoza(

2.5 系統發育分析

為了明確白毫早茶樹的進化地位,本研究選取禾本科的玉蜀黍屬(Zea)、甘蔗屬(Saccharum)、小麥屬(Triticum)及稻屬(Oryza)植物作為外群,基于31個保守蛋白質編碼基因的聯合比對數據矩陣,采用最大似然法(ML)構建系統發育樹。圖2結果顯示,白毫早茶樹與福鼎白毫茶樹(C.sinensiscv.Fu-dingBaihao)形成高度支持的分支[自舉值(boot-strap)=100% 1,并與山茶屬中的茶(C.sinensis,OL450428.1)聚為一類,證實其分類地位屬于山茶屬的核心類群[57-58]。值得注意的是,白毫早茶樹所在的山茶屬(Camellia)與多瓣茶屬(Polyspora)構成姐妹群,而紫莖屬(Stewartia)則位于較遠分支,表明白毫早茶樹與多瓣茶屬植物的遺傳分化程度低于其與紫莖屬植物的分化程度。這一結果與葉綠體基因組的結構保守性特征一致,進一步支持葉綠體蛋白質編碼基因在解析山茶屬近緣物種關系方面的有效性[57,59]

3討論

本研究揭示了白毫早茶樹葉綠體基因組的典型四分體結構(LSC/SSC/IRa/IRb)及其獨特的組成特征,其總長度為 157025bp ,與多數被子植物葉綠體基因組規模相當[30.32.60],而IR區的 G+C 含量為42.95% ,顯著高于LSC區( 35.33% )和SSC區中 30.55% )。值得注意的是,基因注釋結果顯示, trnK- UUU的內含子大小為 2488bp ,遠超多數植物葉綠體基因的內含子長度[如金腰屬(Chrysosplenium)植物的最大內含子長度僅為 523bp[61] ],可能與其特殊的RNA編輯機制相關。

在密碼子偏好性方面, GC3 的均值( 27.59% )顯著低于 GC1(46.85%) 和 GC2(39.50%) ,呈現 GC1gt; GC2gt;GC3 的遞減規律,與雙子葉植物線粒體基因組中觀察到的嘧啶偏好性( T3/C3 優勢)[62]一致。ENC值分布(均值為44.57)及 GC12 與 GC3 呈低相關性(Radj2=0.016 0) ,表明自然選擇的貢獻率顯著高于突變壓力。值得注意的是,最優密碼子中 90% 以 A/Λ T結尾,與普通小麥葉綠體基因組中以 ΔA/T 結尾的密碼子占優勢形成呼應[2],提示翻譯效率優化可能驅動密碼子選擇,這與茶樹CsNRT1. 1 基因中以 N T結尾的密碼子占比高達 85% 的研究結果相互印證,凸顯了山茶科植物葉綠體基因組在密碼子選擇機制上的趨同演化特征[63]

SSR分析結果顯示,在247個簡單重復中,單核苷酸重復占比為 63.56% ( A/T 占絕對優勢),與山楂屬植物葉綠體基因組的SSR類型分布(單核苷酸重復占比為 65.20% )相比更為保守,反映了二者在微衛星序列演化過程中均受到A、T堿基優勢的強烈驅動[64]

4結論

本研究基于白毫早茶樹葉綠體基因組的結構與進化特征研究,為山茶屬系統分類提供了分子證據。白毫早茶樹葉綠體基因組呈現典型四分體環狀結構(總長 157 025bp ),包括大單拷貝區(LSC)、小單拷貝區(SSC)和1對倒置重復區(IR)。 G+C 含量分布具有顯著的區域保守性:IR區( 42.95% )因富含rRNA基因而較高,SSC區( 30.55% )因A/T偏向基因富集而最低。白毫早葉綠體133個功能基因的數量與組成與山茶屬物種高度一致,印證了該屬基因組的進化保守性。SSR分析結果顯示,在單核苷酸重復中單核苷酸A/T重復占比 97.45% ,此類重復可作為高變異位點用于后續群體遺傳研究。密碼子偏好性分析結果表明,自然選擇是主要驅動因素,貢獻率為 91.47% T 90% 最優密碼子以A/T結尾,反映葉綠體基因組受到翻譯選擇壓力影響的進化特征。系統發育分析結果進一步證實,白毫早茶樹與福鼎白毫茶樹形成高度支持的分支(bootstrap Σ=Σ 100% ),支持其歸屬于山茶屬核心類群。本研究從基因組層面揭示了茶樹品種白毫早的分子標記與進化保守性,為山茶屬系統分類及葉綠體功能進化研究提供了數據支撐。

參考文獻:

[1]WANGZJ,CAIQW,WANGY,etal.Comparative analysisof codonbias in the chloroplast genomes of Theaceae species[J]. Frontiersin Genetics,2022,13:824610.

[2] 楊雨青,譚娟,汪芳,等.茶樹葉綠體基因組的研究與應 用進展[J].生物技術通報,2024,40(2):20-30.

[3] 趙洋,劉振,楊培迪,等.密碼子偏性分析方法及茶樹中密 碼子偏性研究進展[J].茶葉通訊,2016,43(2):3-7.

[4] LIW,ZHANGCP,GUO X,etal.Complete chloroplast genome ofCamellia japonica genome structures,comparative and phylogeneticanalysis[J].PLoSOne,2019,14(5):e0216645.

[5] 徐禮羿.茶樹SNP高密度遺傳連鎖圖譜構建與數量性狀候選 基因挖掘[D].武漢:華中農業大學,2019.

[6] CHENZY,LIUQ,XIAO Y,et al.Complete chloroplast genome sequence of Camellia sinensis:genome structure,adaptive evoluics,2023,64(3) :419-429.

[7]王占軍,李豹,姜行舟,等.兩種茶樹全基因組數據的密碼子 偏好性比較分析[J].中國細胞生物學學報,2018,40(12): 2028-2039.

[8]王占軍,吳子琦,王朝霞,等.3個茶樹品種WOX基因家族的進 化及密碼子偏好性比較[J].南京林業大學學報(自然科學 版),2022,46(2):71-80.

[9]KATOH K, STANDLEY D M. MAFFT multiple sequence alignment software version 7:improvements in performance and usability [J].Molecular Biology and Evolution,2013,30(4):772-780.

[10]YU YF,OUYANG Z,GUOJ,et al. Complete chloroplast genome sequence of Erigeron breviscapus and characterization of chloroplast regulatory elements[J].Frontiers in Plant Science,2021,12: 758290.

[11] SOMARATNE Y,GUAN D L,WANG WQ,et al. The complete chloroplast genomes of two Lespedeza species:insights into codon usage bias,RNA editing sites,and phylogenetic relationships in desmodieae(Fabaceae:Papilionoideae)[J].Plants,2O19,9(1): 51.

[12]陸奇豐,黃至歡,駱文華.極小種群瀕危植物廣西火桐、丹霞梧 桐的葉綠體基因組特征[J].生物多樣性,2021,29(5):586- 595.

[13]GRUENSTAEUDL M,GERSCHLERN,BORSCH T.Bioinformatic workflows for generating complete plastid genome sequences-an example from Cabomba (Cabombaceae) in the context of the phylogenomicanalysisof the water-lilyclade[J].Life,2O18,8(3): 25.

[14]吳東洋.以楊柳科為例葉綠體基因組組裝分析流程管理平臺 開發及應用[D].南京:南京林業大學,2020.

[15] ZHANG T W, ZHANG X W,HU S N,et al. An efficient procedure for plant organellar genome assembly,based on whole genome data from the 454 GS FLX sequencing platform[J]. Plant Methods,2011,7:38.

[16]QU XJ,ZOU D,ZHANG RY,et al. Progress,challenge and prospect of plant plastome annotation[J].Frontiers in Plant Science,2023,14:1166140.

[17]GICHIRA A W,LI Z Z,SAINA JK,et al. The complete chloroplast genome sequence of an endemic monotypic genus Hagenia (Rosaceae) :structural comparative analysis,gene content and microsatellite detection[J].Peer J,2017,5:e2846.

[18]PARVATHY ST,UDAYASURIYANV,BHADANA V.Codon usage bias[J].Molecular Biology Reports,2022,49(1):539-565.

[19]CAMIOLO S,MELITO S,PORCEDDU A.New insights into the interplaybetween codon bias determinants in plants[J].DNA Research,2015,22(6) :461-470.

[20]QIAO Z S,LIJQ,ZHANG XL,et al,Genome-wide identification,expression analysis,and subcellular localization of DET2 gene family in Populus yunnanensis[J]. Genes,2024,15(2):148.

[21]HERSHBERG R,PETROV D A. General rules for optimal codon cnoIce[J」. rLos Geneucs,zuuy,J(/):e1UuusJo.

[22]ANGOV E. Codon usage:nature’s roadmap to expression and folding of proteins[J].Biotechnology Journal,2011,6(6):650-659.

[23]STAMATAKIS A.RAxML version 8:a tool for phylogenetic analysisand post-analysis of large phylogenies[J].Bioinformatics, 2014,30(9):1312-1313.

[24]王興春,楊致榮,王敏,等.高通量測序技術及其應用[J].中 國生物工程雜志,2012,32(1):109-114.

[25]WANGYH,WEI QY,XUETY,etal. Comparative and phylogenetic analysis of the complete chloroplast genomes of 1O Artemisia selengensis resources based on high-throughput sequencing[J]. BMC Genomics,2024,25(1) :561.

[26]SONGWC,CHEN ZM,SHI WB,et al.Comparative analysis of complete chloroplast genomes of nine species of Litsea(Lauraceae):hypervariable regions,positive selection,and phylogenetic relationships[J].Genes,2022,13(9):1550.

[27]CHI XF, ZHANG FQ,DONG Q,et al. Insights into comparative genomics,codon usage bias,and phylogenetic relationship of species from Biebersteiniaceae and Nitrariaceae based on complete chloroplast genomes[J]. Plants,2020,9(11) :1605.

[28]CHENQY,CHENCL,WANGB,etal.Complete chloroplast genomes of 11 Sabia samples :genomic features,comparative analysis,and phylogeneticrelationship[J].Frontiers in PlantScience, 2022,13:1052920.

[29]FENG Z, ZHENG Y, JIANG Y,et al. Phylogenetic relationships, selective pressure and molecular markers development of six species in subfamily Polygonoideae based on complete chloroplast genomes[J]. Scientific Reports,2024,14(1) :9783.

[30]KIM B,KIMJ,PARK H,et al. The complete chloroplast genome sequence of Bienertia sinuspersici[J].Mitochondrial DNA PartB, Resources,2016,1(1) :388-389.

[31]CHEN M M, ZHANG M,LIANG Z S,et al. Characterization and comparative analysis of chloroplast genomes in five Uncaria species endemic to China[J].International Journal of Molecular Sciences, 2022,23(19) :11617.

[32]JINGZ,LIWJ,SONGF,et al.Comparative analysis of complete Artemisia subgenus Seriphidium(Asteraceae:anthemideae) chloroplast genomes:insights into structural divergence and phylogenetic relationships[J].BMC Plant Biology,2023,23(1):136.

[33]YAN X K,LIU TJ,YUAN X,et al. Chloroplast genomes and comparative analyses among thirteen taxa within Myrsinaceae s.str. clade(Myrsinoideae,Primulaceae)[J].International Journal of Molecular Sciences,2019,20(18):4534.

[34]ZHAOW,GUOLR,YANGY,et al.Complete chloroplast genome sequences of Phlomis fruticosa and Phlomoides strigosa and comparative analysis of the genus Phlomis sensu lato (Lamiaceae) [J].Frontiers in Plant Science,2022,13:1022273.

[35]KAILA T, CHADUVLA P K,RAWAL H C,et al. Chloroplast genome sequence of clusterbean(Cyamopsis tetragonoloba L.) : genome structure and comparative analysis[J]. Genes,2O17,8(9): 212.

[36] SHI HW, YANG M, MO C M ,et al. Complete chloroplast genomes of two Siraitia Merrill species:comparative analysis,positive selection and novel molecular marker development[J].PLoS One, 2019,14(12) :e0226865.

[37]RAUBESONLA,PEERYR,CHUMLEYTW,et al.Comparative chloroplast genomics:analyses including new sequences from the angiosperms Nuphar advena and Ranunculus macranthus[J]. BMC Genomics,2007,8:174.

[38]SOUZAU JB,NUNES R,TARGUETAC P,et al. The complete chloroplast genome of Stryphnodendron adstringens (Leguminosae - Caesalpinioideae) :comparative analysis with related Mimosoid species[J].Scientific Reports,2019,9(1):14206.

[39]JIANGHH,WASEEMM,WANGY,etal.Development of simple sequence repeat markers for sugarcane from data mining of expressed sequence tags[J].Frontiers in Plant Science,2023,14: 1199210.

[40]ASAF S,WAQAS M,KHAN A L, et al. The complete chloroplast genome of wild rice (Oryza minuta)and its comparison to related species[J].Frontiers inPlant Science,2O17,8:304.

[41]ASAF S,KHAN AL,KHAN A,et al. Unraveling the chloroplast genomes of two Prosopis species to identify its genomic information, comparative analysesandphylogenetic relationship[J]. International Journal of Molecular Sciences,2020,21(9):3280.

[42] CHEN S Y, ZHANG H, WANG X, et al. Analysis of codon usage bias in the chloroplast genome of Helianthus annuus J-O1[J]. IOP Conference Series:Earth and Environmental Science,2O21,792 (1) :012009.

[43]CAI ZQ,PENAFLORC,KUEHL JV,et al.Complete plastid genome sequences of Drimys,Liriodendron,and Piper;implications for the phylogenetic relationships of magnolids[J].BMC Evolutionary Biology,2006,6:77.

[44]WANG X S,WANG YQ,LI SH, et al. Analysis of codon usage bias in the Platycarya chloroplast genome[J]. Tree Genetics and Molecular Breeding,2021,11(1):1-11.

[45]YANG X,LUO X N, CAI X P. Analysis of codon usage pattern in Taenia saginata based on a transcriptome dataset[J].Parasitesamp; Vectors,2014,7 ;527.

[46]WANGR,LAN Z,LUO YJ,et al.The complete chloroplast genome of Stachys geobombycis and comparative analysis with related Stachys species[J].Scientific Reports,2024,14:8523.

[47] CHEN J, MA WQ,HU X W,et al. Synonymous codon usage bias in the chloroplast genomes of 13 oil-tea Camellia samples from South China[J].Forests,2023,14(4):794.

[48]趙月梅,徐其碧,楊貴清,等.艾納香葉綠體基因組密碼子使用 偏性分析[J].西部林業科學,2023,52(3):55-62,77.

[49]XU C,CAI X N,CHEN Q Z,et al. Factors affecting synonymous codon usage bias in chloroplast genome of Oncidium Gower Ramsey [J].Evolutionary Bioinformatics Online,2011,7:271-278.

「50]楊洪升.謝平.李麗麗.等.杠板歸葉綠體基因組密碼子偏好 性分析[J/OL].分子植物育種[2025-03-05].https://link. cnki.net/urlid/46.1068.S.20230314.1714.026.

[51]XUPR,ZHANGLJ,LULP,etal.Patternsingenome-wide codonusage bias in representative speciesoflycophytesand ferns [J].Genes,2024,15(7):887.

[52]FUY,LIANGFS,LICJ,etal.Codon usagebiasanalysis in macronuclear genomes of ciliated protozoa[J].Microorganisms, 2023,11(7):1833.

[53]晁岳恩,吳政卿,楊會民,等.11種植物psbA基因的密碼子偏 好性及聚類分析[J].核農學報,2011,25(5):927-932.

[54]PINGJ,ZHONGXN,WANGT,etal.Structuralcharacterization ofTrivalvaria costata chloroplast genome and molecular evolution of rps12geneinmagnoliids[J].Forests,2023,14:1101.

[55]MORTONBR,SO B G.Codon usage in plastid genes is correlated withcontext,positionwithinthegene,and aminoacid content[J]. Journal ofMolecularEvolution,2000,50(2):184-193.

[56]GAO MQ,HUO X W,LUL T,et al. Analysis of codon usage patterns in Bupleurum falcatum chloroplast genome[J].Chinese Herbal Medicines,2023,15(2):284-290.

[57]YANGJB,YANGSX,LIHT,et al.Comparativechloroplast genomes of Camellia species[J].PLoS One,2013,8(8) :e73053.

[58]HAOBQ,XIAYY,ZHANGZY,etal.Comparative analysisof the complete chloroplast genome sequences of four Camellia species[J].BrazilianJournal ofBotany,2024,47(1) :93-103.

[59]LINP,YIN HF,WANGKL,et al. Comparative genomic analysis uncovers the chloroplast genome variation and phylogenetic relationships of Camellia species[J].Biomolecules,2022,12(10): 1474.

[60]YEXM,HUDN,GUOYP,et al.Complete chloroplast genome ofCastanopsissclerophylla(lindl.)schott:genome structureand comparativeand phylogenetic analysis[J].PLoSOne,2O19,14 (7):e0212325.

[61]WU ZH,LIAOR,YANGTG,et al.Analysis of six chloroplast genomes provides insight into the evolution of Chrysosplenium (Saxifragaceae)[J].BMCGenomics,2020,21(1):621.

[62]張文娟.基于密碼子水平的生物信息學分析及進化研究[D]. 上海:復旦大學,2006.

[63]胡振民,萬青,李歡,等.茶樹CsNRTI.1基因密碼子使用 特性分析[J].江蘇農業學報,2019,35(4):896-903.

[64]趙振寧,孫浩田,宋雨茹,等.山楂屬植物葉綠體基因組特征與 密碼子偏好性分析[J].江蘇農業學報,2023,39(2):504-517.

(責任編輯:徐艷)

收稿日期:2025-04-28

基金項目:省自然科學基金項目(2023JJ50465);省科技創新計劃項目(2024RC8289):市科技創新計劃項目(2023-RC3501);國家級大學生創新訓練項目(S202310553022)

作者簡介:曾文娟(2004-),女,邵陽人,本科生,主要從事茶學研究。(E-mail)2103562420@qq.com。周偉為共同第一作者。

通訊作者:陳致印,(E-mail) 772612626@qq.com

主站蜘蛛池模板: 黄色网站在线观看无码| 久996视频精品免费观看| 日本一区二区不卡视频| 国产精品太粉嫩高中在线观看| 国产三级视频网站| 久久婷婷六月| 国产成人乱无码视频| 九九香蕉视频| 乱人伦视频中文字幕在线| 91亚洲精选| 日本在线视频免费| 久久久久久久蜜桃| 呦女亚洲一区精品| 国产99在线| 亚洲资源站av无码网址| 伊人久久精品无码麻豆精品| 国产一在线观看| 亚洲精品午夜天堂网页| 国产资源免费观看| 免费国产福利| 2021最新国产精品网站| 亚洲精品在线影院| 伊人久久婷婷五月综合97色| a级毛片免费网站| 在线精品亚洲一区二区古装| 亚洲中文字幕在线精品一区| 国产国语一级毛片| 日韩福利在线视频| 久久国产拍爱| 黄色网站在线观看无码| 久996视频精品免费观看| 色婷婷狠狠干| 日韩av高清无码一区二区三区| 亚洲综合中文字幕国产精品欧美| 精品撒尿视频一区二区三区| 91小视频版在线观看www| 91视频首页| 高潮毛片免费观看| 国产人碰人摸人爱免费视频| 亚洲av无码牛牛影视在线二区| 国产精品久线在线观看| 在线精品欧美日韩| 一区二区影院| 色噜噜在线观看| 最新亚洲人成无码网站欣赏网 | 97久久精品人人做人人爽| 2020国产精品视频| 中文成人在线视频| 色婷婷亚洲十月十月色天| 波多野结衣无码中文字幕在线观看一区二区 | 久久久久久久久亚洲精品| 亚洲国产综合自在线另类| 国产午夜人做人免费视频| 欧美综合成人| 午夜不卡视频| 久久黄色视频影| 成人在线不卡| 毛片在线播放网址| 欧美中文字幕一区| jijzzizz老师出水喷水喷出| 中文国产成人精品久久一| 国产91丝袜在线播放动漫| 99久久这里只精品麻豆| 国产精品刺激对白在线| 亚洲男人在线| 婷婷六月天激情| 国外欧美一区另类中文字幕| 一本大道香蕉中文日本不卡高清二区| 91网红精品在线观看| 久久99国产综合精品女同| 国产激爽大片高清在线观看| 欧美激情首页| 欧美精品1区| 亚洲视频影院| 一本久道久综合久久鬼色| 色综合狠狠操| 国产高清在线精品一区二区三区| 亚洲中文字幕日产无码2021| 免费 国产 无码久久久| 国产精品成人啪精品视频| 中国国产A一级毛片| 四虎永久在线|