畜牧兽医学报  2021, Vol. 52 Issue (5): 1293-1306. DOI: 10.11843/j.issn.0366-6964.2021.05.015    PDF    
幼龄和成年太行山羊附睾头差异竞争性内源RNAs机制分析
郭相前, 董复成, 黄鑫芸, 刘文忠, 乔利英, 张春香, 任有蛇     
山西农业大学动物科学学院, 太谷 030801
摘要:旨在筛选出幼龄和成年太行山羊附睾头中差异表达的长链非编码RNAs(lncRNAs)、微小RNAs(miRNAs)和mRNAs,构建太行山羊附睾头中免疫相关基因调控的竞争性内源RNAs(ceRNAs)网络。本研究选取健康状况良好、体重相近的幼龄(2月龄)和成年太行山羊(2周岁)公羊各3只,去势采集其附睾头组织进行全转录组测序,用DESeq2软件筛选出幼龄和成年太行山羊附睾头差异mRNAs、lncRNAs和miRNAs。利用miRanda软件和R-package(reshape2、dplyr、tidyr),基于ceRNA-score原理分析得到差异ceRNAs表达谱,对预测所得具有ceRNAs关系的mRNAs进行GO和KEGG富集分析,并绘制得到太行山羊附睾头免疫相关ceRNAs网络。最后,各随机挑选8个mRNAs、miRNAs和lncRNAs,并通过qRT-PCR验证全转录组测序结果的准确性。根据分析结果,以幼龄太行山羊附睾头为对照,成年山羊附睾头差异表达mRNAs有6 461个,其中上调2 997个、下调3 464个;差异表达lncRNAs共有1 147个,其中703个上调、444个下调;差异表达miRNAs共有182个,其中81个上调、101个下调。共得到具有ceRNAs调控关系的lncRNAs 366个,其中上调213个,下调153个;mRNAs有3 131个,其中1 253个上调,1 878个下调;miRNAs有140个,其中48个上调,92个下调。分析具有ceRNAs机制的基因发现,表达量显著上调的与免疫相关的mRNAs:淋巴细胞抗原6复合位点蛋白G5B(LY6G5B)、脂质运载蛋白9(LCN9)、解整合素金属蛋白酶28(ADAM28)和粘蛋白15(MUC15)基因在成年太行山羊附睾头表达量高,且极显著高于幼龄太行山羊(P < 0.01)。GO和KEGG富集分析表明,具有ceRNAs机制的差异表达基因富集在内质网蛋白加工通路、蛋白质输出通路、粘蛋白型O-聚糖生物合成通路、细胞外基质受体相互作用通路等。qRT-PCR验证结果表明,除chi-miR-320-3p外,其余差异表达的mRNAs、lncRNAs和miRNAs表达趋势与全转录组测序结果一致。附睾头免疫相关ceRNAs网络分析表明,lncRNA-MSTRG.22929.11、lncRNA-MSTRG.57822.5、lncRNA-MSTRG.26758.1、lncRNA-MSTRG.12113.3、lncRNA-MSTRG.59930.2等lncRNAs作为ceRNAs可以调控附睾头免疫相关基因表达。本研究筛选出了幼龄和成年太行山羊附睾头差异ceRNAs,挖掘并绘制了免疫相关的关键ceRNAs网络,这些lncRNAs作为ceRNAs可为太行山羊附睾头免疫调控机制研究提供参考依据。
关键词ceRNA    太行山羊    附睾头    免疫    
Analysis of the Differentially Expressed ceRNAs Mechanism in Caput Epididymis between Goat Kids and Adult Goats
GUO Xiangqian, DONG Fucheng, HUANG Xinyun, LIU Wenzhong, QIAO Liying, ZHANG Chunxiang, REN Youshe     
College of Animal Science, Shanxi Agricultural University, Taigu 030801, China
Abstract: The aim of this study was to explore and screen the differentially expressed long noncoding RNAs (lncRNAs), microRNAs (miRNAs) and mRNAs in caput epididymis between Taihang goat kids and adult Taihang goats, and to construct a competitive endogenous RNAs (ceRNAs) network associated with the regulatory mechanism of immune-related genes in caput epididymis of Taihang goats. Three 2-month-old healthy Taihang goat kids with similar body weight and three 2-year-old healthy adult Taihang goats with similar body weight were selected in this study. After castration, the caput epididymis were collected for whole transcriptome sequencing. The differentially expressed mRNAs, lncRNAs and miRNAs in caput epididymis were screened between goat kids and adult goats using DESeq2. Based on the principle of ceRNA-score, the differentially expressed ceRNAs was obtained by using miRanda software and R-package (reshape2, dplyr, tidyr). GO and KEGG enrichment analysis were used to study the function of mRNAs in the ceRNAs mechanism. The ceRNAs network related with immune function was drawn in caput epididymis of Taihang goats. Finally, 8 mRNAs, 8 miRNAs and 8 lncRNAs were randomly selected to verify the accuracy of the whole transcriptome sequencing results by qRT-PCR. According to the analysis results, compared with caput epididymis of goat kids, there were 6 461 differentially expressed mRNAs in caput epididymis of adult Taihang goats, including 2 997 up-regulated mRNAs and 3 464 down-regulated mRNAs; There were 1 147 differentially expressed lncRNAs, including 703 up-regulated lncRNAs and 444 down-regulated lncRNAs; There were 182 differentially expressed miRNAs, including 81 up-regulated miRNAs and 101 down-re-gulated miRNAs. Three hundred and sixty-six lncRNAs with ceRNAs regulatory relationship were obtained, including 213 up-regulated lncRNAs and 153 down-regulated lncRNAs; There were 3 131 mRNAs with ceRNAs regulatory relationship, including 1 253 up-regulated mRNAs and 1 878 down-regulated mRNAs; There were 140 miRNAs with ceRNAs regulatory relationship, including 48 up-regulated miRNAs and 92 down-regulated miRNAs. Among the genes in the ceRNAs mechanism, lymphocyte antigen 6 complex locus protein G5B (LY6G5B), epididymal-specific lipocalin-9 (LCN9), adisintegrinand metalloproteinase 28 (ADAM28) and Mucin 15 (MUC15) were the significantly up-regulated genes related with immune function, which expression in adult goats were higher than those in goat kids (P < 0.01). GO and KEGG analysis of the DEGs in the ceRNAs mechanism showed that of differentially expressed genes were involved in protein processing in endoplasmic reticulum, protein export, mucin type O-Glycan biosynthesis, ECM-receptor interaction and other pathways related to cell protein synthesis, secretion and transport. The results of qRT-PCR showed that the expression trends of the differentially expressed mRNAs, lncRNAs and miRNAs selected randomly were consistent with the whole transcriptome sequencing results, except for chi-miR-320-3p. In this study, the analysis results of the immune-related ceRNAs network in caput epididymis showed that lncRNA-MSTRG.22929.11, lncRNA-MSTRG.57822.5, lncRNA-MSTRG.26758.1, lncRNA-MSTRG.12113.3, lncRNA-MSTRG.59930.2 and other lncRNAs were related to immunity and sperm maturation in caput epididymis. In this study, the differentially expressed ceRNAs were screened in caput epididymais between Taihang goat kids and adult Taihang goats. The ceRNAs network related to immunity will provides a foundation for the study of lncRNAs as ceRNAs in the regulation mechanism of epididymal immune function in Taihang goats.
Key words: ceRNA    Taihang goat    caput epididymis    immune    

