基于转录组数据挖掘牦牛皱胃发育代谢的关键候选基因

刘益丽,唐 娇,闵 奇,杨 露,王泽宁,胡 莲,赵 迪,江明锋*

(1.西南民族大学畜牧兽医学院,成都 610041;2.西南民族大学青藏高原研究院,成都 610041)

牦牛是长期生活在我国青藏高原上的特有物种,可为牧民提供生产及生活所必需的资源,具有极高的经济价值。犊牦牛出生后,从液体奶到断奶采食草料经历了巨大的生理和代谢变化,消化道系统的结构和功能迅速发育。牦牛皱胃是典型的管状器官[1],是牦牛的第四个胃,又称为真胃,刚出生犊牛的消化吸收功能主要依靠皱胃和小肠来完成[2]。新生犊牛皱胃已有一定程度的发育,具有消化吸收能力,摄食液体乳时,经食管沟可直接进入皱胃消化,利用葡萄糖为机体供能[3-4],犊牛出生2～3周内,皱胃是最主要的营养消化场所;成年时期的皱胃主要是对前胃消化的较小食糜进一步消化吸收[5]。由此可见,皱胃在牦牛的整个生长发育过程中发挥着重要的生理功能。

皱胃早期相关研究主要集中在疾病、引发疾病因素、形态及其影响因素[6-8]。近年研究主要在皱胃相关酶[9-10]、微生物及其作用等方面[11-12]。有研究对犊牛皱胃凝乳酶和胃蛋白酶安全性进行了评估,结果表明皱胃凝乳酶和胃蛋白酶可作为食品酶,用于牛奶加工和奶酪生产[10]。周亚楠等[13]研究发现,哺乳期补饲开食料能正向调控牦牛犊牛的皱胃形态与功能发育,进而提高对营养物质的消化利用率。Nazish等[11]从山羊和绵羊皱胃分离出睾酮科单胞菌、威氏假单胞菌和线虫等,研究3种皱胃真菌对线虫进行生物防治,3种皱胃真菌培养滤液浓度增加,线虫的死亡率增加,对线虫卵孵化具有抑制作用。Kasl等[12]研究发现,含有醋酸钠的口服电解质溶液可限制降低犊牛皱胃内大肠杆菌的体外生长潜力,可能减少皱胃大肠杆菌过度到肠道,从而减少腹泻发病率[14-15]。可见皱胃的深入研究对反刍动物的健康发育和相关产品开发具有重要的意义。

对牦牛复胃的转录组研究主要集中在瘤胃,网胃、瓣胃和皱胃的相关研究较少。本团队谢书琼等[16]前期对牦牛皱胃全长转录组进行了测定和初步分析,但基于不同年龄段牦牛皱胃组织的转录组分析和功能发育研究鲜见报道。因此,本研究在前期基础上,以液体乳为主阶段(1、20日龄),过渡阶段(60日龄),反刍为主阶段(15月龄和3岁)共5个年龄段的牦牛皱胃组织为研究对象,进行RNA-Seq高通量测序及功能分析,筛选参与皱胃发育和代谢调控的信号通路及相关基因,为今后阐明牦牛皱胃发育和营养代谢的调控分子机制提供更多的理论依据和重要参考,为皱胃组织相关产品开发提供新思路。

1 材料与方法

1.1 试验动物与样品采集

动物样本采自四川省阿坝藏族羌族自治州红原县,牦牛培育采用自然放牧方式,犊牛出生后随母哺乳,自由采食。本试验选取5个年龄段,即1日龄、20日龄、60日龄、15月龄与3岁的健康牦牛作为研究对象,每组各3头。按照以下程序人道地处理牦牛:放血前电休克,昏迷时经颈部放血并快速解剖,用灭菌剪刀取皱胃上皮及皮下组织块,剪成约1 cm×1 cm×0.5 cm的小块,用DEPC水冲洗后迅速放入冻存管中,并立即投入液氮中保存。

1.2 皱胃指数测量

