小鼠中年期和老年期睾丸组织转录组分析

张 涛，杨理凯，路宏朝，王 令，刘 欢，左甜甜

(陕西理工大学 生物科学与工程学院，陕西汉中 723001)

老龄化问题已成为世界各国面临的共同课题，中国老龄化速度不断加快。衰老涉及内分泌系统多个方面的改变，包括肾上腺、甲状腺、垂体生长激素、性腺、骨矿物质代谢和糖代谢等[1-2]。Holladne等在1958年发现老年男子的睾丸精索静脉睾酮水平低于青年男子[3]，是衰老降低睾酮产生的有力证据。近年来，学者把由老龄化所导致的血液中雄激素低下为生物化学基础的症候群称之为“部分雄激素缺乏综合征(partial androgen deficiency of agingmale,PADAM)”，临床表现为失眠、健忘、性功能障碍等多种老年性疾病[4-5]。睾丸是精子产生的场所，并分泌雄激素。精子在曲精小管中从精原干细胞到精母细胞，然后通过减数分裂形成精子，该过程也涉及睾丸中支持细胞及间质细胞等。睾丸支持细胞分布在曲精小管，负责生殖细胞的营养与保护，而间质细胞分布于睾丸间质中[6]，在脑垂体分泌的黄体生成素作用下，生成和分泌雄激素——睾酮，促进性腺发育、精子发生，维持男性性征[7]，男性40岁后随着年龄的增长，下丘脑-垂体-睾丸轴功能减退，睾丸功能降低。

睾丸作为雄性的主要标志器官，其生长发育是一个复杂而精密的生理过程，随着年龄老化而出现睾丸内分泌变化与功能衰退，多数学者认为是由于衰老对下丘脑-垂体-睾丸轴功能产生影响而导致，目前已开展大量人与动物睾丸发育[8]、损伤与保护[9]及组织中基因表达差异[10]等研究，但是相关的生理机制和分子机理仍然没有阐述清楚。测序技术的出现为人类从分子水平解析生命的奥秘提供新的途径，转录组测序与分析技术的应用已经有20多年，广泛应用于生物学的研究，且对小鼠不同发育阶段睾丸组织中不同类型细胞的基因表达进行研究[11-13]，但多数研究集中在胚胎期、幼年、青年和成年期，而关于中年期和老年期的睾丸组织的转录组研究未见报道。本研究选取出生后180日龄和360日龄小鼠，分别代表中年期和老年期，采用Illumina Hiseq测序平台对睾丸组织的基因表达情况进行比较分析，为探讨睾丸功能退化的分子机制提供基础数据，以期为老龄化过程中生物学与医学等某些方面科学问题的解决提供新思路。

1 材料与方法

1.1 试验动物

C57BL小鼠购自第四军医大学实验动物中心，实验室长期饲养，选取中年期(180日龄)和老年期(360日龄)的睾丸组织进行研究。

1.2 试验方法

1.2.1 总RNA提取 每只小鼠取 1 个睾丸，按照Trizol Reagent 试剂盒的操作流程提取总RNA，测定质量浓度后按照中年期和老年期 2 个组分别将 3 个样本等量混合构建 cDNA 文库。

1.2.2 测序文库的构建 基于 mRNA 特有的Poly(A)结构进行 mRNA 纯化。通过离子打断的方式将 mRNA 打断为 200～300 bp。采用随机引物和逆转录酶合成 cDNA 第 1 条链，合成第 2 条 cDNA 链时将 dTTP 替换成 dUTP，以提高结果的准确性。cDNA 文库构建完成后，采用 PCR 扩增进行文库片段富集，建立测序文库，文库大小在 300～400 bp，最后通过 Agilent 2100Bioanalyzer 对文库进行质检。

1.2.3 Illumina测序 选择上海派森诺生物科技有限公司 Illumina Hiseq 测序平台对文库进行双末端(Paired End，PE)测序 。

1.2.4 数据处理 采用 Cutadapt(Version 1.2.1)除去3′端的接头，数据平均含量按照 Q20 标准控制。进行比对的参考基因组数据来源为 Ensembl 数据库(http:∥www.ensembl.org/)，参考的基因组信息为Mus_musculus.GRCm38.dna.primary_assembly.fa，数据库版本为 Ensembl 86.38。通过 Bowtie2 建立参考基因组索引，然后使用 Tophat2(http:∥tophat.cbcb.umd.edu/)将过滤后的 Reads 比对到参考基因组上。

