基于转录组测序分析湖羊波浪形羔皮形成的关键基因

Analysis of key genes in Hu wavy pattern lambskin based on transcriptome sequencing

摘要

旨在对湖羊4/4直毛和4/4波浪形羔皮进行转录组测序，以分析影响羔皮性状的关键基因。采集4/4直毛和4/4波浪形羔皮羊的体侧部皮肤组织（n=4），提取总RNA进行转录组测序，筛选差异表达基因（DEGs），进行基因本体论（GO）和京都基因与基因组百科全书（KEGG）富集分析，并进一步开展KEGG通路的基因集富集分析（GSEA）及DEGs的转录因子预测。对7个DEGs进行实时荧光定量PCR验证（qRT-PCR），同时比较DEGs在不同花案面积个体羔皮直毛和波浪处的表达。结果显示，4/4直毛和4/4波浪形羔皮中共筛选出315个DEGs（阈值|log₂FC|≥1且FDR<0.05），其中波浪形羔皮中274个基因下调，41个基因上调。GO分析发现，DEGs主要富集在细胞外结构组织、细胞外基质和整合素结合等3771个条目；KEGG分析显示，DEGs主要富集在细胞黏附分子、PI3K-AKT信号通路和RAS信号通路等209条通路。GSEA显示，DEGs主要与组氨酸代谢、肾素分泌和蛋白质消化吸收等通路相关；与4/4直毛形羔皮相比，富集到的与毛囊有关的RAS、PI3K-AKT和RAP1信号通路在4/4波浪形羔皮中呈下调趋势；转录因子预测发现编码转录因子（TFs）的差异基因共12个，其中上调基因2个，下调基因10个，主要包括SRY相关高迁移率族盒（SOX）、叉头框（FOX）和T盒（T-box）等TFs家族。qRT-PCR结果显示，筛选出的7个与羔皮性状相关的基因表达趋势与转录组测序结果一致，且IGF1和MFAP5基因在不同花案面积个体羔皮波浪处的表达量显著低于直毛处。综上，筛选出的AEBP1、IGF1、MMP27、MFAP5、RPTN、SPI1和THY1可作为影响羔皮性状的重要候选基因，为进一步解析湖羊波浪形羔皮形成的分子机制奠定基础。

Abstract

The aim of this study was to conduct transcriptome sequencing on Hu lambskins with 4/4 straight hairs and 4/4 wavy pattern， and to identify the key genes affecting lambskin traits. Skin tissues （n=4） were collected from the side of the body of 4/4 straight hair and 4/4 wavy pattern Hu lambs. Total RNA was extracted for transcriptome sequencing， and differentially expressed genes （DEGs） were screened and subjected to gene ontology （GO） and kyoto encyclopedia of genes and genomes （KEGG） enrichment analyses. Gene set enrichment analysis （GSEA） of the KEGG pathways was also conducted， and the DEGs encoding transcription factors （TFs） were identified. Seven DEGs were selected for verification of their transcript profiles by real-time fluorescence quantitative PCR （qRT-PCR）， and the transcript profiles of DEGs in the site of straight hair and wavy pattern in individual lambskin with different wavy areas were compared. Results showed that a total of 315 DEGs （threshold |log₂FC|≥1 and FDR<0.05） were screened out by comparing the transcriptomes of 4/4 straight hair and 4/4 wavy pattern lambskins. Of the 315 DEGs， 274 were down-regulated and 41 were up-regulated in wavy pattern lambskin compared with straight hair lambskin. The GO analysis revealed that DEGs were mainly enriched in 3771 items including extracellular structural tissues， extracellular matrix， and integrin binding. The KEGG analysis showed that DEGs were mainly enriched in 209 pathways including the cell adhesion molecules， PI3K-AKT signaling， and RAS signaling pathways. The GSEA results indicated that DEGs were mainly related to pathways such as histidine metabolism， renin secretion， and protein digestion and absorption. The enriched signaling pathways related to hair follicles， including the RAS， PI3K-AKT， and RAP1 signaling pathways， showed a down-regulated trend in the 4/4 wavy pattern lambskin compared with the 4/4 straight hair lambskin. Among the DEGs， 12 encoded TFs， of which two were up-regulated and 10 were down-regulated in wavy pattern lambskin compared with straight hair lambskin. The down-regulated TF DEGs included those encoding SRY-related high-mobility group box （SOX）， forkhead box （FOX）， and T-box （T-box） TFs. The transcript levels of the seven candidate genes determined by qRT-PCR were consistent with those determined from the transcriptome analysis， and confirmed significantly lower transcript levels of IGF1 and MFAP5 in wavy pattern skin regions than in straight hair skin regions in individuals with different pattern coverage areas. In conclusion， AEBP1， IGF1， MMP27， MFAP5， RPTN， SPI1， and THY1 were identified as important candidate genes affecting lambskin traits. These results lay a foundation for further understanding the molecular mechanism of the formation of wavy lambskin in Hu sheep.

Graphical abstract

关键词

