畜牧兽医学报  2023, Vol. 54 Issue (3): 989-1002. DOI: 10.11843/j.issn.0366-6964.2023.03.013    PDF    
引用本文
叶峻宁, 邓铭, 薛慧雯, 柳广斌, 邹娴, 孙宝丽, 郭勇庆, 刘德武, 李耀坤. 影响山羊胎儿肌肉发育mRNA和lncRNA的鉴定与分析[J]. 畜牧兽医学报, 2023, 54(3): 989-1002.  
YE Junning, DENG Ming, XUE Huiwen, LIU Guangbin, ZOU Xian, SUN Baoli, GUO Yongqing, LIU Dewu, LI Yaokun. Identification and Analysis of mRNA and lncRNA Affecting Goat Fetal Muscle Development[J]. Acta Veterinaria et Zootechnica Sinica, 2023, 54(3): 989-1002.  
影响山羊胎儿肌肉发育mRNA和lncRNA的鉴定与分析
叶峻宁1, 邓铭1, 薛慧雯1, 柳广斌1, 邹娴2, 孙宝丽1, 郭勇庆1, 刘德武1, 李耀坤1     
1. 华南农业大学 动物科学学院, 广州 510642;
2. 广东省农业科学院动物科学研究所, 畜禽育种国家重点实验室/广东省畜禽育种与营养研究重点实验室, 广州 510640
收稿日期：2022-08-04
基金项目：广东省现代农业产业技术体系(2019KJ127)；广东省科技特派员驻镇帮镇扶村项目(KTP20210344)
作者简介：叶峻宁(1998-)，男，广东东莞人，硕士，主要从事草食动物分子遗传的研究，E-mail: neilye@stu.scau.edu.cn.
通信作者：李耀坤, 主要从事遗传育种与繁殖及智慧养殖研究, E-mail: ykli@scau.edu.cn.
摘要：旨在分析雷州山羊胎儿70、90及120天背最长肌转录组测序中的差异表达mRNAs和lncRNAs及其靶基因，对差异表达基因进行功能富集分析，挖掘影响雷州山羊胎儿骨骼肌发育的关键候选基因。本研究选用年龄约2.5~3.5岁，饲养状况相同，健康无病，体重接近的雷州母山羊，分为妊娠期70、90和120天3组，每组各3只共9只，提取其胎儿背最长肌的RNA构建文库后使用Illumina HiSeqTM 2500进行转录组测序，以Q < 0.05为阈值筛选出差异mRNAs和lncRNAs并进行生物学功能富集分析及构建靶向关系网络。最后通过荧光定量PCR鉴定差异基因表达水平，验证测序结果的可靠性。结果发现，在M70 vs. M90组中，总共发现185个差异表达mRNAs，其中170个上调，15个下调，62个差异表达lncRNAs，40个上调，22个下调；在M120 vs. M90组中，共发现1 048个差异表达mRNAs，666个上调，382个下调，352个差异表达lncRNAs，175个上调，177个下调。其中TNNT3、TNNI1、TNNT1、TNNI2、ITGA4、ITGA11、ITGB4等基因参与肌肉发育相关的调控通路如肌钙蛋白复合物、肌丝、粘附斑、ECM-受体相互作用、PI3K-Akt信号通路等。根据lncRNA与mRNA的位置关系进行cis靶基因预测，总共筛选到37 464个靶基因，有10个差异表达mRNAs被差异表达lncRNA靶向，其中NDUFB2、BRAF、METTL7B和ITGA7参与粘着斑以及肌动蛋白细胞骨架的调控，METTL7B和ITGA7还参与PI3K-Akt信号通路，lncRNA TCONS_00214300靶向METTL7B和ITGA7。同时，8个差异表达mRNAs和lncRNAs的RT-PCR结果显示，8个基因的表达水平与转录组测序的表达趋势相似，说明测序结果可靠。结果显示，发现了影响雷州山羊胎儿肌肉发育的mRNA与lncRNA，发现了一些差异表达lncRNA和mRNA可能存在调控关系，对了解雷州山羊胎儿肌肉发育的分子调控机制提供了新的参考，也为雷州山羊分子选育提供了理论借鉴。
关键词：雷州山羊    mRNA    lncRNA    背最长肌    
Identification and Analysis of mRNA and lncRNA Affecting Goat Fetal Muscle Development
YE Junning1, DENG Ming1, XUE Huiwen1, LIU Guangbin1, ZOU Xian2, SUN Baoli1, GUO Yongqing1, LIU Dewu1, LI Yaokun1     
1. College of Animal Science, South China Agricultural University, Guangzhou 510642, China;
2. Guangdong Key Laboratory of Animal Breeding and Nutrition, State Key Laboratory of Livestock and Poultry Breeding, Institute of Animal Science, Guangdong Academy of Agricultural Sciences, Guangzhou 510640, China
Corresponding author: LI Yaokun, E-mail：ykli@scau.edu.cn.
Abstract: The study aimed to examine the differentially expressed mRNAs and lncRNAs as well as their target genes in the transcriptome sequencing of Leizhou goat fetuses at 70, 90 and 120 days, and to investigate important candidate genes controlling the development of skeletal muscle in Leizhou goat fetuses, functional enrichment analysis of differentially expressed genes was carried out. Leizhou female goats were separated into 3 groups of 70, 90, and 120 days of gestation, with 3 animals in each group, for a total of 9 animals. The goats were about 2.5-3.5 years old, in the same feeding condition, healthy, and disease-free. The transcriptome was sequenced using Illumina HiSeqTM 2500 after the RNA of fetal longissimus dorsi was extracted in order to establish a library. The differential mRNAs and lncRNAs were screened at a threshold of Q < 0.05, enriched for biological functions, and constructed into a target relationship network. In order to confirm the accuracy of the sequencing results, the expression levels of differential genes were lastly detected using fluorescence quantitative PCR. It was found that in the M70 vs. M90 group, a total of 185 differentially expressed mRNAs were identified, of which 170 were upregulated, 15 were downregulated, and 62 differentially expressed lncRNAs, 40 were upregulated, and 22 were downregulated. In the M120 vs. M90 group, we identified 1 048 differentially expressed mRNAs, 666 were upregulated, 382 were downregulated, 352 differentially expressed lncRNAs, 175 were upregulated, and 177 were downregulated. Among them, TNNT3, TNNI1, TNNT1, TNNI2, ITGA4, ITGA11, and ITGB4 genes were involved in regulatory pathways related to muscle development such as troponin complex, myofilament, adhesion plaque, ECM-receptor interaction, PI3K-Akt signaling pathways, etc. A total of 37 464 target genes were screened according to the positional relationship between lncRNA and mRNA for cis target gene prediction. Ten differentially expressed mRNAs were targeted by differentially expressed lncRNAs, among which NDUFB2, BRAF, METTL7B, and ITGA7 regulated adherent spots and the actin cytoskeleton, and METTL7B and ITGA7 were also engaged in PI3K-Akt signaling pathway. lncRNA TCONS_00214300 targeted METTL7B and ITGA7. Meanwhile, RT-PCR results of 8 differentially expressed mRNAs and lncRNAs showed that the expression levels of the 8 genes were similar to the expression trends of transcriptome sequencing, indicating that the sequencing results were reliable. The findings identify the lncRNAs and mRNAs that regulate fetal muscle development in Leizhou goats, as well as the potential regulatory linkages between some lncRNAs and mRNAs with differential expression. In addition, the studys also offer a theoretical foundation for molecular selection and breeding of Leizhou goats, and a novel source for studying the molecular regulatory mechanism of fetal muscle development in Leizhou goats.
Key words: Leizhou goat    mRNA    lncRNA    longissimus dorsi    

