水貂(Neovison vison)是一种小型珍贵毛皮动物，自引入至今已有60多年的历史，目前我国已成为世界最大的水貂饲养国。但是与国外相比，我国种貂一直存在体型较小、毛绒质量较差的问题，优质种源一直受制于欧洲、北美等地区国家，长期以来陷入了“引种-退化-再引种”的恶性循环，对水貂养殖业造成极大的经济负担。要从根本上改变现状，就必需培育我国自主知识产权的体型大、毛绒品质好的优良品种，所以研究水貂生长发育机制是水貂养殖业发展大势所趋。
近年来，随着大数据的迅猛发展，转录组测序技术已经广泛应用于各项研究中。例如在羊[1]、牛[2]、猪[3]等动物上已经利用转录组技术研究了生长发育，并且得到了关键基因及关联代谢通路。例如Chao等[4]利用转录组发现了29种lncRNAs与绵羊肌肉发育、代谢、细胞增殖和凋亡有关。Li等[5]构建了牛胚胎期和成年期的肌肉组织基因表达图谱，筛选到阶段特异表达基因，并发现GosB基因可以作为改变成肌细胞数量的育种候选基因。Ropka-Molik等[6]利用转录组和miRNAome数据检测筛选到与猪肌肉发育相关的43个基因和18个miRNAs。而对于水貂这种研究基础非常薄弱的特种经济动物来说，组学应用较少。据目前报道，水貂基因组长度为2.4 Gb[7]，利用转录组筛选到毛色相关基因9个[8]，揭示了水貂胚胎着床机制[9], 还未见到转录组学在水貂生长发育方面的研究。
本研究利用转录组测序技术比较分析美国短毛黑水貂快速生长过程(45和90日龄)中胸肌肌肉组织差异表达基因，筛选与生长发育相关的候选基因，为研究水貂生长发育调控机制及培育大体型水貂品种奠定基础。
1 材料与方法 1.1 样本采集试验所用6只健康美国短毛黑水貂公貂为大连名威貂业有限公司提供，分别于45(貂号为219、249、271)、90日龄(貂号为219(2)、249(2)、271(2))猝死后采集胸肌肌肉样品，每个阶段3个重复个体，采集后迅速放入液氮中。
1.2 水貂肌肉组织总RNA提取将6份胸肌样品分别用AMBION动物试剂盒提取总RNA。Nanodrop 2000检测RNA纯度，Agilent 2100检测RNA完整性。检测合格后放入－80 ℃冰箱保存备用。
1.3 水貂肌肉组织转录组测序文库的构建及测序提取样品总RNA，使用DNase消化DNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA；加入打断试剂将mRNA打断成短片段，以打断后的mRNA为模板, 用六碱基随机引物合成一链cDNA，然后合成二链cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100 Bioanalyzer质检合格后。测序工作由上海欧易生物公司使用Illumina HiSeqTM2500测序平台完成。
1.4 转录组测序数据质量评估为了保证数据质量，在信息分析前使用FASTQC软件对原始数据进行质量评估。并且通过NGS QC TOOLKIT v2.3.3[10]软件对数据预处理以去除低质量的片段。
1.5 DE NOVO拼接以及UNIGENE注释虽然水貂的基因组已经公布，但是由于注释信息还不完善，所以本研究采用De novo拼接方法。使用Trinity[11]软件paired-end的拼接方法得到Transcript序列，之后利用TGICL[12]软件聚类去冗余延伸得到一套最终的Unigene，以此作为后续分析的参考序列。然后基于BLAST算法[13]将Unigene序列分别与NR、SWISSPROT和KOG库进行比对，取e < 1×10-5的注释。
1.6 差异显著基因筛选及功能富集分析本研究以45日龄为对照组，筛选显著差异的条件为校正后P < 0.05，|log2FC|>1。差异基因筛选后，对其进行GO和KEGG分析。
1.7 测序结果qRT-PCR验证 1.7.1 引物设计与合成为了验证转录组测序结果准确性，本研究随机选取了6个差异表达基因(ZNF106、TFRC、GDF11、CCNB2、CCL24、FLNB)，参照转录组测序结果序列设计引物，进行qRT-PCR验证。引物序列如表 1所示。
|
表 1 qRT-PCR验证差异表达基因引物信息 Table 1 qRT-PCR primers of 6 selected differentially expressed genes |
采用SBYR Green染料法对测序结果准确性进行验证，以β-actin作为内参基因。反应体系：SYBR Green 10 μL，上下游引物各0.8 μL，cDNA 2 μL，水6.4 μL，共20 μL体积。反应条件：95 ℃ 60 s；95 ℃ 15 s，60 ℃ 15 s，72 ℃ 45 s，35个循环进行。
2 结果 2.1 水貂肌肉组织总RNA质量检测45和90日龄水貂的肌肉组织总RNA质量检测电泳图见图 1，NanoDrop分光光度计和Agilent 2100检测结果见表 2。由结果可以看出样品均达到建库要求。
|
图 1 总RNA电泳图 Fig. 1 Electrophoresis result of the total RNA |
|
|
表 2 总RNA质量检测 Table 2 The detection results of RNA quality |
通过Solexa RNA的paired-end测序，得到了大量的样本数据。对其进行预处理后数据结果统计见表 3，从表中数据可以看出，测序数据质量较高，可用于后续分析使用。
|
|
表 3 质量预处理后数据结果统计 Table 3 Data result statistics after quality pretreatment |
通过对美国短毛黑水貂45与90日龄的肌肉转录组数据分析，以45日龄为对照组，发现45与90日龄共有显著差异表达基因279个(以P < 0.05，|log2FC|>1)，火山图见图 2。
|
图 2 45~90日龄水貂肌肉组织差异表达基因火山图 Fig. 2 The volcano plot of differentially expressed genes in muscle of minks on 45-90 days old |
为验证转录组数据的准确性，采用qRT-PCR方法，从测序结果中随机选取6个差异表达基因，对转录组数据进行验证。选取上调基因ZNF106、TFRC、GDF11，下调基因CCNB2、CCL24、FLNB。结果表明，这6个基因qRT-PCR验证结果与转录组测序结果一致(图 3)，证明了本转录组测序结果的可靠性。
|
图 3 差异基因qRT-PCR验证 Fig. 3 qRT-PCR validation of differentially expressed genes |
对符合筛选条件的279个差异表达基因进行GO富集分析，对其功能进行描述。统计每个GO条目中所包括的差异基因个数，并用超几何分布检验方法计算每个GO条目中差异Unigene富集的显著性。GO条目满足P < 0.05显著富集的共有792条，其中554条为生物过程，包括干细胞的增殖、各器官发育的调控(输尿管、胰腺、淋巴结、肺、肾、小脑、脊髓、心脏、生殖系统)、肌肉的发育、骨骼肌组织发育、成骨细胞分化、增殖、骨矿化调控、骨化的负调控、膜内骨化、成骨细胞的增殖、软骨的发育、破骨细胞的分化、骨小梁的形成等；161条为分子功能，包括氧化还原酶活性、细胞周期调节活性、细胞因子的活性、肌球蛋白结合等；77条为细胞组分，包括线粒体膜、纺锤体微管、核膜等。3种GO条目的前20个见图 4。
|
图 4 GO功能富集 Fig. 4 GO enrichment analysis |
在553条生物学过程GO条目中筛选出与水貂肌肉生长相关的有66条(表 4)，根据GO条目间的层次关系构建GO-trees(图 5)，主要归集为发育过程、代谢过程、生物过程的正调控、生物过程的负调控、生物调节、生物过程调节、细胞的过程、多细胞生物过程、细胞成分组织或生物合成、生物粘附、定位、建立定位、信号、对刺激的反应、单一的生物过程等15类。从这66条GO条目中，共筛选得到差异表达基因44个(表 5)，其中上调20个，下调24个，这些差异基因可能在水貂的生长发育过程中起到调控作用。其中转录因子有2个，为JUNB、ATF3；转化生长因子有2个，为GDF11、TGFβ1；PPAPDC3、ABLIM1、CTSZ、TMSB10、TNNT2、MICAL1、FLNB等基因都是参与肌肉生长发育的相关基因，这些筛选到的差异基因可以作为后续研究水貂生长发育的候选基因。
|
|
表 4 筛选出的GO条目 Table 4 The selected significant enrichment GO terms |
|
蓝色代表level 3的GO条目。红色表示对应的level 2的GO条目 The blue represent the GO terms of level 3, the red represent the GO terms of level 2 图 5 GO条目层次关系树 Fig. 5 GO-tree |
|
|
表 5 筛选的差异表达基因 Table 5 The selected differentially expressed genes |
基于KEGG数据库，对差异表达基因进行了通路注释，结果发现，差异表达基因显著富集于47条信号通路。其中与肌肉生长发育相关的有12条，主要包括p53信号通路、PPAR信号通路、FOXO信号转导通路、酸代谢、氮代谢等。具体如表 6所示。
|
|
表 6 筛选与肌肉生长发育相关的通路 Table 6 The selected KEGG pathways related to muscle growth and development |
转录组研究是基因功能及结构研究的基础和出发点，通过新一代高通量测序，能够全面快速地获得某一物种特定组织或器官在某一状态下的几乎所有转录本序列信息，目前已被广泛应用于各个领域。而对于水貂这种基因组信息不完善且分子基础研究薄弱的非模式物种来说，利用转录组测序技术比较生长过程序列得到差异表达基因，从而全面且系统地挖掘水貂生长发育调控相关基因，是一个经济快速的途径。本课题组多年来对水貂生长性状数据的监测发现，在45日龄断乳到90日龄期间是水貂的快速生长期[14]，所以本研究对水貂45、90日龄的胸肌肌肉组织进行转录组测序分析，以P < 0.05，|log2FC|>1为条件筛选到279个显著差异基因，对这些基因进行GO和KEGG分析。筛选到与水貂肌肉生长相关的GO条目有66条，从这66条GO条目中，共筛选得到差异表达基因44个，其中上调20个，下调24个，这些差异基因可能在水貂的生长发育过程中起到调控作用，可以作为后续研究水貂生长发育的候选基因。
真核生物转录起始过程需要转录因子的参与，以确保在细胞和机体的整个生命周期中，基因在正确的时间、以正确的数量表达于正确的细胞中。根据其作用特点可以分为普遍转录因子和组织细胞特异性转录因子两类。在水貂快速生长期筛选到的转录因子JUNB和ATF3均属于组织特异性转录因子。JUNB基因最初是在哺乳动物细胞中作为一种直接的早期生长反应基因被发现的。它能够抑制细胞增殖和转化，拮抗c-Jun活性[15]，调控肌细胞蛋白质合成不依赖于mTOR信号通路[16]。ATF3(激活转录因子3)是哺乳动物激活转录因子/cAMP应答元件结合(CREB)蛋白家族的一员。它在绝大多数正常的细胞和健康的组织中呈现稳定的的低表达水平，根据受到的刺激不同表现出促进或抑制基因转录[17]。在本研究中，这两种转录因子均呈下调表达趋势，提示，水貂快速生长有可能是在这两个转录因子的参与调控下，使得大量细胞增殖分化导致的。本研究筛选到的2个转化生长因子GDF11、TGFβ1均属于TGF-β超家族，它们均具有调节细胞生长和分化的功能。GDF11是一种GDF8同源蛋白，在肌肉和神经发育过程中同样的调控机制可能被用于控制组织大小，但是该基因表达是否随年龄增长而增加或减少还没有定论[18-19]。在水貂快速生长期该基因是随着年龄增长呈增加的趋势。TGF-β1它是一种分泌蛋白，执行许多细胞功能，包括控制细胞生长、增殖、分化和凋亡。该基因表达随着水貂生长呈现下调趋势，提示该基因在水貂快速生长中起负调控作用。
肌肉组织是由蛋白质的不断增加及肌细胞的不断增殖分化而形成的，该过程极为复杂，相关基因通过复杂的调控网络与信号通路最终促进肌肉的生长发育。本研究筛选出的与水貂肌肉发育相关基因很多已经被证实，如PPAPDC3(NET39)在肌肉形成和维持方面发挥作用, 表现为负调控肌细胞分化，该基因与Tmem38A和WFS1联合敲除几乎完全阻断了肌管的形成[20-21]，研究报道较少。ABLIM1介导肌动蛋白与细胞质靶点之间的相互作用，具有不同亚型的剪接转录变异，还具有负调控破骨细胞分化的功能[22-23]。CTSZ可以降解肌纤维，还参与生长因子、血管增殖调节等，有助于细胞的迁移、增殖及凋亡等调节，是影响猪生长、肉质等性状的一个重要候选基因[24]。TMSB10在细胞骨架的组织中起重要作用，也被认为是一种新的血管生成、软骨形成和神经生长的刺激因子[25]。TNNT是一种在动物体肌肉和器官组织熟化过程中容易降解的原纤维蛋白，其中TNNT2为心脏肌钙蛋白，它是调节心脏肌肉收缩的肌钙蛋白复合体成分之一，负责使复合体与肌凝蛋白结合，协调肌丝的Ca2+浓度依赖性ATP酶活性，从而在心肌细胞核舒张过程中发挥作用，TNNT2基因在脊椎动物心肌和胚胎骨骼肌中表达[26]。MICAL1属于MICAL家族，该家族参与了多个重要的细胞通路，并受肌动蛋白分离因子的调节，进而调控肌丝的组装和突触链接。MICAL家族是细胞分裂过程中肌动蛋白骨架动力学和膜转运的关键调控因子[27-28]。FLNB基因属于纤丝状肌动蛋白基因家族成员之一，该家族基因能够定义细胞形状，调节细胞间作用力及影响信息传递，这些过程是各个器官正常生长发育所必需的。FLNB间接抑制RhoA通路[29]，与人类的身高有关[30]。
在KEGG通路分析结果中与肌肉发育相关的通路有12条，在这些通路中FOXO信号通路对骨骼肌类型分化具有调控作用。PPARs主要通过调控脂代谢相关基因的表达及调控脂肪细胞的分化来发挥在体脂沉积方面的重要作用。PPARα信号通路主要在肝脏和骨骼肌的脂质代谢中发挥重要作用，可以调节过氧化物酶体的β氧化[31]。PPAR信号通路、p53信号通路、FOXO信号转导通路、氮代谢等信号通路在猪[32]、羊[33]肌肉生长转录组中同样被富集到。
4 结论本研究通过对美国短毛黑水貂45与90日龄的肌肉转录组比对分析，筛选到279个显著差异基因，并对其进行GO功能和KEGG分析，建立了较为全面的水貂快速生长时期转录组数据库，丰富了水貂在基因水平上的数据信息, 为研究水貂生长发育机制及培育大体型水貂品种奠定了基础。
| [1] | LIU N, HE J N, YU W M, et al. Transcriptome analysis of skeletal muscle at prenatal stages in polled Dorset versus small-tailed Han sheep[J]. Genet Mol Res, 2015, 14(1): 1085–1095. DOI: 10.4238/2015.February.6.12 |
| [2] | SILVA-VIGNATO B, COUTINHO L L, CESAR A S M, et al. Comparative muscle transcriptome associated with carcass traits of Nellore cattle[J]. BMC Genomics, 2017, 18: 506. DOI: 10.1186/s12864-017-3897-x |
| [3] | SHANG P, WANG Z X, CHAMBA Y, et al. A comparison of prenatal muscle transcriptome and proteome profiles between pigs with divergent growth phenotypes[J]. J Cell Biochem, 2019, 120(4): 5277–5286. DOI: 10.1002/jcb.27802 |
| [4] | CHAO T L, JI Z B, HOU L, et al. Sheep skeletal muscle transcriptome analysis reveals muscle growth regulatory lncRNAs[J]. PeerJ, 2018, 6: e4619. DOI: 10.7717/peerj.4619 |
| [5] | LI H, WEI X F, YANG J M, et al. Developmental transcriptome profiling of bovine muscle tissue reveals an abundant GosB that regulates myoblast proliferation and apoptosis[J]. Oncotarget, 2017, 8(19): 32083–32100. |
| [6] | ROPKA-MOLIK K, PAWLINA-TYSZKO K, ŻUKOWSK K, et al. Examining the genetic background of porcine muscle growth and development based on transcriptome and miRNAome data[J]. Int J Mol Sci, 2018, 19(4): 1208. DOI: 10.3390/ijms19041208 |
| [7] | CAI Z X, PETERSEN B, SAHANA G, et al. The first draft reference genome of the American mink (Neovison vison)[J]. Sci Rep, 2017, 7: 14564. DOI: 10.1038/s41598-017-15169-z |
| [8] | SONG X C, XU C, LIU Z Y, et al. Comparative transcriptome analysis of mink (Neovison vison) skin reveals the key genes involved in the melanogenesis of black and white coat colour[J]. Sci Rep, 2017, 7: 12461. DOI: 10.1038/s41598-017-12754-0 |
| [9] | CAO X Y, XU C, ZHANG Y F, et al. Comparative transcriptome analysis of embryo invasion in the mink uterus[J]. Placenta, 2019, 75: 16–22. DOI: 10.1016/j.placenta.2018.11.004 |
| [10] | PATEL R K, JAIN M. NGS QC toolkit:a toolkit for quality control of next generation sequencing data[J]. PLoS One, 2012, 7(2): e30619. DOI: 10.1371/journal.pone.0030619 |
| [11] | GRABHERR M G, HAAS B J, YASSOUR M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nat Biotechnol, 2011, 29(7): 644–652. DOI: 10.1038/nbt.1883 |
| [12] | PERTEA G, HUANG X Q, LIANG F, et al. TIGR gene indices clustering tools (TGICL):a software system for fast clustering of large EST datasets[J]. Bioinformatics, 2003, 19(5): 651–652. DOI: 10.1093/bioinformatics/btg034 |
| [13] | ALTSCHUL S F, GISH W, MILLER W, et al. Basic local alignment search tool[J]. J Mol Biol, 1990, 215(3): 403–410. DOI: 10.1016/S0022-2836(05)80360-2 |
| [14] |
荣敏, 涂剑峰, 张志明, 等. 明华黑色水貂生长发育规律与生长曲线拟合研究[J]. 黑龙江畜牧兽医, 2016(9): 227–228, 300.
RONG M, TU J F, ZHANG Z M, et al. Analysis of growth and development rules and growth curve fitting of Minghua black minks[J]. Heilongjiang Animal Science and Veterinary Medicine, 2016(9): 227–228, 300. (in Chinese) |
| [15] | PÉREZ-BENAVENTE B, FARRÀS R. Regulation of GSK3β-FBXW7-JUNB axis[J]. Oncotarget, 2013, 4(7): 956–957. |
| [16] |
罗培.转录因子JunB调控肌细胞蛋白质合成机制的研究[D].武汉: 华中农业大学, 2015.
LUO P.Mechanisms in the regulation of protein synthesis of transcription factor JunB in the myoblasts[D].Wuhan: Huazhong Agricultural University, 2015. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10504-1015387648.htm |
| [17] |
陈博雅, 张聪聪, 李玉琳, 等. 成纤维细胞特异敲除转录因子ATF3基因的小鼠模型构建[J]. 心肺血管病杂志, 2018, 37(5): 469–473.
CHEN B Y, ZHANG C C, LI Y L, et al. Construction of a mouse model of fibroblast cell-specific activating transcriptiun factor 3 knockout[J]. Journal of Cardiovascular and Pulmonary Diseases, 2018, 37(5): 469–473. DOI: 10.3969/j.issn.1007-5062.2018.05.022 (in Chinese) |
| [18] | WALKER R G, POGGIOLI T, KATSIMPARDI L, et al. Biochemistry and biology of GDF11 and myostatin:similarities, differences, and questions for future investigation[J]. Circ Res, 2016, 118(7): 1125–1142. DOI: 10.1161/CIRCRESAHA.116.308391 |
| [19] | HARPER S C, BRACK A, MACDONNELL S, et al. Is growth differentiation factor 11 a realistic therapeutic for aging-dependent muscle defects?[J]. Circ Res, 2016, 118(7): 1143–1150. DOI: 10.1161/CIRCRESAHA.116.307962 |
| [20] | ROBSON M I, DE LAS HERAS J I, CZAPIEWSKI R, et al. Tissue-specific gene repositioning by muscle nuclear membrane proteins enhances repression of critical developmental genes during myogenesis[J]. Mol Cell, 2016, 62(6): 834–847. DOI: 10.1016/j.molcel.2016.04.035 |
| [21] | LIU G H, GUAN T L, DATTA K, et al. Regulation of myoblast differentiation by the nuclear envelope protein NET39[J]. Mol Cell Biol, 2009, 29(21): 5800–5812. DOI: 10.1128/MCB.00684-09 |
| [22] | NARAHARA H, SAKAI E, YAMAGUCHI Y, et al. Actin binding LIM 1 (abLIM1) negatively controls osteoclastogenesis by regulating cell migration and fusion[J]. J Cell Physiol, 2019, 234(1): 486–499. DOI: 10.1002/jcp.26605 |
| [23] |
朱祥.肠道病毒71型(EV71)感染相关基因调控网络建模与分析[D].武汉: 湖北工业大学, 2017.
ZHU X.Regulatory network construction and analysis of Enterovirus 71 infection associated genes[D].Wuhan: Hubei University of Technology, 2017. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10500-1017815430.htm |
| [24] |
崔小荣, 于光辉, 李晓玲, 等. 猪CTSZ基因的克隆及真核表达载体的构建[J]. 青岛农业大学学报:自然科学版, 2017, 34(3): 157–163, 168.
CUI X R, YU G H, LI X L, et al. Cloning of porcine CTSZ gene and construction of its eukaryotic expression vector[J]. Journal of Qingdao Agricultural University:Natural Science, 2017, 34(3): 157–163, 168. (in Chinese) |
| [25] | ZHANG W, CHU W H, LIU Q X, et al. Deer thymosin beta 10 functions as a novel factor for angiogenesis and chondrogenesis during antler growth and regeneration[J]. Stem Cell Res Ther, 2018, 9(1): 166. |
| [26] | SHEN J J, JIN J P. Gene regulation, alternative splicing, and posttranslational modification of troponin subunits in cardiac development and adaptation:a focused review[J]. Front Physiol, 2014, 5: 165. |
| [27] | FRÉMONT S, ROMET-LEMONNE G, HOUDUSSE A, et al. Emerging roles of MICAL family proteins-from actin oxidation to membrane trafficking during cytokinesis[J]. J Cell Sci, 2017, 130(9): 1509–1517. DOI: 10.1242/jcs.202028 |
| [28] | VITALI T, MAFFIOLI E, TEDESCHI G, et al. Properties and catalytic activities of MICAL1, the flavoenzyme involved in cytoskeleton dynamics, and modulation by its CH, LIM and C-terminal domains[J]. Arch Biochem Biophys, 2016, 593: 24–37. DOI: 10.1016/j.abb.2016.01.016 |
| [29] | HU J J, LU J, GOYAL A, et al. Opposing FlnA and FlnB interactions regulate RhoA activation in guiding dynamic actin stress fiber formation and cell spreading[J]. Hum Mol Genet, 2017, 26(7): 1294–1304. DOI: 10.1093/hmg/ddx047 |
| [30] | LEI S F, TAN L J, LIU X G, et al. Genome-wide association study identifies two novel loci containing FLNB and SBF2 genes underlying stature variation[J]. Hum Mol Genet, 2009, 18(9): 1661–1669. DOI: 10.1093/hmg/ddn405 |
| [31] |
秦文.PPARα信号通路基因与牦牛肌内脂肪酸含量的相关性研究[D].兰州: 甘肃农业大学, 2014.
QIN W.Association of PPARα signal pathway genes related to fatty acids composition in Longissimus dorsi muscle of yak[D].Lanzhou: Gansu Agricultural University, 2014. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10733-1015021512.htm |
| [32] |
陈伟.莱芜猪和大白猪背最长肌miRNA与mRNA转录组测序及特征分析[D].泰安: 山东农业大学, 2014.
CHEN W.Sequencing and characterization of the miRNAome and transcriptome of Longissimus dorsi muscle between Laiwu and Yorkshire pigs[D].Tai'an: Shandong Agricultural University, 2014. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10434-1014346464.htm |
| [33] |
孙丽敏.乾华肉用美利奴羊与小尾寒羊肉用性能及肌肉组织mRNA-miRNA表达谱整合分析研究[D].长春: 吉林农业大学, 2017.
SUN L M.Meat performance and integrated analysis of the muscle tissue mRNA-miRNA in Qianhua mutton Merino and Small-tail Han sheep[D].Changchun: Jilin Agricultural University, 2017. (in Chinese) paperuri: (b6a4b25443b498719b9a8248d48574fa) |