NPEPPS基因在强直性脊柱炎发病机制中的功能基因组学分析

童筱君，李辰璐，王国芬，朱小春，王晓冰

1.台州市中心医院（台州学院附属医院） 风湿免疫科，浙江 台州 318000；2.温州医科大学附属第一医院 风湿免疫科，浙江 温州 325015

强直性脊柱炎（ankylosing spondylitis，AS）是一种慢性炎症性自身免疫性疾病，常导致脊柱疼痛、僵硬及终末期关节融合。AS的发病率依种族而异，在中国约为0.3%[1]。AS病因未明，被认为与遗传和环境相关[2-3]，尤其是遗传因素在AS发病中起着重要作用，其遗传度超过90%[4-5]。

人类白细胞抗原（human leukocyte antigen B27，HLA-B27）是迄今为止所知的复杂疾病中相关性最强的一个等位基因[6]。国际AS遗传学联盟发现除了HLA-B27外，还有31个不同的基因位点被确定与AS相关，大部分位于主要组织相容性复合体（major histocompatibility complex，MHC）区域。四种编码位点的氨基肽酶在MHC I-肽复合物抗原递呈通路中发挥作用[7]，可能在AS的发病机制中起重要作用。国内外已有大量关于ERAP1 和ERAP2 两个氨基肽酶编码位点的研究，以及将相关研究结果转化为临床治疗的应用[7-8]。而另一个与疾病相关的氨基肽酶编码位点NPEPPS 在AS中的后续研究尚不足。本项目利用免疫芯片的数据和二代RNA-Seq技术对疾病相关NPEPPS基因的遗传变异进行eQTL分析，探讨NPEPPS在AS发病和进展中的作用。

1 对象和方法

1.1 对象 选取温州医科大学附属第一医院风湿免疫科收治的54例AS患者作为AS组，78例健康对照者为对照组。AS组所有患者均符合1984年纽约分类标准[9]。本研究经温州医科大学附属第一医院伦理委员会批准，所有患者签署了知情同意书。

1.2 免疫芯片基因定型 从入选者5 mL全血中分离出外周血单个核细胞（peripheral blood mononuclear cells，PBMCs），提取DNA，用Illumina平台进行单核苷酸多态性（single nucleotide polymorphism，SNP）微阵列的基因定型检测，基因型被定位到人类基因组NCBI Build 37（hg 19）。用已知的基因型预测未知的基因型并对缺失的数据进行补缺（imputation）和SNP的选择。样本的基因型用PLINK格式（ped和map），提供患者ID和基因组每个SNP的碱基。用PLINK提取包含NPEPPS的染色体区域的SNPs（Chr17，hg 19）。产生每个SNP的次要等位基因频率（minor allele frequencies，MAF），并用MAF＞0.3来进行SNP的剔除。Map文件上的SNP位点被转换到hg 19，去除标志SNPs。利用参考基因组的SNP数据，根据邻近的变异位点间的连锁不平衡（linkage disequilibrium，LD）来推断未定型的SNPs。在补缺之前用SHAPEIT和IMPUTE2进行缺失基因型的补缺，用info＜0.5 来去除补缺较差的SNPs。用PLINK计算免疫芯片的基因型，并设定窗包括所有与疾病相关SNPs（rs9901869）相关性R2＞0.1，以便进行后续的分析。用PLINK抽取这些SNPs，将基因型从等位基因对转换为最小等位基因计数（0、1或2，分别代表主要基因型、杂合基因型和次要基因型）。剩余的MAF＞0.01的SNPs被分成不同的组，每组都与补缺的SNPs密切相关（即组内R2=1），这若干组不同的SNPs将进入后面的分析。

1.3 RNA测序方法 将RNA从PBMCs样本中提取出来，反转录成cDNA，用Illumina TruSeq标准的全RNA样本试剂盒预处理，在Illumina平台进行测定。用TopHat version 2.0.6获得的短read定位到人类基因组NCBI Build 37（hg 19）。对齐的read应用于HTSeq，产生每个基因的read计数，这可以作为经过read多样性校正后的基因表达的测定。来源一个基因的不同转录产物，差异来源于转录起始/终止位点，外显子、内含子及非翻译区的包含与否。用Cufflinks suit来汇编read到融合的个体转录文件中，获得参考转录组。用RSEM来将read排列到参考转录序列中去，来计算同种型的表达丰度。

1.4 eQTL分析方法 计算正态化和模型验证：DESeq2是R语言中一种差异基因表达分析的程序包，用于RNA-Seq计数数据的应用。用内部的函数将计数数据正态化，并观察其大小和水平层次，判断序列的多样性。转录本的计数用DESeq2正态化，并提取转录产物。每个同种型的FASTA序列用Cufflinks的gffread程序获得。用NCBI网站（Madden，2002）的基础局部比对工具（BLAST）来跟人类基因组比对，提供外显子不同或者3’或者5’UTR异常的特征信息。GLME模型用于检验基因型对转录表达的效应。用Wilcoxon符号秩和检验来检测转录表达在AS组和对照组间的差别。基因型和转录产物比率用线性模型进行验证，并对疾病状态进行校正。

1.5 统计学处理方法 采用R v3.6.1 for Linux软件包进行统计学分析。2 组间比较用秩和检验，用Wald test计算DESeq2的P值。P＜0.05为差异有统计学意义。

2 结果

2.1 SNP rs9901869 不同基因型对应的NPEPPS中TCONS_00206150转录本和NPEPPS的表达比较 对2组进行免疫芯片基因分型和RNA测序的NPEPPS表达水平分析，纯合等位基因频率和杂合等位基因频率，2组间NPEPPS及其TCONS_00206150转录本表达的差异无统计学意义（P＞0.05），见图1。