山羊是一种重要的农业动物，有着独特的风味，与其他家畜相比，山羊肉含有更多的蛋白质，更低的脂肪和胆固醇[1]。我国不仅是羊肉生产大国，更是羊肉消费大国。羊肉作为我国主要的畜产品之一，近年来产量持续稳定增长且长期居于世界第一。2020年初新冠肺炎疫情对肉羊产业造成了一定冲击，随着疫情防控转向常态化，羊肉产能恢复势头良好。2021年我国肉羊生产继续向好，羊肉产量514万吨，比2020年增加21.7吨，增幅4.4%，同时我国羊肉表观消费量达到554. 9万吨，比2020年增长5.0%[2]。由此而见，我国羊肉消费需求旺盛，虽然羊肉产能逐年增长，但是国内羊肉供给仍然紧张，供需缺口呈常态化。

肌生成是一个复杂且受严格控制的过程，涉及肌生成谱系的胚胎前体、成肌细胞生长和细胞周期退出，单核前体细胞在最后阶段融合生成多核的成熟肌管[3]。先前的研究发现，在肌肉发育中起关键作用的肌源性调节因子(MRFs)，如Myf5、MyoD、Myf6和肌生成素，控制着这一系列活动。Myf5和MyoD参与控制成肌细胞的产生、增殖和寿命，而Myf6和肌生成素则参与分化的最后阶段[4-5]。

长链非编码RNA (long noncoding RNA，lncRNA)，是动植物细胞中广泛存在的一类长度大于200 nt的RNA，具有无蛋白质编码功能、低表达和物种间保守性低的特点[6-7]。目前许多研究表明，lncRNA参与了许多生物进程如组蛋白修饰、基因表达调控、细胞增殖分化等[8-9]。越来越多的研究发现，lncRNA在骨骼肌发育的过程中发挥了重要的作用。lncRNA-CTTN-IT1在湖羊骨骼肌卫星细胞(SMSC)中作为ceRNA竞争性结合miR-29a来调节YAP1，可促进湖羊骨骼肌卫星细胞的增殖和分化以及作为研究湖羊肌肉细胞的新分子标记物[10]。在内蒙古黑牛中，lncRNA-403通过调节KRAS和MYF6，可能通过多通路网络调控方式参与牛骨骼肌卫星细胞的分化[11]。lncRNA-MEG3可以海绵化miRNA-135并影响MEF2C的表达，促进牛骨骼肌分化[12]。

目前，对于山羊肌肉的有关研究大部分仍处于初步阶段，即通过转录组测序寻找相关调控山羊肌肉生长发育的差异基因，Han等[13]通过对1月龄与9月龄的武安黑山羊进行RNA测序，发现了36个差异表达lncRNAs，构建了一个包含4个lncRNAs、3个miRNAs和8个mRNAs的调控网络；Zhan等[14]对多个日龄的简州大耳山羊胎儿进行RNA测序，发现了48个差异表达的反义lncRNAs，它们可能通过负调控靶基因的表达参与肌肉发育；Shen等[15]使用RNA测序比较了辽宁绒山羊和子午岭黑山羊背最长肌的miRNAs，在辽宁绒山羊中鉴定出159个差异表达miRNAs，可能是导致山羊产肉性能差异的重要基因。可能由于我国养羊业发展处于起步阶段，并未形成规模产业化，不少地区仍使用散养模式，山羊良种化程度低，由于育种需要长期投入，回报周期长，难以在短时间内取得经济利益，核心育种体系不完善已成为困扰和制约我国肉羊良种改良的障碍[16]。雷州山羊主要分布在中国广东省雷州半岛，是广东省一种本土山羊品种，具有耐粗饲和耐湿热，抗病能力强，性成熟早，肉质优越等特点[17]。然而，由于近亲繁殖、缺乏系统的本品种选育，加上饲养管理技术落后，品种退化严重，体格明显变小。所以，对雷州山羊的转录图谱进行研究，挖掘影响其肌肉发育的基因与通路势在必行。由于继发性肌生成是胎儿祖细胞融合产生次级纤维的过程，对胎儿骨骼肌发育至关重要，继发性肌生成对于胚胎发育过程中的肌肉生长和成熟(纤维型)至关重要，但出生后肌肉纤维的数量没有净增加[18-19]，所以本研究选择雷州山羊胎儿70、90、120天的背最长肌作为研究对象，采用RNA-Seq技术筛选出组织中影响肌肉发育及分化的相关mRNA和lncRNA及其靶基因，对差异表达基因进行功能富集分析，对今后阐明雷州山羊骨骼肌发育分子机制提供理论依据，为雷州山羊分子育种、开发提供新思路。

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