1.2.5 差异基因筛选 使用 HTSeq 0.6.1p2(http:∥wwwhuber.embl.de/users/anders/HTSeq)统计比对到每一个基因上ReadCount 值，作为基因的原始表达量。采用 RPKM(Reads Per Kilo bases per Million reads)[14]对表达量进行标准化(Normalization)。采用 DESeq(version 1.18.0)[15]对基因表达进行差异分析，筛选差异表达基因条件为：表达倍数差异 |fold change| > 2，显著性P<0.05。采用 R 语言 ggplots2 软件包绘制差异表达基因的火山图和 MA 图。

1.2.6 GO与KEGG分析 将所有检测到的基因映射到GO(Gene Ontology)数据库每一个 Term，计算每个 Term 的差异基因数目，以整个基因组为背景，采用超几何分布计算差异基因显著富集的 Term 及相应功能；在 KEGG(kyoto encyclopedia of genes and genomes)数据库比对差异表达基因，统计各个 KEGG Pathway 各层级的差异表达基因数目，确定差异表达基因主要参与的代谢途径和信号通路。

2 结果与分析

2.1 数据质量检测

由表1可知，2个样品分别获得25 479 636和 29 158 408 条 reads，总碱基数分别为3 847 425 036 bp和 4 402 919 608 bp，其中识别准确率在 99.9% 以上的碱基总数为 3 439 210 859 bp和 3 967 709 437 bp；Q20 值均在 95% 以上，Q30 在 89% 以上，结果表明测序质量好，保证研究的可靠性。为了排除接头和低质量 Reads 对信息分析的干扰，试验进行数据的过滤处理(表2)，试验两个样品分别检测到 25 187 980 个和28 840 576个有效的Reads，比对到参考基因组上的分别有 18 861 721 个和 22 355 657 个，占总读数的 74.88% 和 77.51%。比对结果显示比对效率较好，能比对到功能基因的比例分别为 76.51% 和 79.96%，其中比对到外显子区的分别为 95.54% 和 96.06%，结果如表3所示，说明测序的数据丰富，有效性好。

2.2 差异表达基因的筛选

为系统了解小鼠中年期及老年期睾丸组织基因表达情况，笔者计算检测到的每一个基因的 RPKM 值，如表4所示。转录组研究领域的主流杂志提出，在有参转录组当中，一般认为 RPKM>1的基因是表达的，RPKM 能够很好地反应基因的表达水平。在中年期小鼠的睾丸组织可检测到 17 770个基因的表达，其中 RPKM 值小于 1 的基因占 62.15%，大于 100 的基因占 0.80%；在老年期小鼠的睾丸组织中，检测到 18 073 个基因的表达，其中 RPKM 值小于1的基因占61.45%，大于 100 的基因占 0.74%。火山图可直观反映基因的分布情况，MA图常用来评估文库Normalization的质量，结果表明(图1、图2)中年期和老年期睾丸组织转录组表达水平存在较大差异。以中年期为对照，老年期表现下调基因 122 个，上调基因 185 个(图3) 其中上调基因中可注释的 181 个，下调基因可注释的 104 个。所有差异表达基因中，RPKM>1 的基因共有 57 个，下调 18 个，上调 39 个。

表1 转录组测序数据评价分析Table 1 Evaluation of transcriptome sequencing data

表2 转录组测序数据比对结果基本统计Table 2 Basic statistics of transcriptome sequencing data

表3 转录组数据比对区域分布统计Table 3 Transcriptome data comparison regional distribution statistics

2.3 GO分类

利用 Blast2GO 软件对 307 个差异表达基因进行分析表明，285 个基因可以分为生物过程、细胞组分和分子功能 3 大类别 25 小类，共涉及生物过程的 1 625 个条目，细胞组分的 244 个条目，分子功能的 333 个条目。 差异表达基因在细胞组分分类下主要涉及膜及其组成相关的各种功能，其中最多的是膜的整体组成(图4)。在生物过程大类下基因涉及的生物学功能最复杂，涉及基因较多的是转运、RNA聚合酶Ⅱ启动子转录调控、氧化还原过程、转录调控与DNA模板化和老龄化。在分子功能分类中，氧化还原酶活性所占比例最多，其次是转运蛋白活性、丝氨酸型内肽酶抑制剂活性、转运体活性和裂解酶活性。此外，从图4也可以看出，与生殖发育和衰老相关的生物学过程有细胞增殖调节、免疫系统过程、骨化、老化、丝氨酸型内肽酶抑制剂活性、氧化还原酶活性等。