哺乳动物在出生时附睾处在未分化期,需要再经过一段以细胞继续增殖为主要特征的未分化期[1-2],在该期末大鼠在21日龄基本形成特殊血-附睾屏障[3],之后附睾再经过以细胞分化为主要特征的分化期和以管腔扩大为主要特征的扩张期以提供精子发育成熟所需的微环境[4],精子生成后附睾才能发挥正常的精子成熟与贮存等功能[5]。附睾不同区段血-附睾屏障结构是有差异的,附睾头部的屏障是由紧密连接和多个桥粒相连组成,定位在附睾主细胞连接处,屏障连接深层的表皮钙粘着蛋白mRNA在42日龄大鼠(性成熟)附睾头部显著高于附睾体部和尾部[6],这为以性成熟前和成年时期为代表的两个不同时期附睾头样品作为试验材料研究附睾头部基因表达调控提供了理论依据。现今,基于多组学分析技术分析竞争性内源RNA(ceRNA)机制,探明附睾不同发育期其表达调控机理对雄性繁殖生理的研究有重要意义。ceRNA是一种竞争性内源结合RNA,主要有lncRNA、circRNA、piRNA和假基因等[7-10]。基于睾丸组织的多组学测序数据,Meng等[11]构建了藻毒素-亮氨酸-精氨酸诱导后与睾丸毒理变化显著相关的lncRNAs、circRNAs、piRNAs、miRNAs和mRNAs相互调控的ceRNAs网络。另有研究表明,lncRNA-自噬相关蛋白7(lncRNA recombinant autophagy related protein 7, lncRNA ATG7)与mRNA人肺腺癌转移相关转录本1(metastasis associated in lung denocarcinoma transcript 1, MALAT1)可以通过结合miR-200a调节泛素E3连接酶(membrane-associated RING-CH7, MARCH7)的表达进而调控卵巢癌细胞迁移和侵袭[12]。lncRNA H19表达下调会通过增强miR-124-3p和整合素β 3(integrin beta 3, ITGB3)表达进而抑制异位子宫内膜中的细胞侵袭和增殖[13]。还有研究发现,lncRNA膜联蛋白A2假基因3(lncRNA annexin A2 pseudogene 3, lncRNA ANXA2P3)可以通过竞争性结合miR-613和miR-206抑制mRNA转酮醇酶(transketolase, TKT)的表达引发非梗阻性无精子症[14]。这说明lncRNAs、假基因等作为ceRNA对哺乳动物繁殖有着重要调控作用。本课题组基于基因组学技术将太行山羊附睾分为9个功能区段[15],其中,附睾头I段有大量高表达且具有多重功能的β防御素和免疫调节因子,如β防御素124(beta-defensin124, gDEFB124)、脂质运载蛋白5(lipocalin 5, LCN5)、淋巴细胞抗原6复合体G5B(lymphocyte antigen 6 complex locus protein G5B, LY6G5B)等[15-17],因此,附睾头是附睾具有免疫调控机制的特殊区段。然而,有关附睾头区段免疫的调控机理尚不清楚。

基于多组学技术系统研究不同RNAs之间的互作关系和ceRNAs机制为本研究提供了思路[18-20],因此,本试验以幼龄和成年太行山羊附睾头区段为研究对象,采用高通量测序技术进行lncRNAs、miRNAs测序并分析lncRNAs-miRNAs-mRNAs关系,筛选出差异表达ceRNAs,进一步构建附睾头与免疫调节相关ceRNAs网络。本研究可为附睾发育过程中免疫调控机制研究提供参考,为揭示山羊附睾头免疫调控的分子机制提供重要的理论依据。