试验所使用的的雷州山羊来自湛江雷州市雷州山羊保种场，采用半开放高床羊舍和散栏饲养模式，饲喂条件为青贮玉米秸秆为主，辅以混合精料。混合精料组成：玉米55%、豆粕18%、小麦麸18%、米糠8%、食盐0.3%、小苏打0.6%和预混料0.1%。选择9只健康无病，体重为26~30 kg，年龄为2.5~3.5岁，妊娠期在70、90、120天的雷州母山羊，分为3组，每组3只。屠宰后采集胎儿第10肋骨处背最长肌组织样品，将其装入冻存管并立即放入液氮中，样品带回实验室后放入－80 ℃冰箱保存。

1.2 总RNA提取和建库

采用Trizol法按照说明书提取每个背最长肌组织样品的总RNA。用1%琼脂糖凝胶检测RNA的降解和污染程度，使用纳米分光光度计(IMPLEN，CA，USA)检查RNA纯度，使用Agilent Bioanalyzer 2100系统的RNA Nano 6000检测试剂盒评估RNA的总量和完整性。选择符合建库标准的RNA送至北京诺禾致源科技股份有限公司测序，采用Illumina HiSeqTM2500测序平台进行双端150 bp高通量测序，并委托北京诺禾致源科技股份有限公司对测序结果进行分析。每个样品分别取1 μg总RNA构建文库，根据NEB Next Ultra Directional RNA Library Prep Kit for Illumina(NEB，Ispawich，USA)的操作步骤，将不同的index进行分别标签建库。

1.3 数据质控与组装

Fastq格式的原始数据(raw reads)首先通过perl脚本进行处理。在此步骤中，通过修剪包含适配器、ploy-N或原始数据中低质量的读取，获得干净数据(clean reads)。同时，计算Q20、Q30和GC含量的干净数据。所有下游分析都是基于高质量的干净数据进行的。使用hisat2-2.0.4构建山羊参考基因组(GCF_001704415.1_ARS1_genomic；http://www.ensembl.org/)的索引，并将配对端干净读取映射到参考基因组。STAR使用最大可映射前缀(MMP)的方法，可以为交汇点读取生成精确的映射结果。

1.4 lncRNA筛选

使用Cuffmerge软件对各样本拼接得到的转录本进行合并，并去掉其中链方向不确定的转录本，得到本次测序完整的转录组信息。之后对合并的转录本集合与数据库进行比较，区分出已知转录本与新转录本。并进一步在新转录本中，根据lncRNA(长链非编码RNA)的结构特点(一类转录本长度超过200 nt)以及不编码蛋白的功能特点，设置一系列严格的筛选条件，通过下述5个步骤的筛选得到Novel_lncRNA：1)转录本exon个数筛选：过滤转录本拼接结果中低可信度的单外显子转录本, 选择exon个数大于等于2的转录本；2)转录本长度筛选：选择转录本长度大于200 bp的转录本；3)转录本已知注释筛选：通过Cuffcompare软件，筛除与数据库注释exon区域有重叠的转录本；4)编码潜能筛选：是否具有编码潜能是判断转录本是否为lncRNA的关键条件。针对拼接得到的转录本，综合目前主流的编码潜能分析方法进行该项筛选, 取这些软件分析结果中没有编码潜能的转录本的交集作为本次分析预测得到的候选Novel lncRNA数据集；5)参考HGNC(The HUGO Gene Nomenclature Committee)关于long non-coding RNA的命名指导标准对候选的Novel lncRNA根据其与编码基因的位置关系进行最终筛选及命名，得到本次分析的Novel lncRNA。

1.5 lncRNA和mRNA差异表达分析

使用DESeq2 R包(1.16.1)对3个时间组的背最长肌组织进行lncRNA和mRNA差异表达分析。以校正P值< 0.05为筛选显著差异表达基因标准。

1.6 lncRNA靶基因预测

位置相关性靶基因分析，根据lncRNA与mRNA的位置关系进行cis靶基因预测，筛选的范围为100K以内。

1.7 GO和KEGG功能富集分析