湖羊 / 羔皮 / 毛囊发育 / 羊毛弯曲 / 转录组测序

Key words

Hu sheep / lambskin / hair follicle development / wool crimp / transcriptome sequencing

湖羊（Ovis aries）作为我国地方特色绵羊品种，具有性成熟早、繁殖力高和四季发情等优良性状，其耐粗饲、宜舍饲和适应性强等更是成为占据国内舍饲养殖主导地位不可或缺的优点。湖羊不仅肉质鲜美，而且小湖羊皮（羔皮）更是作为高级衣物材料远销海外。羔皮是出生48 h内的湖羊幼羔宰剥而得，凭借洁白如雪、扑而不散和触感轻柔等特点被誉为“世界软宝石”。湖羊羔皮根据毛纤维形态可分为直毛形和波浪形，其中波浪形羔皮因“S”状弯曲形成的花纹而拥有更高的经济价值^［1］。按花纹宽度可以将波浪形羔皮细分为0.50~1.25 cm的小花、1.25~2.00 cm的中花和>2.00 cm的大花，其中小花最好，中花次之，大花最差^［2］。依据波浪形所占的羔皮面积还可将其分为1/4、2/4、3/4和4/4花案。可根据荐部毛长、花案面积和紧致度等对小湖羊皮进行等级评定。当荐部毛越长、花纹宽且花案不清晰时，等级越低。因此，波浪形比直毛形羔皮品质好，也是湖羊保种选育的方向^［3］。

羊毛的形态变化取决于毛囊的生长发育，毛囊结构包括外根鞘、内根鞘和毛球，由角质形成细胞、毛乳头细胞和黑色素细胞等组成^［4］，位于毛球底部的毛乳头角质形成细胞的分裂和迁移速率影响羊毛的形态。毛囊可分为初级毛囊（primary hair follicle）与次级毛囊（secondary hair follicle），前者结构完整，毛球深达真皮网状层，直径粗大；后者毛球位于真皮浅层，直径细小。研究表明湖羊不同花案羔皮骨形态发生蛋白7（bone morphogenetic protein 7， BMP7）^［5］、基质金属蛋白酶2（matrix metalloproteinase 2， MMP2）^［5］和SRY-box 转录因子18（SRY-box transcription factor 18， SOX18）^［6］的表达量不同，且SOX18可能通过影响毛乳头细胞增殖调控毛囊发育；中卫山羊（Capra hircus）趋化因子配体27（C-C motif chemokine ligand 27， CCL27）、白细胞介素7（interleukin 7， IL7）和Wnt家族成员2（Wnt family member 2， WNT2）基因表达影响羊毛弯曲度^［7］。对中国美利奴羊皮肤组织进行高通量测序，筛选出25个候选miRNAs可能影响羊毛弯曲性状^［8］。小尾寒羊角质形成细胞过表达音猬因子（sonic hedgehog， SHH）基因提高了与毛弯曲相关的角蛋白71（keratin 71， KRT71）、β-连环蛋白（beta-catenin， CTNNB1）和毛透明蛋白（trichohyalin， TCHH）的表达量^［9］。小鼠（Mus musculus）ETS原癌基因2（E26 oncogene 2， Ets2）转录因子功能缺失会影响毛囊形态和排列的发育过程，导致毛发呈现波浪状和触须弯曲^［10］。然而，针对湖羊波浪形羔皮形成的分子机制尚不明晰。

本研究运用转录组测序技术，以直毛和波浪两种不同羔皮性状的湖羊皮肤组织为对象，进行差异表达基因筛选、基因本体论（gene ontology， GO）和京都基因与基因组百科全书（kyoto encyclopedia of genes and genomes， KEGG）等功能富集，并对差异表达基因进行验证，挖掘影响湖羊羔皮性状的关键基因，为深入解析湖羊羔皮性状形成的分子机制以及我国湖羊选育与保种工作提供参考依据。

1 材料与方法

1.1 试验材料

试验于2024年11月在江苏乾宝牧业有限公司进行。选择出生48 h内的羔羊使用皮肤取样器（8.0 mm，美国）采集体侧部皮肤组织。根据羔皮类型选择4/4直毛（S组）和4/4波浪形（W组）羔羊（n=4），采集的皮肤组织生理盐水冲洗后，一部分置于4%多聚甲醛溶液固定，用于苏木素-伊红（hematoxylin-eosin staining， HE）组织学染色分析，另一部分液氮速冻后于-80 ℃冰箱保存，用于转录组测序和实时荧光定量PCR（real-time quantitative PCR， qRT-PCR）分析。同时采集1/4、2/4、3/4花案面积羔皮（n=3）直毛处和波浪处皮肤组织，用于qRT-PCR分析。

1.2 试验方法

1.2.1 HE染色观察

采集的皮肤组织固定24 h后，经70%乙醇→80%乙醇→90%乙醇→100%乙醇Ⅰ→100%乙醇Ⅱ→100%乙醇Ⅲ进行脱水，50%乙醇与50%二甲苯的混合液→二甲苯Ⅰ→二甲苯Ⅱ进行透明，进行石蜡包埋、切片（厚度=5 μm）和染色。将切片置于显微镜（Pannoramic Scan，匈牙利）下扫描拍照，测定初级和次级毛囊的数量和直径，每个样本重复测定3个视野。