1 材料与方法 1.1 试验动物和样品采集

试验动物来源于山西农业大学动物科学学院试验羊场,选取2月龄性成熟前和成年2周岁健康状况良好、体重相近的太行山羊各3只。通过去势手术取同一侧幼龄山羊附睾头(2月龄,KEH)和成年山羊附睾头(2周岁,AEH)组织作为全转录组测序样品,分别放入2 mL冻存管暂时保存在液氮罐,随后置于-80 ℃超低温冰箱保存备用。

1.2 主要试剂

PrimeScrip RT reagent kit with gDNA逆转录试剂盒(TaKaRa, Japan),TB GreenTMPremix Ex TaqTM II(TaKaRa, Japan)定量试剂盒;miRNA定量表达检测使用的M5 miRNA cDNA synthesis Kit加尾法逆转录试剂盒和M5 miRNA qPCR Assay Kit定量试剂盒(含下游引物)均购自北京聚合美生物有限公司;其他无机试剂均购自北京康达伟宏科技有限公司。

1.3 RNA提取、cDNA文库构建及Illumina测序

使用TRIzol试剂分别提取6个样品RNA,经质量控制检测合格后,分别用于smallRNA、lncRNA和mRNA文库的构建。使用NEB Next Ultra small RNA Sample Library Prep Kit for Illumina试剂盒对6个合格样品分别构建各自的small RNA文库。对同批次样品提取的RNA去除rRNA后使用Truseq RNA sample Prep Kit试剂盒进行各自lncRNA和mRNA文库构建。对构建好的文库用Illumina公司的Illumina HiSeq高通量测序平台进行测序,得到raw reads。

1.4 附睾头组织差异mRNAs、lncRNAs和miRNAs筛选

对原始数据进行质量控制后得到clean reads,与山羊参考基因组(ASM170441v1)比对和映射,在服务器上统计reads counts估计表达量,使用StringTie (v1.3.1)和DESeq2 (v1.18.0)软件进行lncRNAs、mRNAs表达定量和差异表达筛选(|log2(FC)|≥1,FDR < 0.01)。利用miRDeep2(v2.0.5)、RNAfold(v2.1.7)和Bowtie(v1.0.0)软件进行已知及新miRNAs鉴定及其表达定量,使用DESeq2软件对差异表达miRNAs进行筛选(|log2(FC)|≥1,FDR < 0.01)。

1.5 差异ceRNAs筛选与功能注释

使用miRanda(v3.3a)软件分别预测差异miRNAs的靶mRNAs和靶lncRNAs,筛选出具有与miRNAs相互调节关系的lncRNAs和mRNAs;运用R语言软件包(reshape2、dplyr、tidyr),根据RNAs表达量,使用皮尔森相关性分析确定差异表达miRNAs分别与差异表达的lncRNAs和mRNAs差异表达的皮尔森相关系数(cor<-0.8,P < 0.01),以及差异表达lncRNAs和差异表达mRNAs的皮尔森相关系数(cor>0.8,P < 0.01),筛选出具有ceRNAs统计学意义和相互调节关系的lncRNAs和mRNAs。基于ceRNA-score原理,利用自己编写的R语言脚本筛选分析得到幼龄和成年太行山羊附睾头差异ceRNAs。再对被GO和KEGG数据库注释过具有ceRNAs关系的mRNAs进行GO和KEGG富集分析,筛选显著富集结果(P < 0.05)。

1.6 实时荧光定量PCR验证测序结果

为了进一步保证测序结果的可靠性,用测序样品提取的RNA为PCR扩增模板,对随机挑选的8个mRNAs、8个miRNAs和8个lncRNAs进行qRT-PCR分析。使用PrimeScrip RT reagent kit with gDNA逆转录试剂盒和TB GreenTMPremix Ex TaqTM II定量试剂盒对目标mRNAs和lncRNAs进行定量,使用M5 miRNA cDNA synthesis Kit加尾法逆转录试剂盒和M5 miRNA qPCR Assay Kit定量试剂盒对目标miRNAs进行定量。相对表达量结果采用2-ΔΔCt法计算得到。mRNAs和lncRNAs定量检测以18S为内参,miRNAs引物信息表见表 1,miRNAs定量检测以U6作为miRNA的内参,上游引物参照定量试剂盒说明书设计,ncRNAs引物信息表见表 2,下游引物为试剂盒通用引物。

表 1 mRNAs实时荧光定量PCR引物序列 Table 1 The primer sequences of mRNAs for qRT-PCR
表 2 ncRNAs实时荧光定量PCR引物序列 Table 2 The primer sequences of ncRNAs for qRT-PCR
2 结果 2.1 幼龄和成年太行山羊附睾头组织全转录组差异表达分析

DESeq2分析结果显示,以幼龄太行山羊附睾头为对照,成年太行山羊附睾头差异表达mRNAs有6 461个,其中2 997个上调、3 464个下调(图 1A);差异表达lncRNAs共有1 147个,其中703个上调、444个下调(图 1B);差异表达miRNAs共有182个,其中81个上调、101个下调(图 1C)。差异表达谱可用于ceRNAs机制分析。

红色代表上调,蓝色代表下调;A、B、C分别代表差异表达mRNAs、lncRNAs和miRNAs;D、E、F分别代表ceRNAs机制中差异表达的mRNAs、lncRNAs和miRNAs Red denote up-regulation, blue denote down-regulation; A, B, C denote the differentially expressed mRNAs, lncRNAs and miRNAs respectively; D, E, F denote the differentially expressed ceRNAs-mRNAs, ceRNAs lncRNAs and ceRNAs miRNAs, respectively 图 1 差异表达RNAs的火山图 Fig. 1 Volcano plot of the differentially expressed RNAs
2.2 mRNAs-lncRNAs-miRNAs差异ceRNAs分析