表4 不同生长时期睾丸组织RPKM值Table 4 FPKM value of testis tissue in different growth stages

图1 差异表达基因火山图Fig.1 Valcano chart of differentially expressed genes

图2 差异表达基因MA图Fig.2 MA chart of differentially expressed genes

图3 中年期与老年期睾丸组织差异表达基因统计Fig.3 Statistics of differentially expressed genes in middle-and old-aged mouse

B1.炎症反应正向调节 Positive regulation of inflammatory response； B2.氧化还原过程 Oxidation-reduction process；B3.骨化 Ossification；B4.老化 Aging；B5.运输 Transport；B6.RNA聚合酶Ⅱ启动子转录的正调节 Positive regulation of transcription from RNA polymerase Ⅱ promoter；B7.转录正调节，DNA模板化 Positive regulation of transcription, DNA-templated；B8.膜转运 Transmembrane transport；B9.免疫系统过程 Immune system process；B10.调节细胞增殖 Regulation of cell proliferation；C1.细胞外外泌体 Extracellular exosome；C2.细胞外区域 Extracellular region；C3.内质网 Endoplasmic reticulum；C4.细胞外空间 Extracellular space；C5.内体 Endosome；C6.质膜的整体成分 Integral component of plasma membrane；C7.膜的整体组成部分 Integral component of membrane；C8.质膜 Plasma membrane；C9.膜 Membrane；M1.丝氨酸型内肽酶抑制剂活性 Serine-type endopeptidase inhibitor activity；M2.氧化还原酶活性 Oxidoreductase activity；M3.转运体活性 Symporter activity；M4.裂解酶活性 Lyase activity；M5.运输活动 Transporter activity；M6.肝素结合 Heparin binding

2.4 KEGG富集分析

通过KEGG 数据库进行差异基因比对与富集，结果发现 301 个基因归纳到 5 个大类的 18 个小类中，共涉及 139 个信号通路，如图5所示。结果表明，这些差异表达基因主要参与的前 5 个代谢通路为信号转导、脂代谢、碳水化合物代谢、内分泌系统和神经系统。KEGG 归类分析后，筛选出 RPKM>1 的 57 个基因，并统计分析这些基因涉及的信号通路，最终检测到 20 个基因涉及 35 个信号通路(表5)。小鼠中年期睾丸组织有 6 个基因显著上调，涉及PPAR 信号通路、脂类调控、补体系统、氮代谢、嗅觉转导、吞噬体、肠免疫网络、抗原加工与展示和细胞粘附分子等 9个信号通路。而在老年期检测到 15 个基因上调，涉及20个信号通路，主要包括花生四烯酸代谢、卵巢类固醇生成、类固醇激素生成、醛固酮的生成与分泌、肾素血管生成素系统、血管平滑肌、肾素分泌、氨基酸代谢、催产素信号通路、促性腺激素释放激素、TNF信号通路、甘油三酯代谢、补体系统、胰岛素、溶酶体和鞘脂类代谢与信号通路。

A.异生素的生物降解和代谢 Xenobiotics biodegradation and metabolism;B.概述 Overview; C.核苷酸代谢 Nucleotide metabolism;D.其他氨基酸的代谢 Metabolism of other amio acids;E.油脂代谢 Lipid metabolism;F.聚糖生物合成和代谢 Glycan biosynthesis and metabolism;G.碳水化合物代谢 Carbohydrate metabolosm;H.氨基酸代谢 Amio acidsmetabolism;I.感官系统 Sensory system;J.神经系统 Nervous system;K.免疫系统 Immune system;L.内分泌系统 Endocrine system;M.循环系统 Circulatory system;N.信号分子和相互作用 Signaling molecules and interaction;O.信号传导 Signal transduction;P.运输和分解代谢 Transport and catabolism;Q.细胞群落 Cellular community;R. 翻译 Translation

3 讨论与结论

3.1 转录组测序可靠性分析