1.2.2 总RNA提取

采用FastPure Cell/Tissue Total RNA Isolation Kit V2（RC112-00，Vazyme，南京）依照说明书提取总RNA。使用NanoDrop 2000分光光度计（Thermo Scientific，美国）鉴定RNA纯度和浓度。随后委托武汉影子基因科技有限公司进行后续测序。

1.2.3 转录组文库构建及测序

使用试剂盒VAHTS Universal V8 RNA-seq Library Prep Kit for MGI（NRM605-02， Vazyme，南京）依照说明书构建转录组文库。取500 ng的总RNA，首先使用VAHTS mRNA Capture Beads（N401-02， Vazyme，南京）富集mRNA，打断至300~400 bp。然后合成cDNA第一条链和第二条链，对双链cDNA进行末端修复及3′-末端加A，通过TA连接的方式加上测序接头。接着通过接头引物序列进行PCR扩增。利用KAPA Pure Beads（KS8002，KAPA BIOSYSTEMS）对PCR产物进行纯化回收和片段分选。使用Agilent 2100生物分析仪（2100，美国）检测文库片段大小和浓度，质检合格后使用Illumina novaseq 6000进行测序，测序策略为PE 150。

1.2.4 新转录本鉴定

使用Stringtie（V 2.2.1）软件通过组装和评估这些未被注释的转录本，利用CNCI（V 2.0）和CPC2（V 0.1）软件进行编码潜能鉴定。

1.2.5 转录组数据分析

使用featureCounts（V 2.0.3）对所有样本进行基因水平定量，使用每百万碱基的转录本（transcripts per kilobase million， TPM）作为基因表达水平的指标。采用R包DESeq 2（V 1.36.0，有重复）或edgeR（V 3.38.4，无重复）进行差异分析，以阈值|log₂FC|≥1且FDR<0.05为标准筛选差异基因。基于R包clusterProfiler（V 4.4.4）对差异基因进行GO和KEGG富集分析，并进一步开展KEGG通路的基因集富集分析（gene set enrichment analysis， GSEA）及转录因子预测。

1.2.6 实时荧光定量PCR验证

从差异基因中选取7个与毛囊发育有关的候选基因进行qRT-PCR分析，验证转录组测序结果的准确性。以GAPDH为内参基因，计算目的基因的相对表达量，通过NCBI数据库查询候选基因的编码序列（coding sequences， CDS）并设计基因的引物（南京擎科生物科技有限公司合成，表1）。将提取的RNA反转录成cDNA，采用ChamQ Blue Universal SYBR qRT-PCR Master Mix（Q312-02， Vazyme，南京）10 μL，上下游引物各0.4 μL，cDNA 4 μL进行扩增。每组设3个生物学重复，每个重复检测3次。反应条件：95 ℃ 30 s；95 ℃ 10 s，60 ℃ 30 s，共40个循环。

1.3 数据统计分析

采用IBM SPSS Statistics 26软件中的t检验对所有结果进行统计学分析。运用Graph Pad Prism 10软件进行绘图。P<0.05表示差异显著；P<0.01表示差异极显著。

2 结果与分析

2.1 直毛和波浪形羔皮毛囊比较

4/4直毛和4/4波浪形羔皮组织形态学观察发现，初级毛囊和次级毛囊构成毛囊群（图1）。测定初级毛囊（P）和次级毛囊（S）直径以及次级毛囊与初级毛囊数量比值（S/P），结果如表2所示。4/4直毛羔皮S/P、初级和次级毛囊直径均极显著大于4/4波浪形（P<0.01）。

2.2 转录组测序数据质控

经检验，送测的8个湖羊皮肤组织样本RNA无明显降解和条带污染，符合测序要求。对原始测序数据进行Q20和Q30质控，得到了40044654~46929166 bp的原始测序读数（raw reads）和40044652~46924146 bp的过滤后读数（clean reads）（表3）。8个样本的Q20、Q30值和GC含量分别为99.40%、97.18%和47.44%及以上，质量合格，可进行下一步测序。

2.3 差异表达基因分析

以|log₂FC|≥1，且错误发现率（false discovery rate， FDR）<0.05为筛选条件共得到315个差异表达基因（differentially expressed genes， DEGs），其中274个基因下调，41个基因上调（图2A）。聚类结果显示，湖羊羔皮性状直毛（S）与波浪形（W）两组组内4个样本的基因表达模式高度相似，但组间表达模式差异明显（图2B）。

2.4 DEGs的GO富集分析