根据lncRNAs、mRNAs与miRNAs之间的靶向关系和表达量相关性,基于ceRNA-score原理,在筛选出的ceRNAs机制中有3 131个mRNAs,其中1 253个上调,1 878个下调(图 1D);有366个lncRNAs可作为ceRNAs,其中213个上调,153个下调(图 1E);有140个miRNAs可作为ceRNAs,其中48个上调,92个下调(图 1F)。

2.3 差异ceRNAs机制中mRNAs的GO与KEGG分析

对具有ceRNAs机制的幼龄和成年太行山羊附睾头差异mRNAs进行GO富集分析(图 2),结果显示,参与分子功能的mRNAs涉及催化活性、转运活性、分子转导活性等;参与生物学过程的mRNAs主要涉及细胞过程、单有机体过程、生物调控、代谢过程等;参与细胞组分的mRNAs位于细胞器、胞内、胞外区等位置。KEGG分析共挖掘出11个显著富集的通路(图 3),在显著富集通路中有内质网蛋白加工通路、蛋白质输出通路、粘蛋白型O-聚糖生物合成通路、细胞外基质受体相互作用通路,这些通路可参与调节附睾头内精子成熟或免疫过程。

红色代表生物过程,绿色代表细胞组分,蓝色代表分子功能 Red denotes biology process, green denotes cellular component, blue denotes molecular function 图 2 ceRNAs机制中差异表达mRNAs的GO富集分析 Fig. 2 GO analysis of differentially expressed mRNAs in the ceRNAs mechanism
图 3 ceRNAs机制中差异表达mRNAs的KEGG富集分析 Fig. 3 KEGG enrichment analysis of differentially expressed mRNAs in the ceRNAs mechanism
2.4 ceRNAs机制中显著变化的lncRNAs、miRNAs及其靶基因

表 3显示了ceRNAs机制中表达量显著上调和下调的前10位lncRNAs及其预测靶基因。从表 3可知,与幼龄太行山羊附睾头lncRNAs表达量相比,成年附睾头中显著上调前10位lncRNAs的log2(FC)> 11.24。在lncRNAs靶向mRNAs中与免疫相关的基因有DNAJB12、RBM47、NRBF2、DENND1B、USO1、ADAM28、MBOAT7、DUSP6和ACER3。

表 3 ceRNAs机制中上调和下调前10的差异lncRNAs Table 3 Top 10 differentially expressed lncRNAs with up- and down-regulation in the ceRNAs mechanism

表 4显示了ceRNAs机制中表达量显著上调和下调的前10位miRNAs及其预测靶基因。从表 4可知,与幼龄太行山羊附睾头lncRNAs表达量相比,成年附睾头中显著上调前10位miRNAs的log2(FC)>4.54。在miRNAs靶向mRNAs中与免疫相关有SCN8A、CEBPBRHBDL2、UNC93A、TMEM9B。

表 4 ceRNAs机制中上调和下调前10的差异miRNAs Table 4 Top 10 differentially expressed miRNAs with up- and down-regulation in the ceRNAs mechanism
2.5 ceRNAs机制中显著变化mRNAs及免疫相关调控网络构建

表 5显示了ceRNAs机制中表达量显著上调和下调前10位mRNAs。与免疫相关的基因有LY6G5B、脂质运载蛋白9(epididymal-specific lipocalin-9, LCN9)、解整合素金属蛋白酶28(adisintegrinand metalloproteinase, ADAM28)和粘蛋白15(mucin 15, MUC15),与雄性繁殖调控相关的有跨膜附睾蛋白1(transmembrane Epididymal Protein 1, TEDDM1)和ADAM7。上述与免疫和繁殖调控相关的基因在成年太行山羊附睾头表达的FPKM值均在300以上,|log2(FC)|>10且极显著高于幼龄太行山羊(P < 0.01)。

表 5 ceRNAs机制中上调和下调前10的差异mRNAs Table 5 Top10 differentially expressed mRNAs with up- and down-regulation in the ceRNAs mechanism

以表达量显著上调前5位中与免疫相关的基因LY6G5B、LCN9、ADAM28和MUC15为核心,绘制附睾头中调控免疫功能ceRNA网络(图 4)。图 4显示,lncRNA-MSTRG.22929.11、lncRNA-MSTRG.57822.5、lncRNA-MSTRG.26758.1、lncRNA-MSTRG.12113.3、lncRNA-MSTRG.59930.2等lncRNA可靶向miR-495-3p、miR-455-5p和miR-218等miRNAs进而调控LY6G5B、LCN9、ADAM28和MUC15的表达。

红色菱形代表lncRNAs,绿色方形代表mRNAs,蓝色圆形代表miRNAs Red diamonds denote lncRNAs, green squares denote mRNAs, blue circles denote miRNAs 图 4 免疫相关关键的ceRNAs网络 Fig. 4 Network for immune-related key ceRNAs
2.6 全转录组测序结果的qRT-PCR验证

为了验证全转录组测序结果的准确性,对幼龄和成年太行山羊附睾头全转录组差异表达谱中随机选择的8个mRNAs、8个lncRNAs和8个miRNAs进行验证(图 5),结果表明,除chi-miR-320-3p外,其余随机选取的幼龄和成年太行山羊附睾头差异表达mRNAs、lncRNAs和miRNAs的荧光定量结果与全转录组测序结果的差异倍数虽然存在偏差,但是两者之间的上、下调表达趋势一致。