转录组测序技术为基因表达及调控研究和生物重要性状候选基因的鉴定提供新的途径，尤其在不同物种、不同发育阶段及不同组织之间挖掘差异表达基因具有重要作用[16]。本研究通过 Illumina Hiseq 测序平台分析检测小鼠中年期和老年期睾丸组织的 cDNA 文库，获得小鼠睾丸组织在中、老年期大量的基因表达信息与变化差异，结果表明，选择的试验材料测序质量好，各样品的 Q20 均在 95% 以上，Q30 均在 89%以上；在小鼠中年期和老年期睾丸组织分别检测到 25 187 980 个和 28 840 576 个可用的 Reads，比对效率分别达到 74.88% 和 77.51%；测序质量高于前期小鼠睾丸组织 6 日龄 (52.87%) 4 周龄(56.70%)和 10 周龄 (62.35%) 的转录组测序结果[13]，也比其他物种如鸭卵巢组织(70.92%和72.70%)[16]和绵羊(68.05%和71.82%)[17]等组织的高。

3.2 中年期和老年期睾丸组织基因差异表达分析

自从转录组分析技术出现以来，科学家就将该技术应用到小鼠睾丸的基因表达研究，如抑制消减杂交(SSH)、表达序列标签(expressed sequence tag,EST)、cDNA微阵列和基因表达系列分析(serial analysis of gene expression,SAGE)等技术。通过这些技术构建小鼠的睾丸体细胞基因表达谱，鉴定出大量小鼠睾丸组织唯一表达标签，发现一些在睾丸生精细胞、支持细胞和间质细胞特异表达的基因，并在小鼠不同发育阶段检测到许多与精子生成有关的差异表达基因[11-12]。高通量转录组测序技术比传统转录组分析技术更能敏感地检测样本间的基因差异表达。龚未等[13]采用超高通量新一代 Technologies SOLiD 测序平台分析 6 日龄、4 周龄和 10 周龄小鼠睾丸组织的基因表达数据,鉴定出 18 837 个基因，同时鉴定不同发育阶段差异表达的基因。但目前利用高通量转录组测序技术分析小鼠中年期和老年期睾丸组织的基因差异表达未见报道。

本研究对小鼠中年期和老年期睾丸组织cDNA文库进行分析，共检测到 307 个表达差异的基因，可注释的 271 个，上调 173 个，下调 98 个，其中 264 个基因获得 GO 功能注释，187 个基因清楚涉及 KEGG 代谢通路。由于小鼠作为模式生物研究的比较深入，比对数据库中的参考基因比较全面，所以本研究中检测到差异表达基因功能在 GO、KEGG 中被注释的比例比其他物种覆盖度高。但是，也有一些基因没有被注释到，也没有归入到相应的 GO 和 KEGG 功能分类，分析原因可能是目前小鼠睾丸组织的研究多数集中于胚胎发育和青年期，而中年期和老年期的研究相对较少，而那些到中老年期后表达的一些特异基因没有被发现，也有可能是转录组测序及拼接序列较短，获得注释信息的可能性就降低。

表5 差异表达基因涉及的信号通路Table 5 Differentially expressed genes involved in signaling pathways

3.3 差异表达基因的功能分析

GO 归类分发现在细胞组分大类中涉及的基因比较集中，主要是膜相关类别，说明在小鼠中老年期和老年期睾丸生理活动中的差异表达基因主要与膜的生物学功能相关，同时也发现一些基因涉及细胞外的活动，该结果与龚未等[13]通过转录组分析技术研究小鼠出生后的基因差异表达研究是一致，说明细胞外活动相关基因贯穿小鼠睾丸整个生理生化过程。而在生物过程中差异表达基因涉及的范围比较分散和多样，各组分中相关基因数量较少，且差异不明显，这与其前期关于小鼠睾丸组织和其他物种的研究基本一致，推断生物过程在睾丸组织中是最基本的代谢途径。对分子功能大类的分析发现，氧化还原酶活性所占比例最多，可能与老龄化密切相关，随着年龄增大，需要大量的氧化还原酶发挥抗衰老的作用。此外，也检测到明显与生殖发育、性功能及衰老相关的差异表达基因。本研究结果证实膜相关、细胞外、转运和氧化还原等类别在中老年睾丸组织中存在较大差别的基因表达，推测膜相关、细胞外、转运和氧化还原在睾丸功能退化的生理活动中可能具有重要的作用。