准确测量每头牦牛体重,屠宰后迅速取出完整复胃,将瘤胃、网胃、瓣胃和皱胃相连部位结扎,将4个胃室分开,清空内容物,准确测得4个胃室重量并记录。皱胃指数为皱胃重(g)/复胃重(g),数据以“均数±标准差(Mean± SD)”表示,使用Microsoft Excel 2007与Origin9.0软件进行数据采集、记录、运算与分析。两组数据间的比较用t检验方法,多组数据间的差异分析使用ANOVA检验方法。P<0.05表示有显著性差异,P<0.01表示有极显著性差异。

1.3 主要试剂及仪器

1.3.1 主要试剂 mirVanaTMmiRNA ISOlation Kit 试剂盒购自Invitrogen公司(AM1561);TruSeq Stranded mRNA Sample Prep Kit购自Illumina公司(RS-122-2101);Agencourt AMPure XP核酸纯化试剂盒购自Beckman Coulter公司(A63881);QubitTMRNA HS Assay Kit 购自Thermo公司(Q32852);Qubit dsDNA Assay Kit购自Thermo公司(Q328520);Bioanalyzer 2100 RNA-6000 Nano Kit购自Aglient(5067-1511);Bioanalyzer 2100 DNA-1000 Kit购自Aglient(5067-1504);Super Script II Reverse Transcriptase购自Invitrogen公司(18064014);PrimeScript RT reagent Kit with gDNA Eraser购自TaKaRa公司(RR047B); TB Green Premix Ex TaqTMII (Tli RNaseH Plus) 购自TaKaRa公司(RR820A)。

1.3.2 主要仪器 台式离心机购自Eppendorf公司(Centrifuge 5418R);PCR 仪购自Bio-rad公司(MyCycler);紫外分光光度计购自Thermo公司(NanoDrop 2000);实时荧光定量PCR仪购自Roche公司(LightCycler 96);凝胶成像分析系统购自天能公司(Tanon 2500);定量仪购自Invitrogen公司(Qubit2.0);磁力架购自Invitrogen公司(Magneticstand-96);Bioanalyzer生物分析仪购自Aglient(2100)。

1.4 RNA提取与 cDNA文库构建

根据mirVanaTMmiRNA ISOlation Kit试剂盒操作流程提取15个牦牛皱胃组织的总RNA,利用Nanodrop 2000超微量分光光度计检测总RNA的浓度、纯度(OD260 nm/OD280 nm和OD260 nm/OD230 nm比值)和RNA完整值(RNA integrity number,RIN),随后用1%琼脂糖凝胶电泳进行检测,将质检合格的RNA置于-80 ℃超低温冰箱中保存备用。使用Illmina® TruSeq Stranded mRNA样品制备试剂盒构建15个测序样本的cDNA文库。主要步骤包括:将样品总RNA用DNase消化DNA,用Oligo(dT)磁珠富集mRNA;用打断试剂将mRNA打断成短片段,以打断后的mRNA短片段为模板,用六碱基随机引物合成一链cDNA,利用二链合成反应体系合成二链cDNA并用试剂盒纯化双链cDNA;对纯化的双链cDNA进行末端修复、加A尾并连接测序接头,片段大小选择和PCR扩增;利用Agilent 2100生物分析仪检测构建文库的大小、浓度及质量,得到合格cDNA文库。

1.5 转录组测序与质量控制

本试验中15个样本的建库测序及部分数据分析由上海欧易生物医学科技有限公司完成。基于Illumina二代高通测序平台(Illumina HiSeq X ten),采用PE150测序方式对文库进行测序,获得原始序列数据(raw reads)。为了得到可用于后续分析的高质量reads,使用Trimmomatic软件[17]进行质控并去除接头(adaptor),过滤掉低质量碱基及N碱基得到高质量的clean reads,统计clean reads的Q30和GC含量。用 Hisat2[18]将clean reads与牦牛参考基因组(GCF_000298355.1)进行序列比对,通过基因组比对率来评估样本的情况。