图 5 RNA-Seq差异表达mRNAs和ncRNAs的qRT-PCR验证 Fig. 5 qRT-PCR verification for the differentially expressed mRNAs and ncRNAs results of RNA-Seq
3 讨论

哺乳动物出生后附睾需经过未分化期、分化期和扩张期3个阶段发育后才能发挥正常的生理功能,公鼠分别在出生后21 d完成细胞增殖[21],此时血-附睾屏障也基本形成,以阻止血液中的免疫细胞、免疫蛋白等进入附睾,为附睾细胞差异分化重新构建附睾特异免疫系统提供微环境[22-23]。近年来,研究发现了大量与雄性繁殖相关的lncRNAs和miRNAs等非编码RNAs[24-26],但非编码RNAs对附睾独特免疫系统形成的调控机制仍不明晰,尚需进一步研究探明。本研究采用多组学技术比较分析了幼龄和成年太行山羊附睾头显著差异表达lncRNAs-miRNAs-mRNAs的ceRNAs机制。

本研究基于幼龄和成年太行山羊附睾头多组学分析发现,具有ceRNA调控作用的差异基因显著富集在催化活性、转运活性和分子转导活性等GO分子功能条目上。其中,已有研究证明,表 3中的lncRNA-MSTRG.33179.3靶基因DUSP6[27]、lncRNA-MSTRG.17938.1靶基因USO1[28]、lncRNA-MSTRG.57822.6靶基因PGM1[29]、lncRNA-MSTRG.57822.5靶基因ATP6V0B[30]表 4中的chi-miR-410-3p靶基因UNC93A[31]、chi-miR-376d靶基因RHBDL2[32]、chi-miR-154b-5p靶基因SLC12A5[33]等蛋白均具有转运活性、分子转导或催化功能,这些结果佐证了本研究GO富集分析结果的准确性。

具有ceRNAs调控机制的显著上调的LY6G5B、ADAM28、MUC15和LCN9在成年山羊附睾头中高表达,这些基因均与炎症反应和免疫应答相关[34-41],说明lncRNA作为ceRNA对这些附睾免疫相关基因的表达具有极其复杂的调控过程,同时,这些高度活跃的免疫相关基因蛋白的合成、分泌以及转运需要内质网蛋白加工通路、蛋白质输出通路、粘蛋白型O-聚糖生物合成通路的参与,这支持了本研究KEGG富集分析的结果。

ADAM28是ADAM家族成员,包含ADAM家族共有的金属蛋白酶结构域。金属蛋白酶是调控细胞外基质微环境蛋白降解的主要执行蛋白,同时在细胞-基质的黏连和细胞-细胞的黏连中也具有重要作用[35]。Miyamae等[36]发现, ADAM28定位于附睾上皮细胞,同时,ADAM28可通过结合补体C1q(complement 1q, C1q)诱导支气管上皮细胞死亡和自噬,而C1q表达异常通常见于急性炎症感染。本研究转录组测序结果表明,成年山羊附睾头ADAM28表达量极显著高于幼年,同时发现C1q在附睾头中也有表达,因此推测,山羊性成熟后附睾头ADAM28可能在调控血-附睾屏障上附睾上皮细胞之间的黏连、附睾微环境基质稳态以及附睾上皮细胞死亡和自噬中发挥重要作用,可能是附睾头炎症反应和免疫应答的关键分子。

MUC15是膜型黏蛋白家族成员,对于保护组织黏膜和促进上皮细胞的更新和分化具有重要作用。研究表明,MUC15可以通过激活p-ERK使MAPK信号通路活化[37],参与机体炎症反应和应激反应。也有研究表明,MUC15可以与表皮生长因子受体(epidermal growth factor receptor, EGFR)相互作用调控PI3K/Akt信号通路进而调控基质金属蛋白酶的表达[38],同时,PI3K/Akt信号通路也可以调控白细胞介素-8、白细胞介素-1β、肿瘤坏死因子-α等促炎因子的表达[39-40]。本研究中成年山羊附睾头MUC15表达量显著高于幼年的,但转录组测序结果并未发现附睾头有ERK基因的表达,而EGFR在附睾头中有表达,同时,78个基因富集于PI3K/Akt信号通路,这提示MUC15可能通过与EGFR受体结合来调节PI3K/Akt信号通路进而对附睾头微环境与附睾头炎症反应发挥影响,为本课题组下一步研究提供了思路。

LY6G5B是LY6超家族的成员,在中性粒细胞、子宫内膜细胞和肠上皮细胞中均发现有LY6家族蛋白的表达,其蛋白通过GPI锚定在细胞表面,中性粒细胞Ly6G蛋白可以调节中性粒细胞向炎症部位的迁移[34],介导炎症部位炎症反应。本研究结果表明,附睾头LY6G5B基因在成年附睾头中高度表达,说明LY6G5B在附睾头免疫中发挥着至关重要的作用。LCN9基因编码的蛋白是附睾特异脂质运载蛋白9,是lipocalin家族的成员,在附睾头中高度表达,可运送疏水性小分子到特定细胞,参与调节免疫反应与炎症应答[41],但其在山羊附睾头的功能及作用机制仍需进一步探明。