对S和W组中的差异表达基因进行GO富集分析，发现共富集在3771个GO条目。其中富集在生物过程（biological process， BP）的条目有3123条，主要包括细胞外结构组织（33个基因）、细胞外基质组织（30个基因）和细胞间黏附（23个基因）。富集在细胞成分（cellular component， CC）的条目有283条，主要包括细胞外基质（33个基因）、含胶原的细胞外基质（17个基因）和胶原三聚物（36个基因）。富集在分子功能（molecular function， MF）的条目有365条，主要包括整合素结合（16个基因）、细胞黏附分子结合（21个基因）和肝素结合（15个基因）。根据P值将富集到的GO功能在BP、CC和MF 3个分类中进行排序，选取TOP 15的GO条目绘制柱状图（图3A）。此外，DEGs还富集在与毛囊发育有关的GO功能，如毛囊形态发生（2个基因）、毛囊发育的调节（1个基因）和毛囊发育（2个基因）。

2.5 DEGs的KEGG通路富集分析

对DEGs进行KEGG通路富集分析，共富集到209条通路。根据P值将Top 20的通路绘制气泡图（图3B），主要包括细胞黏附分子（18个基因）、ECM-receptor交互（10个基因）和磷脂酰肌醇-3-激酶-蛋白激酶B（phosphatidylinositol 3-kinase-protein kinase B， PI3K-AKT）信号通路（17个基因）。此外，还发现一些可能与毛囊发育有关的通路，包括无翅型MMTV整合位点（Wnt）（3个基因）、RAS（7个基因）以及丝裂原活化蛋白激酶（mitogen-activated protein kinase， MAPK）信号通路（5个基因）等。针对PI3K-AKT信号通路进行具体分析，发现共富集到差异基因17个，其中上调基因1个，下调基因16个，主要参与PI3K-AKT信号通路的3种途径（图4）：一种是通过生长因子-受体酪氨酸激酶-胰岛素受体底物1-磷脂酰肌醇3-激酶IA型，一种是通过细胞外基质-整合素α亚基、整合素β亚基、黏着斑激酶-磷脂酰肌醇3-激酶IA型，还有一种是G蛋白βγ亚基-磷脂酰肌醇3-激酶IB型。

2.6 基因集富集分析（GSEA）

对KEGG数据集进行GSEA，显示富集到331条通路，包括组氨酸代谢、肾素分泌和蛋白质消化吸收等。还富集到与毛囊发育有关的如RAS（273个基因，核心基因88个，图5A）、PI3K-AKT（411个基因，核心基因135个，图5B）和RAP1（278个基因，核心基因128个）信号通路，这些通路在波浪形羔皮中呈下调趋势。

2.7 差异表达转录因子分析

转录因子（transcription factor， TF）可促进功能基因的表达和信号传递。在直毛和波浪形羔皮的转录组测序中，发现编码TFs的差异基因共12个，其中上调基因2个，下调基因10个，主要包括Y染色体性别决定区（sex-determining region of the Y chromosome， SRY）相关高迁移率族盒（SRY-related HMG-box， SOX）、叉头框（forkhead box， FOX）和T盒（T-box）等转录因子家族。

2.8 DEGs的实时荧光定量PCR验证

将转录组测序检测到的差异基因与实验室前期全基因组重测序注释的基因进行比较，筛选出7个可能影响羔皮性状的候选基因AEBP1、IGF1、MMP27、MFAP5、RPTN、SPI1和THY1进行实时荧光定量PCR验证（图6A）。结果发现这7个基因在4/4直毛和4/4波浪形的表达模式与转录组测序结果一致，在波浪形羔皮中的表达均下调，表明转录组测序结果数据的可靠性。对其中的2个基因IGF1和MFAP5在不同花案面积羔皮的直毛与波浪处皮肤组织中的表达量进行检测（图6B），结果显示2个基因在直毛处的表达量均显著高于波浪处（P<0.05）。

3 讨论

3.1 毛囊发育与毛发弯曲

毛囊的发育与毛发的弯曲密切相关。本研究4/4直毛羔皮的次级毛囊与初级毛囊数量的比值（S/P）大于4/4波浪形羔皮，与145 d胎龄中国美利奴羊胎羊^［11］和初生内蒙古绒山羊^［12］相近，可能说明次级毛囊主要维持直毛的形成，S/P值的降低导致没有足够的次级毛囊支撑使得毛发弯曲形成波浪形羔皮。研究表明，羊毛弯曲主要源于毛囊有丝分裂细胞增殖的不对称性；此外，还受到毛乳头中角质形成细胞分裂和迁移速率的影响。当细胞分裂和迁移速率过慢时，角质化过程可能会发生在毛乳头内部，造成毛球受力不均，从而导致毛发弯曲^［13］。本研究中4/4直毛的初级毛囊和次级毛囊直径大于4/4波浪形个体，与前人对小花与大花的组织学观察结果一致，大花的初级毛囊和次级毛囊直径均大于小花^［5］，提示直毛由于S/P值高和毛囊直径大能够支撑毛发直挺，而波浪形羊毛由于S/P值低和小直径毛囊导致毛发细软容易弯曲。这与滩羊二毛期羊毛越细卷曲度越大结果类似^［14］。同样，在中国美利奴羊^［15］和新西兰两用羊^［16］中，羊毛平均纤维曲率与平均纤维直径均呈负相关。因此，毛囊直径的大小以及毛囊类型影响湖羊的羊毛弯曲状态，导致波浪形羔皮形成。