基因本体(Gene Ontology，GO)是一个提供Controlled Vocablary对生物体内基因和产物属性全面描述的国际标准化的基因功能分类体系。使用GOseq R包对差异表达lncRNAs和mRNAs进行GO富集分析，校正P < 0.05的GO项被认为是基因集GO功能的显著富集。京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)是一个集合了基因组、化学等系统功能信息对基因组进行处理的数据库，用于从分子级数据中推断出生物系统(如细胞、生物体和生态系统)的高级功能，特别是由基因组测序和其他高通量试验技术产生的大规模分子数据集(http://www.genome.jp/kegg/)。使用KOBAS来测试KEGG途径中差异表达基因的统计富集。

1.8 lncRNA与mRNA互作网络构建

使用Gephi(0.93)软件对M120 vs. M90组的差异表达lncRNAs靶向的差异表达mRNAs构建连接网络。

1.9 RT-PCR验证分析

从差异表达lncRNAs和mRNAs中各随机选择4个基因，利用RT-PCR进行转录组表达谱验证分析。PCR反应体系20 μL：10 μL SYBR® Green Master Mix，上、下游引物各0.5 μL，1 μL cDNA和8 μL ddH2O。反应条件：95 ℃ 10 min；95 ℃ 5 s，60 ℃ 20 s，72 ℃ 1 min，40个循环。引物设计参考NCBI的山羊相关基因序列，GAPDH作为内参基因，用Primer 5.0软件设计，并由华大基因合成引物(表 1)。

表 1 引物序列 Table 1 Primer sequences
2 结果 2.1 测序数据的质量及比对信息

处理测序数据，去除测序raw reads中低质量的reads，最终得到妊娠70日组(M70-1、M70-2、M70-3)、妊娠90日组(M90-1、M90-2、M90-3)和妊娠120日组(M120-1、M120-2、M120-3)的clean reads分别为10 269 092、10 421 322、10 079 615、10 168 621、10 315 943、10 626 230、10 699 782、10 367 621和12 235 445，9个样品中的GC含量为48.95%~51.66%，Q30的比例在95.92%以上(表 2)。

表 2 9个山羊胎儿背最长肌样本测序数据质量 Table 2 Quality of sequencing data of 9 goat fetal longissimus dorsi samples
2.2 差异表达mRNA和lncRNA鉴定

使用edgeR软件进行差异表达分析，使用层次聚类方法对样本的表达值进行聚类，在M70 vs. M90组中，总共发现185个差异表达mRNAs，其中170个上调，15个下调；62个差异表达lncRNAs，40个上调，22个下调(图 1A、图 1B)(Q < 0.05)。在M120 vs. M90组中，共发现1 048个差异表达mRNAs，666个上调，382个下调；352个差异表达lncRNAs，175个上调，177个下调(图 1C、图 1D)(Q < 0.05)。

A. 差异表达mRNAs火山图(M70 vs. M90)；B. 差异表达mRNAs火山图(M120 vs. M90)；C. 差异表达lncRNAs火山图(M70 vs. M90)；D.差异表达lncRNAs火山图(M120 vs. M90) A. Differentially expressed mRNAs volcano map (M70 vs. M90); B. Differentially expressed mRNAs volcano map (M120 vs. M90); C. Differentially expressed lncRNAs volcano map (M70 vs. M90); D. Differentially expressed lncRNAs volcano map (M120 vs. M90) 图 1 雷州山羊胎儿不同阶段差异表达的mRNAs和lncRNAs Fig. 1 The number of mRNAs and lncRNAs expressed in different stages of the Leizhou goat fetus
2.3 差异表达mRNA功能分析

为了解差异表达基因涉及的生物学功能，对差异表达基因进行GO条目和KEGG通路的功能富集分析。GO富集结果显示，在M70 vs. M90组中，差异基因主要富集在生物学进程、分子功能和细胞组成的3大类的1 358个GO条目中，其中有75个显著富集(P < 0.05)。主要参与细胞外区域、细胞氨基酸分解代谢过程、内肽酶活性调节、有机酸代谢过程、催化活性的负调节等(图 2A)。在M120 vs. M90组中，差异表达基因显著富集在22个GO条目中，包括钙离子结合、肌钙蛋白复合物、肌丝、细胞黏附、生物黏附、肌原纤维、肌动蛋白丝结合等(P < 0.05)(图 2B)。其中，TNNT3、TNNI1、TNNT1、TNNI2都富集于7个与肌肉有关的细胞成分中。KEGG富集结果显示，在M70 vs. M90组，差异基因显著富集在12个通路中(P < 0.05)，包括补体与凝血级联、苯丙氨酸代谢、氨基酸生物合成、代谢途径等(图 3A)。在M120 vs. M90组中，差异基因富集在ECM-受体相互作用、黏着斑、细胞黏附分子、心肌收缩、p53信号通路、蛋白质消化和吸收等25个信号通路中(P < 0.05)(图 3B)。ITGA4、ITGA11、ITGB4同时富集于黏附斑、ECM-受体相互作用、PI3K-Akt信号通路及多个与心肌相关的通路中。

A. 差异表达mRNAs富集分析中前22个GO条目(M70 vs. M90)；B. 差异表达mRNAs富集分析中前22个GO条目(M120 vs. M90)。纵坐标代表GO条目；横坐标代表基因数量。柱子长度代表在条目中的基因数。不同颜色代表P值，颜色越红P值越小。BP. 生物进程；CC. 细胞成分；MF. 分子功能。下同 A. The top 22 GO terms in the enrichment analysis of differentially expressed mRNAs (M70 vs. M90); B. The top 22 GO terms in the enrichment analysis of differentially expressed mRNAs (M120 vs. M90).The vertical coordinates represent GO terms; the horizontal coordinates represent the number of genes. The bars length represent the number of genes in the term. Different colours represent P-value; the redder the colour the smaller the P-value. BP. Biological process; CC. Cellular component; MF. Molecular function. The same as below 图 2 差异表达mRNAs GO富集分析 Fig. 2 GO enrichment analysis of differentially expressed mRNAs
A. 差异表达mRNAs富集分析中前20条KEGG通路(M70 vs. M90)；B. 差异表达mRNAs富集分析中前20条KEGG通路(M120 vs. M90)。纵坐标表示KEGG通路；横坐标表示富集因子。圆的大小代表在通路的基因数，圆越大基因数越多。不同颜色代表Q值，颜色越红，Q值越小。下同 A. The top 20 KEGG pathways in the enrichment analysis of differentially expressed mRNAs (M70 vs. M90); B. The top 20 KEGG pathways in the enrichment analysis of differentially expressed mRNAs (M120 vs. M90). The vertical coordinate indicates the KEGG pathway; the horizontal coordinate indicates the rich factor. The size of the circle represents the number of genes in the pathway, the larger the circle the more genes. The different colours represent Q value, the redder the colour, the smaller the Q value. The same as below 图 3 差异表达mRNAs KEGG富集分析 Fig. 3 KEGG enrichment analysis of differentially expressed mRNAs
2.4 差异表达lncRNA靶基因预测及功能分析

位置相关性靶基因分析，根据lncRNA与mRNA的位置关系进行cis靶基因预测，筛选的范围为100K以内。总共筛选到37 464个靶基因，在M70 vs. M90组中，差异表达的lncRNAs靶向了60个mRNAs，无差异表达mRNA；在M120 vs. M90组中，有161个mRNAs被差异表达的lncRNAs靶向，有10个mRNAs为差异表达mRNAs。对差异表达lncRNAs的靶基因进行GO富集分析，在M70 vs. M90组中，靶基因主要富集于核酸结合转录因子活性、蛋白质转运、蛋白质定位、mRNA和miRNA分解代谢、肌钙蛋白复合物、肌丝等(P < 0.05)(图 4A)；在M120 vs. M90组中，靶基因主要富集于G蛋白偶联受体信号通路、感觉感知、对化学刺激的反应、信号转导、细胞表面受体信号通路等(P < 0.05)(图 4B)。KEGG富集分析结果表明，在M70 vs. M90组，差异靶基因共富集到11个信号通路中，主要富集于嗅觉转导(P < 0.05)(图 5A)；在M120 vs. M90组中，主要富集于嗅觉转导、cAMP信号通路、钙信号通路、PI3K-Akt信号通路、肌动蛋白细胞骨架的调控等，其中NDUFB2、BRAF、METTL7B和ITGA7参与黏着斑以及肌动蛋白细胞骨架的调控，METTL7B和ITGA7还参与PI3K-Akt信号通路(P < 0.05)(图 5B)(表 3)。

A.差异表达lncRNAs富集分析中前19个GO条目(M70 vs. M90)；B.差异表达lncRNAs富集分析中前22个GO条目(M120 vs. M90) A. The top 19 GO terms in the enrichment analysis of differentially expressed lncRNAs (M70 vs. M90); B. The top 22 GO terms in the enrichment analysis of differentially expressed lncRNAs (M120 vs. M90) 图 4 差异表达lncRNAs GO富集分析 Fig. 4 GO enrichment analysis of differentially expressed lncRNAs
A. 差异表达lncRNAs富集分析中前11条KEGG通路(M70 vs. M90)；B. 差异表达lncRNAs富集分析中前20条KEGG通路(M120 vs. M90) A. The top 11 KEGG pathways in the enrichment analysis of differentially expressed lncRNAs (M70 vs. M90); B. The top 20 KEGG pathways in the enrichment analysis of differentially expressed lncRNAs (M120 vs. M90) 图 5 差异表达lncRNAs KEGG富集分析 Fig. 5 KEGG enrichment analysis of differentially expressed lncRNAs
表 3 M120 vs. M90组与肌肉发育相关差异表达基因汇总 Table 3 Summary of differentially expressed genes related to muscle development in M120 vs. M90 group
2.5 lncRNA与mRNA互作网络分析

使用Gephi软件将M120 vs. M90组中差异表达的lncRNAs与其靶向的差异表达mRNAs构建互作网络。METTL11B被18个lncRNAs靶向，TCONS_00214300靶向METTL7B和ITGA7，METTL7B位于TCONS_00214300上游，ITGA7位于下游(图 6)。

被靶向的基因越多，圆越大 The more genes that are targeted, the larger the circle 图 6 M120 vs. M90组差异表达lncRNAs靶向差异表达mRNAs互作图 Fig. 6 Interaction diagram of differentially expressed lncRNAs targeted differentially expressed mRNAs in M120 vs. M90 group
2.6 测序结果RT-PCR验证

通过荧光定量对随机挑选的8个差异表达mRNAs和lncRNAs进行验证。结果表明，随机挑选的8个mRNAs和lncRNAs的荧光定量结果与转录组测序结果的表达趋势相似(图 7)，证实雷州山羊胎儿背最长肌组织转录组测序结果具有可靠性，可用于后续研究。

图 7 差异表达基因RT-PCR验证 Fig. 7 RT-PCR verification of differentially expressed genes
3 讨论

胎儿期是骨骼肌发育的关键时期，肌纤维的数量主要在胎儿时期增加，出生后主要表现为骨骼肌的增大[20]。因此，研究胎儿时期骨骼肌的发育，对提高家畜的产肉量是非常有意义的。本研究选用雷州山羊妊娠70、90、120天的胎儿作为研究对象，对其背最长肌进行RNA-Seq测序，筛选与肌肉发育及分化有关的基因。在M70 vs. M90与M120 vs. M90两组间发现了一些显著差异的信号通路和生物学功能与肌肉有关，包括黏附斑(focal adhesion)、磷脂酰肌醇3-激酶/丝氨酸-苏氨酸蛋白激酶(phosphoinositide-3-kinase/serine-threonine kinase, PI3K-Akt) 信号通路、ECM-受体相互作用(extracellular matrix-receptor interaction)等，同时结合生物学功能分析发现了一些与雷州山羊胎儿背最长肌肌肉发育及分化可能相关的mRNAs和lncRNAs。

在差异表达的mRNAs中发现ITGA4、ITGA11、ITGB4富集于黏附斑、ECM-受体相互作用、PI3K-Akt信号通路及多个与心肌相关的通路中。细胞外基质(extracellular matrix，ECM)广泛存在于动物全身组织中，是一种由细胞分泌组装形成的复杂网络，不仅提供嵌入细胞的物理支架，调节许多细胞进程包括生长、分化、迁移、生存等[21-22]，还是生长因子的重要来源，通过提供许多生长、信号因子的附着和释放，达到帮助细胞间信号转导的功能[23]。细胞外基质是黏着斑和PI3K-Akt通路的重要组成部分。黏着斑调节骨骼肌发育，并参与几个关键信号通路如WNT、MAPK和PI3K-Akt[24]，同时与ECM之间的联系共同调节细胞运动、增殖和分化等多种细胞内途径[25]。在研究中发现PI3K-Akt通路可以促进骨骼肌增殖与分化[26-27]，是公认的骨骼肌生长和肥大的正调节剂[28]。

ITGA4、ITGA11、ITGB4属于整联蛋白家族，充当着细胞外基质和细胞骨架之间的桥梁，检测细胞微环境的变化，使细胞能根据外部环境做出反应[29]，在细胞增殖、迁移、分化和肿瘤转移侵袭中发挥着重要作用[30]。ITGA11在小鼠C2C12细胞系中介导细胞黏附到胶原蛋白I和IV，并在胶原蛋白上形成焦点接触，参与细胞迁移和胶原重组[31]。在胃癌组织中，ITGA11通过PI3K-Akt通路影响胃癌细胞，敲低ITGA11可抑制胃癌细胞增殖、侵袭和迁移[32]。ITGB4在肝细胞癌组织中表达水平增高，沉默ITGB4可以抑制细胞增殖、集落形成能力和细胞侵袭性[33]。目前，有关ITGA4、ITGA11、ITGB4的研究主要与癌症相关，鲜有与骨骼肌相关的研究。通过了解相关富集信号通路的功能以及ITGA11在小鼠C2C12细胞中的研究，猜测ITGA4、ITGA11、ITGB4可能参与到山羊骨骼肌细胞的增殖和迁移，其具体功能还有待深入研究。

差异表达lncRNA的靶基因富集分析中，发现NDUFB2、BRAF、METTL7B和ITGA7参与黏着斑以及肌动蛋白细胞骨架的调控通路，METTL7B和ITGA7也参与PI3K-Akt通路。同时，METTL7B和ITGA7是lncRNA TCONS_00214300的靶基因。

BRAF编码一种属于丝氨酸/苏氨酸蛋白激酶RAF家族的蛋白质。这种蛋白质在调节MAP激酶/ERK信号通路中起作用，影响细胞分裂、分化和分泌。在有关小鼠肌肉的研究中，BRAF在小鼠肢体骨骼肌的形成中起着关键作用，作为成肌细胞迁移的有效诱导剂，通过直接激活PAX3促进了C2C12成肌细胞迁移[34]。BRAFL579VCFC小鼠模型中，小鼠表现出肌纤维过小，BRAF L579V错义蛋白可抑制胚胎肌生成与成肌细胞分化，同时改变p38 MAPK通路的活性从而损害骨骼肌的形成和功能[35]。

ITGA7参与各种细胞与细胞，细胞与细胞外基质(ECM)间的相互作用，包括细胞的生长、存活与凋亡[36]。整合素alpha-7是由ITGA7编码形成的蛋白质，是一种异二聚体整合的膜蛋白，由α和β链组成，介导细胞迁移、形态发育、分化和转移[37]。ITGA7是骨骼肌中的层粘连蛋白受体，将细胞外基质连接到内部肌动蛋白细胞骨架[38]。ITGA7在之前的研究中已被证明为肌肉干细胞(muscle stem cells, MuSC)表面标志物，影响肌源性祖细胞的生成和肌肉分化[39-41]。

METTL7B是甲基转移酶样蛋白家族的成员，被认为与高尔基体蛋白的甲基化和脂质代谢有关[42]，参与各种恶性肿瘤的发展、侵袭与迁移[43]，目前暂未见与肌肉生长发育有关的研究。NDUFB2编码的蛋白属于线粒体呼吸链复合物I(NADH-泛醌氧化还原酶)的45个亚基中的一个亚基，在线粒体内表达量最高[44]，其编码的蛋白具有NADH脱氢酶活性和氧化还原酶活性，在电子从NADH转移到呼吸链过程中起着重要的作用[45]。目前，国内外关于NDUFB2的研究十分有限，对肌肉细胞的影响尚不清楚。因此，METTL7B与NDUFB2对肌肉生长发育的影响还有待后续深入研究。

4 结论

本研究对雷州山羊胎儿70、90及120天的背最长肌进行高通量测序，对差异表达基因进行GO和KEGG功能分析与信号通路富集发现，ITGA4、ITGA11、ITGB4、TNNT3、TNNI1、TNNT1、TNNI2等可通过参与肌钙蛋白复合物、肌丝、黏附斑、ECM-受体相互作用、PI3K-Akt等信号通路和生物代谢途径影响山羊肌肉发育。结合lncRNA靶基因预测发现，NDUFB2、BRAF、METTL7B和ITGA7参与粘着斑以及肌动蛋白细胞骨架的调控通路，并且METTL7B和ITGA7可与lncRNA TCONS_00214300互作，从而对山羊肌肉发育产生作用。本研究结果为了解雷州山羊胎儿肌肉发育的分子调控机制提供了新的参考，也为雷州山羊分子选育提供了理论借鉴。

参考文献
[1]
TEIXEIRA A, SILVA S, RODRIGUES S. Advances in sheep and goat meat products research[J]. Adv Food Nutr Res, 2019, 87: 305-370.
[2]
甘春艳, 李军, 金海. 2021年我国肉羊产业发展概况、未来发展趋势及建议[J]. 中国畜牧杂志, 2022, 58(3): 258-263.
GAN C Y, LI J, JIN H. Overview of the development of the meat sheep industry in China in 2021, future development trends and suggestions[J]. Chinese Journal of Animal Science, 2022, 58(3): 258-263. (in Chinese)
[3]
SHI J A, LI W D, LIU A F, et al. MiRNA sequencing of embryonic myogenesis in chengkou mountain chicken[J]. BMC Genomics, 2022, 23(1): 571. DOI:10.1186/s12864-022-08795-z
[4]
ANTONIOU A, MASTROYIANNOPOULOS N P, UNEY J B, et al. miR-186 inhibits muscle cell differentiation through myogenin regulation[J]. J Biol Chem, 2014, 289(7): 3923-3935. DOI:10.1074/jbc.M113.507343
[5]
POWNALL M E, GUSTAFSSON M K, EMERSON C P JR. Myogenic regulatory factors and the specification of muscle progenitors in vertebrate embryos[J]. Annu Rev Cell Dev Biol, 2002, 18: 747-783. DOI:10.1146/annurev.cellbio.18.012502.105758
[6]
李耀坤, 许香萍, 陆婷婷, 等. 影响山羊产羔数lncRNA的筛选及功能分析[J]. 畜牧兽医学报, 2020, 51(8): 1853-1865.
LI Y K, XU X P, LU T T, et al. Screening and functional analysis of lncRNA affecting litter size in goats[J]. Acta Veterinaria et Zootechnica Sinica, 2020, 51(8): 1853-1865. (in Chinese)
[7]
孟珊, 杨阳, 李睿霄, 等. lncRNA-6617调控猪肌内前体脂肪细胞分化的筛选与功能研究[J]. 畜牧兽医学报, 2022, 53(6): 1712-1722.
MENG S, YANG Y, LI R X, et al. Screening and functional study of lncRNA-6617 regulating porcine intramuscular preadipocytes differentiation[J]. Acta Veterinaria et Zootechnica Sinica, 2022, 53(6): 1712-1722. (in Chinese)
[8]
WAS N, SAUER M, FISCHER U, et al. lncRNA Malat1 and miR-26 cooperate in the regulation of neuronal progenitor cell proliferation and differentiation[J]. RNA, 2022. DOI:10.1261/rna.079436.122
[9]
WU W T, CAO L, JIA Y F, et al. Emerging roles of miRNA, lncRNA, circRNA, and their cross-talk in pituitary adenoma[J]. Cells, 2022, 11(18): 2920. DOI:10.3390/cells11182920
[10]
WU T Y, WANG S H, WANG L H, et al. Long Noncoding RNA (lncRNA) CTTN-IT1 elevates skeletal muscle satellite cell proliferation and differentiation by acting as ceRNA for YAP1 through absorbing miR-29a in Hu sheep[J]. Front Genet, 2020, 11: 843. DOI:10.3389/fgene.2020.00843
[11]
ZHANG X J, CHEN M M, LIU X F, et al. A novel lncRNA, lnc403, involved in bovine skeletal muscle myogenesis by mediating KRAS/Myf6[J]. Gene, 2020, 751: 144706. DOI:10.1016/j.gene.2020.144706
[12]
LIU M, LI B, PENG W W, et al. LncRNA-MEG3 promotes bovine myoblast differentiation by sponging miR‐135[J]. J Cell Physiol, 2019, 234(10): 18361-18370. DOI:10.1002/jcp.28469
[13]
HAN H Y, WANG X W, LI W T, et al. Identification and characterization of lncRNAs expression profile related to goat skeletal muscle at different development stages[J]. Animals, 2022, 12(19): 2683. DOI:10.3390/ani12192683
[14]
ZHAN S Y, XUE Y N, YANG L, et al. Transcriptome analysis reveals long non-coding natural antisense transcripts involved in muscle development in fetal goat (Capra hircus)[J]. Genomics, 2022, 114(2): 110284. DOI:10.1016/j.ygeno.2022.110284
[15]
SHEN J Y, HAO Z Y, LUO Y Z, et al. Deep small RNA sequencing reveals important miRNAs related to muscle development and intramuscular fat deposition in Longissimus dorsi muscle from different goat breeds[J]. Front Vet Sci, 2022, 9: 911166. DOI:10.3389/fvets.2022.911166
[16]
刘恩民, 卢增奎, 乐祥鹏. 我国养羊业现状及未来发展思考[J]. 中国畜牧业, 2018(9): 34-35.
LIU E M, LU Z K, LE X P. Reflections on the current situation and future development of sheep farming in China[J]. China Animal Industry, 2018(9): 34-35. (in Chinese)
[17]
邱金戈, 宋佳杰, 刘德武, 等. 华南地区舍饲川中黑山羊与雷州山羊的生产繁殖性能及死亡状况研究[J]. 中国畜牧杂志, 2021, 57(2): 206-208, 214.
QIU J G, SONG J J, LIU D W, et al. Study on the reproductive performance and mortality of housed Chuanzhong black goats and Leizhou goats in South China[J]. Chinese Journal of Animal Science, 2021, 57(2): 206-208, 214. (in Chinese)
[18]
ZHU M J, FORD S P, MEANS W J, et al. Maternal nutrient restriction affects properties of skeletal muscle in offspring[J]. J Physiol, 2006, 575(Pt 1): 241-250.
[19]
WANG H, NOULET F, EDOM-VOVARD F, et al. Bmp signaling at the tips of skeletal muscles regulates the number of fetal muscle progenitors and satellite cells during development[J]. Dev Cell, 2010, 18(4): 643-654. DOI:10.1016/j.devcel.2010.02.008
[20]
KARUNARATNE J F, ASHTON C J, STICKLAND N C. Fetal programming of fat and collagen in porcine skeletal muscles[J]. J Anat, 2005, 207(6): 763-768. DOI:10.1111/j.1469-7580.2005.00494.x
[21]
MUNCIE J M, WEAVER V M. The physical and biochemical properties of the extracellular matrix regulate cell fate[J]. Curr Top Dev Biol, 2018, 130: 1-37.
[22]
FRANTZ C, STEWART K M, WEAVER V M. The extracellular matrix at a glance[J]. J Cell Sci, 2010, 123(Pt 24): 4195-4200.
[23]
BIRCH H L. Extracellular matrix and ageing[J]. Subcell Biochem, 2018, 90: 169-190.
[24]
ROMER L H, BIRUKOV K G, GARCIA J G N. Focal adhesions: paradigm for a signaling nexus[J]. Circ Res, 2006, 98(5): 606-616.
[25]
LASSITER D G, NYLÉN C, SJÖGREN R J O, et al. FAK tyrosine phosphorylation is regulated by AMPK and controls metabolism in human skeletal muscle[J]. Diabetologia, 2018, 61(2): 424-432.
[26]
SHEN X X, WEI Y H, YOU G S, et al. Circular PPP1R13B RNA promotes chicken skeletal muscle satellite cell proliferation and differentiation via targeting miR-9-5p[J]. Animals (Basel), 2021, 11(8): 2396.
[27]
YOSHIDA T, DELAFONTAINE P. Mechanisms of IGF-1-mediated regulation of skeletal muscle hypertrophy and atrophy[J]. Cells, 2020, 9(9): 1970.
[28]
GLASS D J. PI3 kinase regulation of skeletal muscle hypertrophy and atrophy[M]//ROMMEL C, VANHAESEBROECK B, VOGT P K. Phosphoinositide 3-Kinase in Health and Disease. Berlin: Springer, 2010: 267-278.
[29]
WU P C, WANG Y Y, WU Y J, et al. Expression and prognostic analyses of ITGA11, ITGB4 and ITGB8 in human non-small cell lung cancer[J]. PeerJ, 2019, 7: e8299.
[30]
FANG T Y, YIN X, WANG Y F, et al. Lymph node metastasis-related gene ITGA4 promotes the proliferation, migration, and invasion of gastric cancer cells by regulating tumor immune microenvironment[J]. J Oncol, 2022, 2022: 1315677.
[31]
TIGER C F, FOUGEROUSSE F, GRUNDSTRÖM G, et al. α11β1 integrin is a receptor for interstitial collagens involved in cell migration and collagen reorganization on mesenchymal nonmuscle cells[J]. Dev Biol, 2001, 237(1): 116-129.
[32]
ZHANG H J, ZHANG L, LU M. Inhibition of integrin subunit alpha 11 restrains gastric cancer progression through phosphatidylinositol 3-kinase/Akt pathway[J]. Bioengineered, 2021, 12(2): 11909-11921.
[33]
LI X L, LIU L, LI D D, et al. Integrin β4 promotes cell invasion and epithelial-mesenchymal transition through the modulation of Slug expression in hepatocellular carcinoma[J]. Sci Rep, 2017, 7: 40464.
[34]
SHIN J, WATANABE S, HOELPER S, et al. BRAF activates PAX3 to control muscle precursor cell migration during forelimb muscle development[J]. eLife, 2016, 5: e18351.
[35]
MAEDA Y, TIDYMAN W E, ANDER B P, et al. Ras/MAPK dysregulation in development causes a skeletal myopathy in an activating BrafL597V mouse model for cardio-facio-cutaneous syndrome[J]. Dev Dyn, 2021, 250(8): 1074-1095.
[36]
WU Z Y, KONG X Y, WANG Z H. Integrin α7 knockdown suppresses cell proliferation, migration, invasion and EMT in hepatocellular carcinoma[J]. Exp Ther Med, 2021, 21(4): 309.
[37]
DESGROSELLIER J S, CHERESH D A. Integrins in cancer: biological implications and therapeutic opportunities[J]. Nat Rev Cancer, 2010, 10(1): 9-22.
[38]
HELLER K N, MONTGOMERY C L, SHONTZ K M, et al. Human α7 integrin gene (ITGA7) delivered by adeno-associated virus extends survival of severely affected dystrophin/utrophin-deficient mice[J]. Hum Gene Ther, 2015, 26(10): 647-656.
[39]
BURKIN D J, WALLACE G Q, MILNER D J, et al. Transgenic expression of α7β1 integrin maintains muscle integrity, increases regenerative capacity, promotes hypertrophy, and reduces cardiomyopathy in dystrophic mice[J]. Am J Pathol, 2005, 166(1): 253-263.
[40]
DOMINICI C, RICHARD S. Muscle stem cell polarity requires QKI-mediated alternative splicing of Integrin Alpha-7 (Itga7)[J]. Life Sci Alliance, 2022, 5(5): e202101192.
[41]
ROONEY J E, GURPUR P B, YABLONKA-REUVENI Z, et al. Laminin-111 restores regenerative capacity in a mouse model for α7 integrin congenital myopathy[J]. Am J Pathol, 2009, 174(1): 256-264.
[42]
JIANG Z P, YIN W, ZHU H C, et al. METTL7B is a novel prognostic biomarker of lower-grade glioma based on pan-cancer analysis[J]. Cancer Cell Int, 2021, 21(1): 383.
[43]
SONG H B, LIU D C, WANG L W, et al. Methyltransferase like 7B is a potential therapeutic target for reversing EGFR-TKIs resistance in lung adenocarcinoma[J]. Mol Cancer, 2022, 21(1): 43.
[44]
HENDRICKSON S L, LAUTENBERGER J A, CHINN L W, et al. Genetic variants in nuclear-encoded mitochondrial genes influence AIDS progression[J]. PLoS One, 2010, 5(9): e12862.
[45]
LAUGHLIN T G, BAYNE A N, TREMPE J F, et al. Structure of the complex I-like molecule NDH of oxygenic photosynthesis[J]. Nature, 2019, 566(7744): 411-414.

(编辑   郭云雁)