图 5绘制的免疫相关关键ceRNAs网络中,miR-495-3p、miR-455-5p和miR-218已被证实与免疫相关。miR-495-3p是与炎性反应相关的miRNA,其靶向Toll样受体基因(toll like receptor 4, TLR4),调控NF-κB通路抑制炎性因子的产生[42]。miR-455-5p直接靶标髓样分化因子(myeloid differentiation factor 88, MyD88)和Rel,在多发性硬化症中发挥抗炎作用[43]。miR-218通过下调2,3双加氧酶(indoleamine 2, 3-dioxygenase 1, IDO1)抑制宫颈癌细胞的免疫逃逸[44]。本研究的附睾头转录组测序数据中也有TLR4、MyD88、RelIDO1基因,另外miR-495-3p、miR-455-5p和miR-218靶向LY6G5B、LCN9、ADAM28和MUC15免疫相关基因,说明miR-495-3p、miR-455-5p和miR-218通过这些靶标基因发挥免疫调节作用。本试验筛选出的ceRNAs均为生物信息学分析预测结果,还需进一步进行试验验证。

近年来, 研究发现了lncSmad3、lncEGFR和LncLsm3b等大量与免疫相关的lncRNAs[45],这些研究表明,lncRNAs对炎症反应和免疫应答具有重要作用,本研究分析结果也说明,lncRNAs可做为附睾头免疫应答调控分子。由于lncRNAs可通过介导表观遗传、ceRNAs机制、miRNAs前体等多种方式调控机体生命活动[46],ceRNA机制只是其中一种,因此,本研究筛选所得具有ceRNAs机制中lncRNAs的其他作用有待进一步挖掘。另外,本试验从幼龄和成年太行山羊附睾头差异ceRNAs机制表达谱中除筛选得到免疫相关关键ceRNAs调控网络之外,该机制中其他基因还参与能量代谢、抗氧化系统、蛋白合成等相关生物学过程,推测这可能与附睾头发育及附睾精子成熟等有关,具体调控机制还需进一步试验验证。

4 结论

本研究利用高通量测序技术及生物信息学分析获得了3 131个幼龄和成年太行山羊的附睾头差异表达ceRNAs,具有ceRNAs机制的基因显著富集通路与免疫功能相关。成年附睾头高表达的免疫相关基因LY6G5B、LCN9、ADAM28和MUC15具有复杂ceRNAs调控机制,这为进一步研究非编码RNAs对太行山羊附睾头免疫调控作用提供了参考依据。