图1 SNP rs9901869不同基因型对应的NPEPPS中TCONS_00206150转录本和NPEPPS的表达水平

2.2 2组NPEPPS 不同的转录本的表达情况 12种NPEPPS 不同的转录本的表达密度见图2；2组不同的转录本的表达水平各不相同，其中2组在TCONS_00206150、TCONS_00206153和TCONS_00206137的表达水平间差异有统计学意义（P＜0.05），见图3。

2.3 SNP rs9901869不同基因型下NPEPSS基因及其转录本表达水平比较 SNP rs9901869 不同基因型所对应的NPEPSS基因及其转录本的表达水平各不相同；其中部分转录本，如EFCAB13、TCONS_00206150、TCONS_00206157、TCONS_00206153、TCONS_00206137和TCONS_00206130 2组的NPEPSS基因及其转录本表达差异有统计学意义（P＜0.05），见图4。

2.4 SNP rs180677251 基因型的NPEPPS 基因及其转录本的表达比较 对主导SNP rs9901869 和附近1 MB范围内的基因的eQTL分析，发现SNP rs180677251不同的基因型所对应的NPEPPS基因及其转录本TCONS_00206137的表达水平差异有统计学意义（P＜0.05），见图5。

3 讨论

图2 12种NPEPPS 不同转录本的表达密度图

图3 2组NPEPPS 不同转录本的表达水平比较

到目前为止，尚未发现与健康对照组相比，AS患者疾病相关的氨基肽酶表达有显著变化的证据[10-11]。而NPEPPS 是编码嘌呤霉素敏感性氨基肽酶的基因。NPEPPS 在大多数组织中以胞质氨基肽酶的形式表达，但在大脑中也发现了一种与膜相关的形式[12-13]，是目前已知的唯一能剪切polyQ序列的胞浆氨基肽酶[14]。鼠模型研究发现神经退行性疾病中该基因的过表达能使蛋白聚集减少，增强自噬作用[15]。而在AS患者的回肠活检中发现自噬的证据，并被证实是IL-23 在肠道表达的驱使者[16]。NPEPPS是否通过自噬作用或者其他机制作用于细胞内肽的处理，仍需要深入探讨。至今为止，AS基因表达研究尚未发现NPEPPS存在基因与基因之间的相互作用[17]。虽然已经发现一个与AS疾病相关的位于NPEPPS基因的关键SNP rs9901869[18]，但是这个SNP对NPEPPS的基因表达是否具有调控作用，以及是否存在其他与AS疾病相关的SNP能影响NPEPPS基因表达，尚不明确。故本项目通过免疫芯片基因分型技术和RNA测序技术，探讨AS患者和健康对照者NPEPPS基因上的SNP rs9901869的基因型与NPEPPS基因表达水平之间的关系，发现AS疾病风险等位基因的基因型与NPEPPS基因及其不同转录本的表达水平相关。进一步揭示了NPEPPS的基因位点多态性影响了NPEPPS的表达水平。鉴于NPEPPS在疾病发病中的潜在作用，推测NPEPPS的eQTL调控可能与AS的发病相关。

图4 不同SNP rs9901869 的基因型下NPEPSS 基因及其转录本的表达比较

图5 SNP rs180677251 基因型的NPEPPS 基因及其转录本的表达比较

将基因型信息和基因表达数据进行配对分析，有助于阐明遗传变异对一个基因转录模式的影响，即表达数量性状基因座。当这些变异也能显示出疾病相关性时，它们将有助于证明单个SNP在疾病发病中的作用。如果这些SNPs干扰了那些被细胞器识别并参与mRNA的产生和DNA基序的加工过程，将会影响附近基因的表达（cis-eQTLs），或其他相互关联通路上的基因表达（trans-eQTLs）[19-20]。随着二代测序技术的进展，RNA测序（RNA-Seq）为eQTL分析提供了一个绝好的平台，使高通量的基因组至转录检测都得以实现。跟传统的基因芯片不同，RNAseq可以提供更多单基因源转录产物的编码序列，而这个差异可能是由于剪接位点不同、外显子的包含与否或3’和5’非翻译区的交替所引起[21]。这些信息可以检测到影响mRNA剪接前的那些变异（sQTLs），而这些变异可能造成同一基因在细胞中表达出完全不同功能的异构体。本研究发现在AS和健康对照者中NPEPSS不同的转录本的表达水平各不相同，其中TCONS_00206150、TCONS_00206153和TCONS_00206137的表达水平差异显著。SNP rs9901869 及其附近的SNP rs180677251的不同基因型所对应的NPEPSS基因及其转录本的表达水平各不相同，进一步证实AS中基因型可能影响基因表达的推论。

本研究的缺陷是进行基因研究的样本数相对较少，虽然这些样本能够检测出大量SNPs，但是仍需要大量样本的研究证实本研究结果。本研究中发现AS和健康对照者中AS相关SNP rs9901869与NPEPPS表达的直接关系仍无法确定，仍需要进一步通过细胞和（或）动物模型进一步验证。

通过对AS和正常人免疫基因的研究，发现除了HLA-B27以外可能还有其他更复杂的遗传学背景和环境因素参与疾病的发生和进展，比如说与免疫炎症性疾病相关的SNP有关，以及还有一些与MHC I-肽复合物抗原递呈通路相关的基因位点变异有关，进一步揭示了AS发病中参与的重要基因位点和关键信号通路，有助于探讨AS发病的潜在机制，为今后的靶向治疗提供了依据。