1.6 主成分分析

利用HTSeq-count[19]软件获取落到各个样本中基因的reads 数目,使用cufflinks[20]软件对蛋白编码基因的FPKM[21]表达量值定量,用R语言绘制主成分分析(principal component analysis,PCA)二维图。

1.7 差异表达基因筛选及功能富集分析

利用HTSeq-count[19]软件获得比对到每个样本中蛋白编码基因的reads数,使用DESeq (2012) R package[22]的函数estimateSizeFactors对数据进行标准化,使用函数nbinomTest计算差异倍数(Fold Change)和差异的显著性(P-value)。选取P<0.05,log2(Fold Change)>2 的差异表达基因(DEGs)进行GO和KEGG[23]通路富集分析。以1日龄样本为参照,利用GO Seq分别对20日龄、60日龄、15月龄和3岁龄的牦牛皱胃DEGs进行GO功能注释,基于GO数据库功能,按照细胞组分(cellular component,CC),生物过程(biological process,BP)和分子功能(molecular function,MF)将4个年龄段的DEGs进行分类,以P<0.05为标准筛选得到显著富集的GO条目和KEGG信号通路。

1.8 qRT-PCR验证

为进一步验证牦牛皱胃RNA-Seq数据的可靠性,随机选取5个差异表达基因,利用Primer BLAST设计引物,采用实时荧光定量PCR验证基因的表达水平。将测序返样RNA用PrimeScript RT reagent Kit with gDNA Eraser逆转录试剂盒反转录为cDNA,以GAPDH作为内参基因,利用TB Green® Premix Ex TaqTMII (Tli RNaseH Plus)试剂盒在LightCycler 96 系统上检测。荧光定量PCR反应体系为20 μL:2×TB Green® Premix Ex TaqTMII 10 μL,上、下游引物(10 μmol·L-1)各0.5 μL,皱胃组织cDNA 2 μL,RNase Free ddH2O 7 μL。PCR 程序: 95 ℃ 30 s;95 ℃ 5 s,60 ℃ 30 s,40个循环;95 ℃ 5 s, 60 ℃ 1 min;重复3次。采用2-△△ct法分析数据,用t-检验统计分析相对表达量,数据结果表示为“平均数士标准差(Mean±SD)”,P<0.05为差异显著。对应的引物序列见表1。

表1 基因及对应引物序列Table 1 Gene and the corresponding primer sequences

2 结 果

2.1 不同年龄段牦牛皱胃重量和指数分析

对牦牛5个年龄段15个样本的复胃重量和皱胃指数进行分析,如表2所示,在同一年龄段,1日龄皱胃最重,瘤胃次之,皱胃重显著高于网胃和瓣胃重(P<0.05);20日龄皱胃重显著高于网胃和瓣胃重(P<0.05),与瘤胃重差异不显著(P>0.05);60日龄瘤胃重最大,显著高于皱胃、网胃和瓣胃重(P<0.05);15月龄瘤胃重显著高于瓣胃(P<0.05),皱胃与网胃重最小,差异不显著(P>0.05);3岁瘤胃重量最大,显著高于瓣胃(P<0.05),皱胃重大于网胃重,差异不显著(P>0.05)。如图1所示,1日龄、20日龄和60日龄牦牛皱胃重量差异不显著(P>0.05),15月龄与3岁牦牛的皱胃重显著增加(P<0.05),皱胃重整体呈现增加趋势;1日龄、20日龄牦牛的皱胃指数差异不显著(P>0.05);与1日龄相比,15月龄与3岁牦牛的皱胃指数显著下降(P<0.05),皱胃指数与皱胃重量变化趋势相反,整体呈现下降趋势。

表2 不同年龄段牦牛复胃重量的变化Table 2 Changes of compound stomach weight of yaks at different ages

横坐标表示年龄,N=3/组;左纵坐标表示皱胃重量,右纵坐标表示皱胃指数 The abscissa represents the age, N=3/group, the ordinate on the left side represents the weight of abomasum, and the ordinate on the right side represents the abomasal index图1 不同年龄组牦牛皱胃重量与指数分析Fig.1 Analysis of weight and index of abomasum in yak at different ages