遗传学已经阐明，表型的产生有主效基因作用、微效多基因假说及其与环境互作效应等机制。现代生物学研究证明生物体中不同类型的生物学功能由多个基因相互协调进行调控，形成一个精密的调控网络，通过信号通路富集分析有助于系统而全面地解析基因功能。转录组分析可以检测到不同物种或者不同时间与空间的遗传表达差异，实现差异表达基因的信号通路富集分析。本研究检测到与小鼠中年期和老年期睾丸组织转录组表达差异相关的 139 个信号通路，其中信号转导、脂代谢、碳水化合物代谢、内分泌系统和神经系统等 5 个代谢通路涉及的差异表达基因较多，同时与癌症和免疫学疾病相关的信号通路富集一些基因。有研究报道，在有参照转录组分析当中，一般认为 RPKM>1 的基因是表达的，本研究着重将差异表达基因中 RPKM>1 的进行分析，发现中年期和老年期小鼠睾丸组织在脂代谢、GnRH 信号通路、固醇类激素生成、TNF 信号通路、补体系统、胰岛素和溶酶体相关等信号通路存在较大差异。本研究中涉及脂肪合成与代谢的有脂肪酸结合蛋白4( FABP4)、内皮脂肪酶(Lipg)、蜡样脂褐质沉着症神经元5( Cln5)和鞘磷脂磷酸二酯酶1酸性溶酶体( Smpd1)。 FABP4 在脂肪、卵巢和睾丸组织中广泛分布，能结合多种疏水性化合物，参与脂肪酸的合成[18]，该基因的缺失可增加脂肪酸浓度，促使脂肪分解[19]。本研究中发现年小鼠睾丸组织中 FABP4 的表达水平高于老年小鼠，说明睾丸组织随着年龄的增加脂肪分解能力提高，贮存脂肪能力减少，为激素合成与精子生成提供的能量就不足。Lipg在内皮细胞合成，冠状动脉、肝脏、肺、睾丸、肾脏和卵巢多种组织中表达，主要功能是参与脂蛋白代谢，其磷脂酶活性强于甘油三酯酶，对富含磷脂的高密度脂蛋白的水解作用更强，在动物试验和人体内均发现内皮脂酶的水平与高密度脂蛋白(HDL-C)的水平呈负相关[20]。也有研究发现Lipg与炎症反应也有密切的联系，可促进细胞间黏附因子1的表达及单核细胞对血管内皮粘的附作用，这都与老龄化有密切的关系。本研究在小鼠中年期和老年期均检测到该基因表达，且老年小鼠睾丸组织中Lipg基因表达水平显著高于中年小鼠，结果与其组织器官的研究基本一致。

睾丸雄激素的合成由下丘脑分泌的 GnRH 调节，其通过促进垂体分泌促性腺物质影响性激素生成而影响雄性动物生殖活动[21]。本研究中 GnRH 信号通路和 TNF 信号通路涉及基质金属蛋白酶(MMP-14)，睾丸组织中老年鼠要显著高于中年鼠，有研究报道 MMP-14 异常表达与多种恶性肿瘤密切相关[22-23]，说明老龄化可增加患恶性肿瘤的风险。前列腺素D合成酶(Prostaglandin D synthase，PGDS)催化不稳定的中间产物前列腺素 H2(PGH2) 发生异构化生成，其分为脑型脂质运载蛋白 PGDS(lipocalin-PGDS，LPGDS) 和生血型 PGDS(hematopoietic PGDS，hPGDS)[24]。LPGDS在中枢神经系统、雄性生殖器官和心脏等组织分布，具有转运亲脂性物质和催化 PGD2 合成的双重功能，在脑脊液、精浆及血浆内都可检测到[25-26]。关于LPGDS如何影响雄性生殖能力的研究还没有清楚解析，但推测其携带肾上腺素、类固醇、维生素A等越过血-睾屏障帮助精子细生成与成熟。本研究中LPGDS基因在老年期小鼠睾丸组织的表达水平显著高于中年期，可以作为其功能研究的一个切入点。此外，补体系统、胰岛素、溶酶体、氨基酸代谢和嗅觉转导系统等检测到一些差异表达基因，其在睾丸生理生殖活动中的作用还需进一步研究。

本研究对中年期和老年期小鼠睾丸组织进行转录组分析，检测小鼠睾丸组织不同年龄阶段基因表达的差异及相关的信号通路。该研究为丰富小鼠睾丸组织转录组信息提供中老年期的数据，也为探讨睾丸功能退化及其在老龄化过程中生物学与医学等某些方面的研究奠定基础。