3.2 影响湖羊直毛和波浪形羔皮性状的差异基因

本研究315个DEGs中波浪形羔皮个体的下调表达基因为274个，上调表达基因为41个。Sun等^［5］在3对湖羊小花和大花皮肤组织的转录组测序中筛选出在大花皮肤组织上调表达基因分别为100、55和54个，下调表达基因分别为37、82和83个。相比于小花和大花的转录组测序，本研究检测到了更多的DEGs。GO功能富集到毛囊形态发生和发育调节等与毛囊有关的条目。KEGG富集到与毛囊发育有关的通路，包括PI3K-AKT、Wnt以及RAS通路等。Chen等^［17］使用PI3K-AKT特异性抑制剂处理小鼠，结果抑制了表皮干细胞和皮肤衍生前体，影响毛囊再生，表明PI3K-AKT信号通路是毛囊发育的重要通路。本研究发现富集到的PI3K-AKT信号通路的基因主要通过3条途径，且富集到此通路的差异基因大部分下调，表明PI3K-AKT信号通路影响毛囊的发育进而导致羊毛弯曲。Zhou等^［18］研究发现切样同源盒1（cut like homeobox 1， CUX1）可以促进湖羊真皮细胞增殖并影响Wnt/β-catenin信号通路关键基因的表达；Urban等^［19］研究发现RAS信号通路上游和下游元件突变会导致心面皮肤综合征和Costello综合征，患者会出现头发异常，证实RAS信号通路在毛发生长中起到重要作用，本研究结果与上述结果基本一致。因此，湖羊波浪形羔皮的形成可能受到PI3K-AKT、Wnt以及RAS信号通路的影响。

常规的KEGG富集分析依赖于显著上调或下调的基因，容易遗漏部分不显著但具有重要生物学意义的基因。GSEA分析能有效解决这种问题，其不需要明确的阈值，根据整体的趋势更有效地挖掘基因的信息^［20］。本研究对KEGG数据集进行GSEA分析，富集到331条通路，主要包括组氨酸代谢、肾素分泌和蛋白质消化吸收等。此外，富集到的一些与毛囊有关的如RAS、PI3K-AKT和RAP1信号通路呈下调趋势，提示毛囊可能处于功能抑制状态，导致毛囊周期停滞（如退行期或静止期）、细胞增殖减弱，进而影响羔皮性状。

转录因子（TFs）作为蛋白质分子可以调节基因转录，也可以直接影响生物表型^［21］。本研究发现来自SOX、FOX和T-box等多个转录因子家族的差异基因。目前已有研究证明SOX^［22］、FOX^［23］和T-box^［24］等家族中的部分基因与毛囊的生长发育有着密切联系。其中FOX家族在器官的发育中起到重要作用，家族成员叉头框蛋白N1（forkhead box N1， FoxN1）可以促进羊毛弯曲的形成^［23］，叉头框D2反义RNA（forkhead box D2 antisense RNA 1，FOXD2-AS1）也被认为是可以治疗斑秃的靶点^［25］。因此，可以推测这些转录因子可能在湖羊不同羔皮性状的形成中起到重要作用。

3.3 影响湖羊直毛和波浪形羔皮性状的候选基因

结合前期实验室的全基因组重测序结果对DEGs进行筛选，最终得到7个与羔皮性状有关的基因。其中MFAP5和RPTN在湖羊中已有相关报道。研究人员对湖羊小花和直毛的毛乳头细胞进行单细胞测序，发现MFAP5和RPTN在小花中呈下调趋势^［26］，本研究qRT-PCR试验结果与其一致。本研究还发现不同花案面积个体中MFAP5在波浪形羔皮处的表达量低于直毛处。MFAP5富含丝氨酸和苏氨酸，在纤维化和细胞外基质重塑中起到重要作用。Lee等^［27］利用原位空间转录组学研究发现MFAP5在脂膜肌和间质成纤维细胞下缘表达，表明MFAP5可能通过局部的细胞外基质弹性和收缩力导致毛发弯曲；此外，成纤维细胞是组成毛乳头的重要成分，因此MFAP5可能通过影响毛乳头细胞进而影响毛发弯曲导致羔皮性状改变。RPTN作为一种蛋白质编码基因，与角蛋白中间丝的形成有关，在皮肤和头发根髓鞘细胞表面的形成中起到重要作用^［28-29］。在对毛发性状的研究方面，RPTN在欧洲人和东亚人中的突变频率高于非洲人，因此被认为可能促进欧洲人和东南亚人的直发形成^［30］，这也提示RPTN的表达量上升可能促进直毛的形成，本研究中直毛处RPTN表达量高于波浪处的结果与其一致，因此，MFAP5和RPTN表达下调可能促进毛发弯曲进而形成波浪形羔皮。