2.2 不同年龄段牦牛皱胃总RNA质量与测序数据质量控制

15个皱胃组织样本的总RNA浓度均大于300 ng·μL-1,RIN值大于7.5,OD260 nm/OD280 nm在2.0～2.2之间,28S/18S核糖体RNA条带清晰,RNA完整性较好,RNA质量检验合格,可进一步用于后续试验分析。利用Agilent 2100生物分析仪检测cDNA文库的片段大小与质量符合要求,cDNA文库有效浓度大于4 nmol·L-1,表明cDNA文库构建成功。

如表3所示,15个皱胃组织样本的有参转录组测序获得过滤后的clean reads数目为741.36 M,过滤后得到的测序量106.19 G,有效碱基百分比均在93%以上,Q30均在95%以上,GC含量在50%左右,比对到参考基因组的reads比率均在94%以上,表明15个样本的测序质量均较高,数据准确可靠,所选参考基因组的组装满足后续分析要求。

表3 牦牛皱胃组织转录组数据质量检测分析Table 3 Quality analysis of transcriptome sequencing data of abomasum tissues in yak

2.3 不同年龄段牦牛皱胃基因表达水平分析

对不同年龄段牦牛皱胃组织基因表达水平进行分析,5个年龄组间的基因表达水平与基因密度均存在差异。如图2A所示,1日龄、20日龄和60日龄3个年龄组组间的基因表达水平差异较小,15月龄组的基因表达水平低于其他组,3岁组整体基因表达量高于其他组;如图2B所示,20日龄组基因密度高于其他组。

A.不同年龄组牦牛皱胃基因表达水平箱线图,横坐标为年龄,纵坐标为log10(FPKM+1),且展示了最大值、上四分位数、中位数、下四分位数和最小值。B.不同年龄组牦牛皱胃基因密度分布图,横坐标为log10(FPKM),纵坐标为基因密度 A. Boxplots of gene expression levels in abomasum of yaks at different ages, with the abscissa of age and the ordinate of log10(FPKM+1), showing the maximum, upper quartile, median, lower quartile and minimum values. B. Distribution of abomasal gene density in yaks at different ages, log10(FPKM)in abscissa and gene density in ordinate图2 不同年龄组牦牛皱胃基因表达水平箱线图和基因密度分布图Fig.2 Boxplots of gene expression levels and distribution of gene density in abomasum of yaks at different ages

2.4 不同年龄段牦牛皱胃差异表达基因分析

为进一步探讨在5个年龄阶段牦牛皱胃组织DEGs的表达模式,本研究采用了两种比较方法[24],即封闭组 (20日龄vs.1日龄,60日龄vs.1日龄,15月龄vs.1日龄和3岁vs.1日龄),相继时间组 (20日龄vs.1日龄,60日龄vs.20日龄,15月龄vs.60日龄和3岁vs.15月龄)。如表4所示,以1日龄组为对照,随年龄增加DEGs逐渐增多,在20日龄、60日龄、15月龄和3岁组中分别鉴定到1 310、1 715、1 931和2 199个DEGs,且4个对比组中上调基因的数量均高于下调基因(表4),其特有DEGs分别为 341、242、395和565个,4个组共有基因565个(图3A),包括GKN1、CXCL17、SCNN1B、SCNN1G和CCL5等基因。以前一个时间点为对照,20日龄、60日龄、15月龄和3岁组分别鉴定到1 310、861、569和597个DEGs,20日龄vs.1日龄组的DEGs明显高于其他3个对比组,其特有DEGs分别为 950、419、250和284个,4个组共有基因9个(图3B),分别是MMP9、CDH26、GSTA3、TMEM61、ADGRG5、LOC102269721、LOC102287938、LOC102276015和LOC102271445。

表4 不同年龄段牦牛皱胃组织差异表达基因数目统计Table 4 Statistics of the number of DEGs in abomasum tissues of yaks at different ages

