幼龄和成年太行山羊附睾头差异竞争性内源RNAs机制分析

郭相前，董复成，黄鑫芸，刘文忠，乔利英，张春香，任有蛇

(山西农业大学动物科学学院，太谷 030801)

哺乳动物在出生时附睾处在未分化期，需要再经过一段以细胞继续增殖为主要特征的未分化期[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,lncRNAATG7)与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,lncRNAANXA2P3)可以通过竞争性结合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 TaqTMII(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 TaqTMII定量试剂盒对目标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机制分析。

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)。

红色代表上调，蓝色代表下调；A、B、C分别代表差异表达mRNAs、lncRNAs和miRNAs；D、E、F分别代表ceRNAs机制中差异表达的mRNAs、lncRNAs和miRNAsRed 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.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的差异lncRNAsTable 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、CEBPB、RHBDL2、UNC93A、TMEM9B。

表4 ceRNAs机制中上调和下调前10的差异miRNAsTable 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的差异mRNAsTable 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，蓝色圆形代表miRNAsRed 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、Rel和IDO1基因，另外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对太行山羊附睾头免疫调控作用提供了参考依据。