本研究中筛选出的AEBP1、IGF1、MMP27和SPI1在其他品种的羊中已被证实影响毛囊的生长发育。对藏绒山羊细绒和粗绒蛋白质组学以及皮肤组织转录组学的研究发现，AEBP1是影响羊绒性状的潜在候选基因^［31］。AEBP1是编码羧肽酶A的蛋白家族成员之一，作为转录抑制因子在脂肪生成和平滑肌细胞分化中发挥作用。研究发现，AEBP1相关的Ehlers-Danlos综合征患者存在脱发症状，表明AEBP1可能和毛囊的生长发育存在联系^［32］。IGF1被认为是直接或者间接参与毛囊生长发育的基因^［33］，其在皮肤中的表达量会随着绵羊的生长周期而改变，且存在品种差异性^［34］。IGF1和FGF1间接参与美利奴羊和小尾寒羊毛囊生长发育和纤维性状的调控^［35］。MMP27是基质金属蛋白酶家族的一员，该家族参与胚胎发育、繁殖和组织重塑以及关节炎中细胞外基质的分解等。研究表明，绒山羊中褪黑激素通过下调活化B细胞核因子κB和激活蛋白1转录因子，抑制MMP27等基因表达，同时上调组织金属蛋白酶抑制剂-1/-2/-3的表达，从而延缓皮肤衰老、优化次级毛囊微环境，最终促进次级毛囊活化影响其生长发育^［36］。SPI1是一种蛋白质编码基因，通过编码ETS结构域转录因子在骨髓和B淋巴细胞发育过程中激活基因表达。研究发现lncRNA018392通过募集转录因子SPI1促进集落刺激因子1受体转录，加速绒山羊毛囊真皮成纤维细胞周期，促进增殖并抑制凋亡，而真皮成纤维细胞可以通过调控毛囊干细胞影响毛囊发育和生长^［37-38］。此外，SPI1还在雄激素性脱发中起关键作用^［39］。以上研究表明，AEBP1、IGF1、MMP27和SPI1可能会通过影响毛囊的生长发育进而影响湖羊羔皮性状。

THY1编码位于细胞表面的糖蛋白以及免疫球蛋白超家族成员，不仅参与多种类型的细胞黏附和细胞通讯，而且被广泛用作造血干细胞的标志物。THY1虽然在羊的毛囊发育中未有研究，但对小鼠的毛发研究发现，THY1作为未成熟间充质细胞的标记物，与波形蛋白、CD34和巢蛋白等共同识别出位于玻璃膜外、具有干细胞性质且无增殖活性的A3阳性纺锤形间充质细胞^［40］，这提示THY1可能通过调控毛乳头形成、细胞外基质重塑参与毛囊发育、周期性再生及毛发弯曲的调控。Yes相关蛋白（yes-associated protein， YAP）是调控表皮增殖的关键转录因子，而THY1作为细胞-基质和细胞-细胞互作的重要调节分子，可通过对皮肤稳态控制和再生调控中的YAP活性作用，间接调控毛囊的生长与再生，从而影响毛囊的生长发育^［41］。因此，本研究中的THY1可能通过影响毛囊的生长发育进而作用于湖羊羔皮性状。

4 结论

本研究对4/4直毛和4/4波浪形羔皮进行转录组测序分析，筛选出差异基因315个，其中波浪形羔皮中274个基因下调，41个基因上调。差异基因富集在3771个GO条目和209条KEGG通路，其中PI3K-AKT、Wnt以及RAS信号通路等可能是影响羔皮性状的重要调控途径，SOX、FOX和T-box等是影响羔皮性状的重要转录因子。AEBP1、IGF1、MMP27、MFAP5、RPTN、SPI1、THY1是影响羔皮性状的重要候选基因。研究结果可为进一步解析湖羊波浪形羔皮性状的形成机制奠定基础。

参考文献

[1] Li Z N， Zhang Z L， Du J Z， et al. Morphological investigation on the formation of skin wavy pattern in the Huyang lambs. Scientia Agricultura Sinica， 1982， 15（3）： 73-82.

[2] 李志农， 张仲倫， 杜建中， 湖羊羔皮花纹形成的形态学观察. 中国农业科学， 1982， 15（3）： 73-82.

[3] Chen R H. Identification and grading of Hu lambskin. Journal of Economic Animal， 1980（1）： 40-41.

[4] 程瑞禾. 小湖羊皮的鉴定与分级. 毛皮动物饲养， 1980（1）： 40-41.

[5] Chen F X.Tanning technology of chekiang lamb skin. China Leather， 2005， 34（5）： 40， 46.

[6] 程凤侠. 小湖羊皮鞣制工艺. 中国皮革， 2005， 34（5）： 40， 46.

[7] Prost-Squarcioni C. Histology of skin and hair follicle. Médecine Sciences， 2006， 22（2）： 131-137.

[8] Sun W， Ni R， Yin J F， et al. Genome array of hair follicle genes in lambskin with different patterns. PLoS One， 2013， 8（7）： e68840.

[9] He M L， Lv X Y， Cao X K， et al. SOX18 promotes the proliferation of dermal papilla cells via the Wnt/beta-catenin signaling pathway. International Journal of Molecular Sciences， 2023， 24（23）： 16672.