图3 不同年龄段牦牛皱胃差异表达基因韦恩图,4个封闭组(A),4个相继组(B)Fig.3 Venn diagram of DEGs in abmasum in yaks with 4 closed (A) and 4 consecutive (B) groups at different ages

2.5 不同年龄段牦牛皱胃差异表达基因的功能注释及分类

对封闭组的DEGs进行GO富集分析,在20日龄、60日龄、15月龄和3岁组中分别有1 191、2 578、1 117和2 835条显著性条目。按P值排序筛选BP、CC和MF三大类的前10条显著性GO条目。如图4所示,不同年龄组的前10条显著性GO条目中有共有的,也有各年龄阶段特有的条目。在生物过程中,总共19个条目,共有条目12个,其中20日龄、60日龄、15月龄和3岁组共有3个条目,分别是免疫反应、炎症反应和适应性免疫反应;20日龄、15月龄和3岁年龄组特有条目分别有1、4和2个,20日龄富集到补体激活的调节显著性条目,60日龄无特有条目,15月龄分别富集到T细胞受体信号通路、细胞溶质钙离子浓度的正向调节、Fc-ε受体信号通路等显著性条目,3岁组是细胞趋化性和干扰素-γ介导的信号通路显著性条目。在细胞组分中,总共14个条目,共有条目10个,其中20日龄、60日龄、15月龄和3岁组共有5个条目,分别富集于质膜外侧、细胞外区、细胞外间隙和质膜等显著性条目;20日龄、15月龄和3岁组各有1个特有条目,分别富集于胶原三聚体、钠通道复合物、丝足和含胶原蛋白的细胞外基质显著性条目。在分子功能中,总共17个条目,共有条目10个,其中20日龄、60日龄、15月龄和3岁4个年龄段共有6个条目,分别富集于抗原结合、C-C趋化因子、C-C趋化因子受体活性、趋化因子结合等显著性条目;20日龄、60日龄、15月龄和3岁4个年龄段特有条目分别为3、1、2和1个,20日龄组富集到双链RNA结合、免疫球蛋白受体结合等条目,60日龄组富集到乙醇脱氢酶(NADP+)活性条目,15月龄组是甲状腺激素结合和跨膜信号受体活性条目,3岁组富集到具有拉伸强度的细胞外基质结构成分条目。

2.6 不同年龄段牦牛皱胃差异表达基因KEGG富集分析

对封闭组的差异表达基因分别进行KEGG功能富集分析,在20日龄、60日龄、15月龄和3岁组中分别有100、102、100和111条显著通路。按P值排序筛选4个年龄组的前30条显著性KEGG通路(图5)。结果显示,每个年龄组前30条显著通路同样存在共有的,也有各年龄阶段特有的,总共43条通路,共有通路33个,其中20日龄、60日龄、15月龄和3岁4个年龄段共有20个条目,包括吞噬体、同种异体移植排斥反应和细胞因子-细胞因子受体相互作用;20日龄、60日龄、15月龄和3岁4个年龄段特有信号通路分别是4、1、3和4条,如20日龄组富集到白细胞跨内皮迁移通路,15月龄组富集到亚油酸代谢、谷胱甘肽代谢和糖尿病并发症中的AGE-RAGE信号通路,3岁组富集到ECM受体相互作用、细胞色素P450对外源性物质的代谢、药物代谢-细胞色素P450等信号通路。

2.7 牦牛皱胃发育和代谢候选基因筛选