参考文献
[1] JIANG F X, TEMPLE-SMITH P, WREFORD N G. Postnatal differentiation and development of the rat epididymis: a stereological study[J]. Anat Rec, 1994, 238(2): 191–198. DOI: 10.1002/ar.1092380205
[2] DACHEUX J L, BELLEANNÉE C, JONES R, et al. Mammalian epididymal proteome[J]. Mol Cell Endocrinol, 2009, 306(1-2): 45–50. DOI: 10.1016/j.mce.2009.03.007
[3] AGARWAL A, HOFFER A P. Ultrastructural studies on the development of the blood-epididymis barrier in immature rats[J]. J Androl, 1989, 10(6): 425–431. DOI: 10.1002/j.1939-4640.1989.tb00132.x
[4] GATTI J L, CASTELLA S, DACHEUX F, et al. Post-testicular sperm environment and fertility[J]. Anim Reprod Sci, 2004(82-83): 321–339.
[5] FRENETTE G, GIROUARD J, SULLIVAN R. Comparison between epididymosomes collected in the intraluminal compartment of the bovine caput and cauda epididymidis[J]. Biol Reprod, 2006, 75(6): 885–890. DOI: 10.1095/biolreprod.106.054692
[6] CYR D G, ROBAIRE B, HERMO L. Structure and turnover of junctional complexes between principal cells of the rat epididymis[J]. Microsc Res Technol, 1995, 30(1): 54–66. DOI: 10.1002/jemt.1070300105
[7] DONG R L, LIU J W, SUN W, et al. Comprehensive analysis of aberrantly expressed profiles of lncRNAs and miRNAs with associated ceRNA network in lung adenocarcinoma and lung squamous cell carcinoma[J]. Pathol Oncol Res, 2020, 26(3): 1935–1945. DOI: 10.1007/s12253-019-00780-4
[8] FENG C, LI Y X, LIN Y, et al. CircRNA-associated ceRNA network reveals ErbB and Hippo signaling pathways in hypopharyngeal cancer[J]. Int J Mol Med, 2019, 43(1): 127–142.
[9] STRANIERO L, RIMOLDI V, SAMARANI M, et al. The GBAP1 pseudogene acts as a ceRNA for the glucocerebrosidase gene GBA by sponging miR-22-3p[J]. Sci Rep, 2017, 7(1): 12702. DOI: 10.1038/s41598-017-12973-5
[10] LYU K X, LI Y, XU Y, et al. Using RNA sequencing to identify a putative lncRNA-associated ceRNA network in laryngeal squamous cell carcinoma[J]. RNA Biol, 2020, 17(7): 977–989. DOI: 10.1080/15476286.2020.1741282
[11] MENG X N, PENG H R, DING Y Z, et al. A transcriptomic regulatory network among miRNAs, piRNAs, circRNAs, lncRNAs and mRNAs regulates microcystin-leucine arginine (MC-LR)-induced male reproductive toxicity[J]. Sci Total Environ, 2019, 667: 563–577. DOI: 10.1016/j.scitotenv.2019.02.393
[12] HU J, ZHANG L, MEI Z, et al. Interaction of E3 ubiquitin ligase MARCH7 with long noncoding RNA MALAT1 and autophagy-related protein ATG7 promotes autophagy and invasion in ovarian cancer[J]. Cell Physiol Biochem, 2018, 47(2): 654–666. DOI: 10.1159/000490020
[13] LIU S P, QIU J J, TANG X Y, et al. LncRNA-H19 regulates cell proliferation and invasion of ectopic endometrium by targeting ITGB3 via modulating miR-124-3p[J]. Exp Cell Res, 2019, 381(2): 215–222. DOI: 10.1016/j.yexcr.2019.05.010
[14] 周璇, 朱永通, 褚庆军, 等. 长链非编码RNA作为竞争性内源RNA在非梗阻性无精子症中的作用机制探讨[J]. 中华医学杂志, 2019, 99(35): 2761–2767.
ZHOU X, ZHU Y T, CHU Q J, et al. Effects and mechanism of lncRNA serving as ceRNA in non-obstructive azoospermia[J]. National Medical Journal of China, 2019, 99(35): 2761–2767. DOI: 10.3760/cma.j.issn.0376-2491.2019.35.009 (in Chinese)
[15] 杜海燕. 山羊附睾头β防御素家族表达特点及gBD 124功能分析[D]. 晋中: 山西农业大学, 2018.
DU H Y.The expression characteristics of β-defensin family from caprine epididymis and function analysis of gBD 124[D]. Jinzhong: Shanxi Agricultural University, 2018.(in Chinese)
[16] 任有蛇, 郭丽娜, 张春香, 等. 山羊Lcn 5的表达特点及其在繁殖器官中定位[J]. 畜牧兽医学报, 2015, 46(5): 711–718.
REN Y S, GUO L N, ZHANG C X, et al. Expression characteristics of Lcn 5 and its localization in reproduction organ of bucks[J]. Acta Veterinaria et Zootechnica Sinica, 2015, 46(5): 711–718. (in Chinese)
[17] 张春香, 张国林, 郭丽娜, 等. 基于高通量转录组测序的山羊睾丸和附睾头差异表达基因分析[J]. 畜牧兽医学报, 2014, 45(3): 391–401.
ZHANG C X, ZHANG G L, GUO L N, et al. Study on differentially expressed genes between caprine testis and epididymis caput based on transcriptomes with high-throughput RNA-seq technology[J]. Acta Veterinaria et Zootechnica Sinica, 2014, 45(3): 391–401. (in Chinese)
[18] WANG D Y, CHEN Z J, ZHUANG X N, et al. Identification of circRNA-associated-ceRNA networks involved in milk fat metabolism under heat stress[J]. Int J Mol Sci, 2020, 21(11): 4162. DOI: 10.3390/ijms21114162
[19] WANG X G, YANG C H, GUO F, et al. Integrated analysis of mRNAs and long noncoding RNAs in the semen from Holstein bulls with high and low sperm motility[J]. Sci Rep, 2019, 9(1): 2092. DOI: 10.1038/s41598-018-38462-x
[20] LIU Y F, SUN Y Y, LI Y L, et al. Analyses of long non-coding RNA and mRNA profiling using RNA sequencing in chicken testis with extreme sperm motility[J]. Sci Rep, 2017, 7(1): 9055. DOI: 10.1038/s41598-017-08738-9
[21] ROBAIRE B, HINTON B T.Chapter 17-the epididymis[M]//PLANT T M, ZELEZNIK A J.Knobil and Neill's Physiology of Reproduction.San Diego: Academic Press, 2015: 691-771.
[22] VOISIN A, WHITFIELD M, DAMON-SOUBEYR-AND C, et al. Comprehensive overview of murine epididymal mononuclear phagocytes and lymphocytes: unexpected populations arise[J]. J Reprod Immunol, 2018, 126: 11–17. DOI: 10.1016/j.jri.2018.01.003
[23] WIJAYARATHNA R, HEDGER M P. Activins, follistatin and immunoregulation in the epididymis[J]. Andrology, 2019, 7(5): 703–711.
[24] ZHANG F J, ZHANG Y, LV X L, et al. Evolution of an x-linked miRNA family predominantly expressed in mammalian male germ cells[J]. Mol Biol Evol, 2019, 36(4): 663–678. DOI: 10.1093/molbev/msz001
[25] YANG H, WANG F, LI F Z, et al. Comprehensive analysis of long noncoding RNA and mRNA expression patterns in sheep testicular maturation[J]. Biol Reprod, 2018, 99(3): 650–661. DOI: 10.1093/biolre/ioy088
[26] LI Y, WANG H Y, QIN Y M, et al. Deep sequencing reveals microRNA signature is altered in the rat epididymis following bilateral castration[J]. Genes Genom, 2019, 41(7): 757–766. DOI: 10.1007/s13258-019-00803-z
[27] LEE B, SAHOO A, SAWADA J, et al. MicroRNA-211 modulates the DUSP6-ERK5 signaling axis to promote BRAFV600E-driven melanoma growth in vivo and BRAF/MEK inhibitor resistance[J]. J Invest Dermatol, 2020. DOI: 10.1016/j.jid.2020.06.038
[28] HEO Y, YOON H J, KO H, et al. Crystal structures of Uso1 membrane tether reveal an alternative conformation in the globular head domain[J]. Sci Rep, 2020, 10(1): 9544. DOI: 10.1038/s41598-020-66480-1
[29] BACKE P H, LAERDAHL J K, KITTELSEN L S, et al. Structural basis for substrate and product recognition in human phosphoglucomutase-1(PGM1) isoform 2, a member of the α-D-phosphohexomutase superfamily[J]. Sci Rep, 2020, 10(1): 5656. DOI: 10.1038/s41598-020-62548-0
[30] MEO-EVOLI N, ALMACELLAS E, MASSUCCI F A, et al. V-ATPase: a master effector of E2F1-mediated lysosomal trafficking, mTORC1 activation and autophagy[J]. Oncotarget, 2015, 6(29): 28057–28070. DOI: 10.18632/oncotarget.4812
[31] KIM T H, KIM D, LEE Y, et al. Expression of UNC93A induced by CpG-DNA-liposome complex in mice[J]. J Korean Soc Appl Biol Chem, 2014, 57(3): 281–287. DOI: 10.1007/s13765-014-4050-z
[32] NOY P J, SWAIN R K, KHAN K, et al. Sprouting angiogenesis is regulated by shedding of the C-type lectin family 14, member A (CLEC14A) ectodomain, catalyzed by rhomboid-like 2 protein (RHBDL2)[J]. FASEB J, 2016, 30(6): 2311–2323. DOI: 10.1096/fj.201500122R
[33] FUKUDA A, WATANABE M. Pathogenic potential of human SLC12A5 variants causing KCC2 dysfunction[J]. Brain Res, 2019, 1710: 1–7. DOI: 10.1016/j.brainres.2018.12.025
[34] LEE P Y, WANG J X, PARISINI E P, et al. Ly6 family proteins in neutrophil biology[J]. J Leukoc Biol, 2013, 94(4): 585–594. DOI: 10.1189/jlb.0113014
[35] NARITA D, SECLAMAN E, URSONIU S, et al. Increased expression of ADAM12 and ADAM17 genes in laser-capture microdissected breast cancers and correlations with clinical and pathological characteristics[J]. Acta Histochem, 2012, 114(2): 131–139. DOI: 10.1016/j.acthis.2011.03.009
[36] MIYAMAE Y, MOCHIZUKI S, SHIMODA M, et al. ADAM28 is expressed by epithelial cells in human normal tissues and protects from C1q-induced cell death[J]. FEBS J, 2016, 283(9): 1574–1594. DOI: 10.1111/febs.13693
[37] 刘丽云, 龚健, 徐晋珩, 等. 甲状腺乳头状癌BRAF基因突变及BRAF、MUC15、RKIP蛋白表达的临床意义[J]. 中国临床研究, 2017, 30(11): 1450–1454, 1458.
LIU L Y, GONG J, XU J H, et al. Clinical significance of BRAF gene mutation and BRAF, MUC15 and RKIP proteins expressions in papillary thyroid carcinoma[J]. Chinese Journal of Clinical Research, 2017, 30(11): 1450–1454, 1458. (in Chinese)
[38] WANG R Y, CHEN L, CHEN H Y, et al. MUC15 inhibits dimerization of EGFR and PI3K-AKT signaling and is associated with aggressive hepatocellular carcinomas in patients[J]. Gastroenterology, 2013, 145(6): 1436–1448. DOI: 10.1053/j.gastro.2013.08.009
[39] CHANG M C, LEE J J, CHEN Y J, et al. Lysophosphatidylcholine induces cytotoxicity/apoptosis and IL-8 production of human endothelial cells: related mechanisms[J]. Oncotarget, 2017, 8(63): 106177–106189. DOI: 10.18632/oncotarget.22425
[40] ZHENG W H, FENG Z H, LOU Y T, et al. Silibinin protects against osteoarthritis through inhibiting the inflammatory response and cartilage matrix degradation in vitro and in vivo[J]. Oncotarget, 2017, 8(59): 99649–99665. DOI: 10.18632/oncotarget.20587
[41] LI X Q, ZHAN X N, LIU S G, et al. Cloning and primary characterizations of rLcn9, a new member of epididymal lipocalins in rat[J]. Acta Biochim Biophys Sin (Shanghai), 2012, 44(10): 876–885. DOI: 10.1093/abbs/gms072
[42] 张帆, 余国玺, 王力锋. miR-495-3p对LPS刺激人牙周膜细胞增殖及炎性因子表达的影响[J]. 口腔生物医学, 2017, 8(2): 86–89.
ZHANG F, YU G X, WANG L F. Effect of miR-495-3p on the proliferation of hPDLC and the production of inflammatory cytokines induced by LPS[J]. Oral Biomedicine, 2017, 8(2): 86–89. (in Chinese)
[43] TORABI S, TAMADDON M, ASADOLAHI M, et al. miR-455-5p downregulation promotes inflammation pathways in the relapse phase of relapsing-remitting multiple sclerosis disease[J]. Immunogenetics, 2019, 71(2): 87–95. DOI: 10.1007/s00251-018-1087-x
[44] ZHU L, TU H D, LIANG Y M, et al. MiR-218 produces anti-tumor effects on cervical cancer cells in vitro[J]. World J Surg Oncol, 2018, 16(1): 204. DOI: 10.1186/s12957-018-1506-3
[45] 雷艳婷, 李兆瀛, 孙博, 等. LncRNA和RBP在免疫调节中的研究进展[J]. 中国免疫学杂志, 2020, 36(15): 1896–1901.
LEI Y T, LI Z Y, SUN B, et al. Research progress of LncRNA and RBP in immune regulation[J]. Chinese Journal of Immunology, 2020, 36(15): 1896–1901. DOI: 10.3969/j.issn.1000-484X.2020.15.022 (in Chinese)
[46] LAURENT G S, WAHLESTEDT C, KAPRANOV P. The landscape of long noncoding RNA classification[J]. Trends Genet, 2015, 31(5): 239–251. DOI: 10.1016/j.tig.2015.03.007