[10] Li X B， Liu Z F， Liu Y， et al. Mining genes related to wool bending of Zhongwei Goat based on WGCNA and GSEA. Acta Veterinaria et Zootechnica Sinica， 2022， 53（9）： 2930-2943.

[11] 李晓波， 刘占发， 刘悦， 基于WGCNA与GSEA方法挖掘调控中卫山羊羊毛弯曲相关基因. 畜牧兽医学报， 2022， 53（9）： 2930-2943.

[12] Yu L J， Zhang Y H， Lazate·Ainiwaer， et al. Screening microRNAs related to wool crimp in Chinese Merino sheep. Xinjiang Agricultural Sciences， 2021， 58（3）： 573-580.

[13] 于丽娟， 张艳花， 拉扎特·艾尼瓦尔， 中国美利奴羊羊毛弯曲相关microRNA的筛选. 新疆农业科学， 2021， 58（3）： 573-580.

[14] Dong Y J， Hao X J， Wu J Q， et al. Exploration of the effect of SHH on wool bending through Krox20 regulation of IGFBP5 expression based on sheep keratinocytes. Acta Veterinaria et Zootechnica Sinica， 2023， 54（6）： 2365-2375.

[15] 董亚洁， 郝晓静， 吴晋强， 基于绵羊角质形成细胞探究SHH通过Krox20调节IGFBP5的表达影响羊毛弯曲. 畜牧兽医学报， 2023， 54（6）： 2365-2375.

[16] Yamamoto H， Flannery M L， Kupriyanov S， et al. Defective trophoblast function in mice with a targeted mutation of Ets2. Genes & Development， 1998， 12（9）： 1315-1326.

[17] Yang H Y L， Yang Y L， Fang C H， et al. Fetal skin hair follicle development and morphology of Chinese Merino （Xinjiang Junken type）. Southwest China Journal of Agricultural Sciences， 2021， 34（1）： 202-207.

[18] 杨涵羽璐， 杨永林， 方晨辉， 中国美利奴羊（新疆军垦型）胎儿期皮肤毛囊结构及形态学观察. 西南农业学报， 2021， 34（1）： 202-207.

[19] Zhang Y J， Yin J， Li C Q， et al. Study on development of skin and hair follicle from fetal Inner Mongolian Arbas cashmere goats. Acta Veterinaria et Zootechnica Sinica， 2006， 37（8）： 761-768.

[20] 张燕军， 尹俊， 李长青， 内蒙古阿尔巴斯绒山羊胎儿期皮肤毛囊发生发育规律研究. 畜牧兽医学报， 2006， 37（8）： 761-768.

[21] Cloete E， Khumalo N P， Ngoepe M N. The what， why and how of curly hair： a review. Proceedings of the Royal Society A： Mathematical， Physical and Engineering Sciences， 2019， 475（2231）： DOI： 10.1098/rspa.2019.0516.

[22] Li X L， Tao J Z， Ding W， et al. Analysis of wool traits and its relationship in Tan sheep during the er-mao period. China Animal Husbandry & Veterinary Medicine， 2019， 46（2）： 442-448.

[23] 李向龙， 陶金忠， 丁伟， 滩羊二毛期羊毛性状分析及其相关性研究. 中国畜牧兽医， 2019， 46（2）： 442-448.

[24] Yang H， Yang Y L， Liu S R， et al. Correlation analysis of the physical and chemical properties of wool among strains of Chinese Merino （Xinjiang Junken type）. China Animal Husbandry & Veterinary Medicine， 2011， 38（6）： 37-40.

[25] 杨华， 杨永林， 刘守仁， 中国美利奴羊（新疆军垦型）品系间羊毛理化性能及性状相关分析. 中国畜牧兽医， 2011， 38（6）： 37-40.

[26] Pickering N K， Blair H T， Hickson R E， et al. Genetic relationships between dagginess， breech bareness， and wool traits in New Zealand dual-purpose sheep. Journal of Animal Science， 2013， 91（10）： 4578-4588.

[27] Chen Y， Fan Z M， Wang X X， et al. PI3K/Akt signaling pathway is essential for de novo hair follicle regeneration. Stem Cell Research & Therapy， 2020， 11（1）： 144.

[28] Zhou H， Huang S N， Lv X Y， et al. Effect of CUX1 on the proliferation of Hu sheep dermal papilla cells and on the Wnt/β-catenin signaling pathway. Genes， 2023， 14（2）： 423.

[29] Urban J， Qi L， Zhao H， et al. Comparison of hair manifestations in cardio-facio-cutaneous and costello syndromes highlights the influence of the RAS pathway on hair growth. Journal of the European Academy of Dermatology and Venereology， 2020， 34（3）： 601-607.

[30] Wen F， Kumari N， Taylor Vi J G. RNA-seq and gene set enrichment analysis （GSEA） in peripheral blood mononuclear cells （PBMCs）. Methods in Molecular Biology， 2025， 2880： 179-192.