本研究以1日龄为对照,筛选20日龄、60日龄、15月龄和3岁组皱胃组织DEGs,根据显著性高低,分别从4个年龄组前30个DEGs中筛选影响牦牛皱胃发育的相关候选基因,见表5。结果显示,在20日龄、60日龄、15月龄和3岁的年龄组中均包含共有基因,如KRT20;在20日龄、15月龄和3岁的年龄组中均包含共有基因,如GKN1、CXCL17;在60日龄、15月龄和3岁的年龄组中均包含共有基因,如JCHAIN、SCNN1B、SCNN1G;60日龄和3岁组均包含CCL5、B2M;另外4个组也有一些独有基因,如20日龄中的GPR34、PTGS2、LAPTM5,60日龄组中的SCNN1B,15月龄组中的IGF2BP3,3岁龄组中的IGF2BP2,这些基因参与牦牛皱胃的不同发育阶段,可能是与牦牛皱胃发育相关的重要候选基因。另外,本研究富集到参与皱胃的葡萄糖代谢、葡萄糖转运、脂肪酸转运、肽转运、胃蛋白酶和溶菌酶活性相关的信号通路,分别筛选到12、12、16、4、1和11个候选基因,见表6。其中GANAB、PRKCSH、GBA2、SLC2A1、SLC2A3、SLC2A4、CPT1B、GOT2等是与牦牛皱胃的代谢和营养吸收相关的重要候选基因。

A-C. 4个年龄组DEGs在BP、CC、MF三大类中前10个GO条目。气泡大小表示DEGs的数量(计数),气泡颜色表示错误发现率(FDR P-value),x轴为-lg(P-value)值,y轴为GO条目 A-C. The top 10 GO terms in molecular function, cellular component and biological process, respectively. The bubble size represents the number (count) of DEGs, the bubble color represents false discovery rate (FDR P-value), the abscissa represents -lg(P-value), and the ordinate displays GO term图4 差异表达基因GO富集气泡图Fig.4 Bubble charts of GO enrichment of differentially expressed genes

2.8 差异表达基因验证

对随机选取的5个差异基因进行qRT-PCR验证,结果如图6所示,以1日龄为对照,ABAC1基因在20日龄和15月龄上调;ACAT1基因在15月龄上调,但在20日龄、60日龄和3岁组下调;G6PC3、GOT2和SLC38A2基因在20日龄、60日龄、15月龄和3岁组均下调。上述所挑选的差异表达基因的qRT-PCR变化规律与转录组测序数据表达规律基本一致,进一步表明测序结果具有可靠性。

3 讨 论

皱胃发育与多种因素密切相关,如饲养方式[25]、日粮[5,13]等,对反刍动物的营养物质消化吸收和生长发育具有重要影响。近年来,已有研究对牛[26-28]、羊[5]、骆驼[29]等的皱胃组织进行形态鉴定和生理功能研究,但这些研究主要集中在皱胃组织形态、免疫方面等,对分子水平的功能研究较少。目前关于牦牛皱胃发育的分子机理尚不清楚,需进一步研究,为牦牛的营养调控和健康生长提供理论参考和新思路。本研究中,随年龄增长牦牛皱胃重整体呈增长趋势,但20和60日龄的皱胃重差异不显著,且60日龄有下降趋势,可能是因为20～60日龄是自然放牧犊牦牛生长发育的一个重要阶段,正在经历由以液体奶为主到以牧草为主的过渡,瘤胃快速增长以适应这一变化,但网胃、瓣胃和皱胃重变化均较小。15月龄和3岁的皱胃重显著增长,但皱胃指数逐渐下降,提示随着年龄增加皱胃组织发育更加完善,但功能方面可能会发生一定程度的变化。1日龄牦牛皱胃重在复胃重中占比最大,显著高于其他3个前胃(P<0.05),可侧面反映牦牛皱胃可能成熟比较早,且新生牦牛皱胃在复胃中发挥着重要作用。20日龄至15月龄期间,皱胃重占比相较于前一个年龄均显著下降(P<0.05),侧面反映20日龄后牦牛的能量提供不再完全由皱胃组织主导,而部分由快速发育的瘤胃进行VFA代谢所提供的营养替代,与本实验室唐娇等[30]前期对5个年龄段牦牛瘤胃的转录组研究结果相呼应。

表5 不同年龄组牦牛皱胃的前30个差异表达基因Table 5 The top 30 differentially expressed genes in abomasum of yaks in different age groups

表6 参与皱胃主要代谢功能的候选通路与基因Table 6 The pathways and genes involved in major metabolic functions of abomasum

图6 差异表达基因的实时荧光定量PCR验证Fig.6 qRT-PCR verification of DEGs

牦牛皱胃胃壁是由浆膜、肌层、黏膜层和黏膜下层构成,其中黏膜又可细分为上皮、固有层和黏膜肌层[31]。研究发现,牦牛皱胃的腺区固有层均分布着大量弥散淋巴细胞和孤立淋巴小结,使其具有比其他反刍动物更强的黏膜免疫功能[1],这可能与其长期生长在高原环境密切相关。KEGG富集分析显示,与1日龄相比,在20日龄、60日龄、15月龄和3岁4个年龄段的DEGs共同显著富集到20条信号通路,这些通路多与牦牛皱胃组织发育及免疫调控相关,包括吞噬体、细胞粘附分子、细胞因子-细胞因子受体、趋化因子信号通路、肠道IgA生成通路、Thl和Th2细胞分化、Th17细胞分化信号通路等。研究表明,Th17细胞释放IL-17A、IL-17F、IL-21、IL-22和IL-9,并通过它们与肠上皮细胞相互作用,从而诱发胃肠道炎症[32],皱胃与炎症性肠病之间是否具有联系、是否可通过皱胃探究牦牛炎症性肠病的发病机制,以及Th17细胞分化途径是否与其相关值得进一步研究。20日龄、60日龄、15月龄和3岁4个年龄段特有信号通路分别有4、1、3和4条,20日龄组是白细胞跨内皮迁移通路,15月龄组是亚油酸、谷胱甘肽代谢和糖尿病并发症中的AGE-RAGE信号通路,3岁组是ECM受体相互作用、细胞色素P450对外源性物质的代谢、药物代谢-细胞色素P450, 研究发现细胞色素P450外源性物质的代谢途径和ECM受体相互作用途径在维持瘤胃上皮完整性和健康中发挥作用[33],可推断这些代谢通路可能在牦牛皱胃的健康发育过程中发挥着重要作用。

根据差异表达基因的显著性高低,分别从4个年龄组前30个DEGs中筛选影响牦牛皱胃发育相关的候选基因,其中GKN1、CXCL17、SCNN1B、SCNN1 G、CCL5和IGF2BP3等基因可能与牦牛皱胃发育相关。胃动蛋白1(gastrokine-1,GKN1)是胃动蛋白家族成员之一,该基因能在胃粘膜中特异性表达一种新型自分泌/旁分泌蛋白。研究发现,GKN1具有保护胃窦粘膜,并能促进损伤胃窦粘膜的修复和增殖分化,促进胃粘膜上皮细胞分裂与移动、维持上皮层完整性和抑癌等功能[34-36]。趋化因子17(CXC chemokine ligand 17,CXCL17)是一种与粘膜内环境稳定相关的新型CXC趋化因子,该基因在支气管和胃部均有表达,且前者表达高于胃部[37-38]。研究发现CXCL17参与胃癌的发生[39],能促进胃癌细胞的增殖[40]。另外,在对幽门螺旋杆菌根除者的胃内菌群特征研究中发现,CXCL7可能影响胃内菌群结构[41],与胃部健康息息相关。钠通道上皮1亚基β(sodium channel epithelial 1 subunit beta,SCNN1B)位于染色体16p12.2-p12.1,编码640个氨基酸[42-43]。研究发现,SCNN1B基因能诱导细胞凋亡、阻滞细胞周期、抑制癌细胞侵袭、增殖与迁移和促进细胞粘附,其过表达可抑制胃癌细胞上皮间质转化[44-45]。钠通道上皮1亚基γ(sodium channel epithelial 1 subunit gamma,SCNN1G)主要在胃癌细胞中的细胞质和细胞膜中,研究表明SCNN1G能显著抑制胃癌细胞增殖、侵袭与迁移,该基因的表达与WNT/β-catenin信号通路紧密相关[46]。C-C基序趋化因子配体5(C-C motif chemokine ligand 5,CCL5)是CC亚家族趋化因子基因,具有调节正常T细胞表达和分泌的功能,同时能促进炎性细胞因子的释放,参与免疫调节和炎症过程[47];CCL5也能提高活化T细胞的葡萄糖摄取能力和ATP水平[48]。并在瘤胃上皮的基底和棘突层定位到CCL5,在瘤胃内环境稳定中发挥重要作用[30,49]。胰岛素样生长因子II mRNA结合蛋白3(insulin like growth factor 2 mRNA binding protein 3,IGF2BP3)是一种RNA结合蛋白,能参与细胞极化、迁移、增殖和分化等生物过程[50-51]。综上,GKN1、CXCL17、SCNN1B、SCNN1G、CCL5和IGF2BP3基因可能在牦牛皱胃组织的发育过程中发挥着关键调控作用。