[31] Spitz F， Furlong E E M. Transcription factors： from enhancer binding to developmental control. Nature Reviews Genetics， 2012， 13（9）： 613-626.

[32] Christiano A M. Hair follicle epithelial stem cells get their sox on. Cell Stem Cell， 2008， 3（1）： 3-4.

[33] Yan R Q. Research of FoxN1 on crimping fomation of wool. Taiyuan： Shanxi Agricultural University， 2019.

[34] 闫瑞琴. FoxN1对羊毛弯曲形成作用的研究. 太原： 山西农业大学， 2019.

[35] Grisanti L， Clavel C， Cai X Q， et al. Tbx18 targets dermal condensates for labeling， isolation， and gene ablation during embryonic hair follicle formation. Journal of Investigative Dermatology， 2013， 133（2）： 344-353.

[36] Liu R Q， Liu L D， Xu J D， et al. Identification of potential hub genes in alopecia areata. Experimental Dermatology， 2024， 33（10）： e70002.

[37] Wang S H， Wu T Y， Sun J Y， et al. Single-cell transcriptomics reveals the molecular anatomy of sheep hair follicle heterogeneity and wool curvature. Frontiers in Cell and Developmental Biology， 2021， 9： 800157.

[38] Lee H， Kim S Y， Kwon N， et al. Single-cell and spatial transcriptome analysis of dermal fibroblast development in perinatal mouse skin： dynamic lineage differentiation and key driver genes. Journal of Investigative Dermatology， 2024， 144（6）： 1238-1250.

[39] Krieg P， Schuppler M， Koesters R， et al. Repetin （Rptn）， a new member of the “fused gene” subgroup within the S100 gene family encoding a murine epidermal differentiation protein. Genomics， 1997， 43（3）： 339-348.

[40] Huber M， Siegenthaler G， Mirancea N， et al. Isolation and characterization of human repetin， a member of the fused gene family of the epidermal differentiation complex. Journal of Investigative Dermatology， 2005， 124（5）： 998-1007.

[41] Pospiech E， Lee S D， Kukla-Bartoszek M， et al. Variation in the RPTN gene may facilitate straight hair formation in Europeans and East Asians. Journal of Dermatological Science， 2018， 91（3）： 331-334.

[42] Zhao B R， Wu C L， Sammad A， et al. The fiber diameter traits of Tibetan cashmere goats are governed by the inherent differences in stress， hypoxic， and metabolic adaptations： an integrative study of proteome and transcriptome. BMC Genomics， 2022， 23（1）： 191.

[43] Angwin C， Ghali N， van Dijk F S. Case report： Two individuals with AEBP1-related classical-like EDS： Further clinical characterisation and description of novel AEBP1 variants. Frontiers in Genetics， 2023， 14： 1148224.

[44] Damak S， Su H， Jay N P， et al. Improved wool production in transgenic sheep expressing insulin-like growth factor 1. Nature Biotechnology， 1996， 14（2）： 185-188.

[45] Jia B， Xi J F， Zhang S Y， et al. The developmental patterns of GH-r， IGF-1 and IGF-IR gene expression in sheep skin. Hereditas （Beijing）， 2006， 28（9）： 1078-1082.

[46] Fu J Q， Zhang X Y， Wang D， et al. Analysis of the long non-coding and messenger RNA expression profiles in the skin tissue of super Merino and small-tailed Han sheep. Current Issues in Molecular Biology， 2024， 46（9）： 9588-9606.

[47] Diao X G， Duan C H， Yao L Y， et al. Melatonin promotes the development of secondary hair follicles in adult cashmere goats by activating the Keap1-Nrf2 signaling pathway and inhibiting the inflammatory transcription factors NFkappaB and AP-1. International Journal of Molecular Sciences， 2023， 24（4）： 3403.

[48] Jin M， Fan W Y， Lv S H， et al. LncRNA018392 promotes the proliferation of Liaoning cashmere goat skin fibroblasts by upregulating CSF1R through binding to SPI1. Molecular Biology Reports， 2024， 51（1）： 920.

[49] Liu Y， Guerrero-Juarez C F， Xiao F， et al. Hedgehog signaling reprograms hair follicle niche fibroblasts to a hyper-activated state. Developmental Cell， 2022， 57（14）： 1758-1775.

[50] Premanand A， Rajkumari B R. Bioinformatic analysis of gene expression data reveals Src family protein tyrosine kinases as key players in androgenetic alopecia. Frontiers in Medicine， 2023， 10： 1108358.

[51] Ichikawa C， Izawa T， Juniantito V， et al. Rat hair follicle-constituting cells labeled by a newly-developed somatic stem cell-recognizing antibody： a possible marker of hair follicle development. Histology and Histopathology， 2013， 28（2）： 257-268.

[52] Sedov E， Koren E， Chopra S， et al. THY1-mediated mechanisms converge to drive YAP activation in skin homeostasis and repair. Nature Cell Biology， 2022， 24（7）： 1049-1063.

基金资助

国家重点研发计划(2022YFD1300205)