牦牛皱胃可消化的营养物质包含葡萄糖、脂质、蛋白质等,本研究富集到参与皱胃的葡萄糖代谢、葡萄糖转运、脂肪酸转运、肽转运、胃蛋白酶和溶菌酶活性相关的信号通路,分别筛选到12、12、16、4、1和11个候选基因。其中GANAB、PRKCSH、GBA2、SLC2A1、SLC2A3、SLC2A4、CPT1B、GOT2、SLC15A1等是与牦牛皱胃代谢和营养吸收相关的重要候选基因。葡萄糖酰神经酰胺酶2(glucosylceramidase beta 2,GBA2)基因编码一种胆汁酸β-葡萄糖苷酶,可催化胆汁酸3- o -葡萄糖苷作为内源性化合物的水解,在碳水化合物的运输和代谢中起作用[52]。溶质载体家族2成员1(solute carrier family 2 member 1,SLC2A1)又称GLUT1,是反刍动物促进葡萄糖摄入的主要转运蛋白[53]; 溶质载体家族27成员4(solute carrier family 27 member 4,SLC27A4)又称FATP4,是动物机体调控脂肪酸摄取的重要因子,能促进机体对脂质的吸收[54];溶质载体家族15成员1(solute carrier family 15 member 1,SLC15A1)又称PepT1,参与小肽的跨膜运转,可直接促进机体对蛋白质的消化吸收[55]。研究发现,补饲精料更有利于犊牦牛对葡萄糖、脂肪、蛋白质的消化吸收,同时发现牦牛皱胃GLUTI、PepT1和FATP4基因相对表达量均显著升高[13,56],说明这些基因在牦牛皱胃的营养代谢和吸收过程中扮演着关键的角色。目前国内外对于牦牛皱胃营养转运载体的研究极少,有待进一步探究。

4 结 论

本研究对液体乳为主阶段(1、20日龄)、过渡阶段(60日龄)、反刍为主阶段(15月龄和3岁龄)共5个年龄段的牦牛皱胃组织进行了转录组测序,并对复胃相关参数进行测定。在整个发育过程中,瘤胃重的增长最快,其次是瓣胃、网胃和皱胃,且皱胃指数整体呈下降趋势。GO注释分析和KEGG富集分析结果表明,差异表达基因多参与免疫反应、吞噬体、细胞因子-细胞因子受体相互作用、ECM受体相互作用、细胞色素P450对外源性物质的代谢、药物代谢-细胞色素P450等信号通路。这些信号通路有利于牦牛在高原环境中更好促进皱胃细胞的增殖分化,提高免疫力,更有利于皱胃的生长发育。对牦牛皱胃营养代谢通路进一步筛选,得到GANAB、PRKCSH、GBA2、SLC2A1、SLC2A3、SLC2A4、CPT1B、GOT2和SLC15A1等与皱胃营养吸收和代谢相关的重要候选基因。上述这些基因可能参与牦牛皱胃的发育和营养代谢过程,可为进一步理解牦牛皱胃的发育和营养代谢机制提供理论基础。