线粒体核糖体蛋白基因内含子间最佳匹配片段的相对位置分布

The Relative Position Distributions of the Optimal Matched Segments Between Introns of Mitochondrial Ribosomal Protein Genes

摘要

对内含子间相互作用的研究是揭示基因表达调控作用的重要手段。为研究第一内含子与其他内含子间相互作用特征，并考虑内含子间的环状匹配特征包含了这种相互作用的重要信息，以人类、小鼠、红鳍东方鲀、黑腹果蝇和秀丽隐杆线虫五个物种的线粒体核糖体蛋白基因序列为样本，研究第一内含子与其他内含子间的环状匹配特征。采用Smith-Waterman方法，通过局部比对得到每个基因的第一内含子与其他内含子反向互补序列的最佳匹配片段，同时，将这些片段按8-mer CG模体的分类方式分为不同类型的片段。将基因组中所有第一内含子序列的长度标准化处理后，计算最佳匹配片段和各类片段CG对在第一内含子上的相对位置，并分析其物种差异性。结果显示，最佳匹配片段的相对位置分布没有显示出物种差异性，而其CG对的相对位置分布显示了物种 差异性，且1CG类片段的CG对最可几位置随物种进化呈定向进化趋势。研究结果表明，这些最佳匹配片段可能是具有某些组织和功能的生物学元件，而CG模体体现了这些功能元件的重要特征。

Abstract

To investigate the interaction characteristics between the first intron and other introns and interpret the important information of the interaction contained in the circular matching feature between introns， the mitochondrial ribosomal protein gene sequences of five species， including Homo sapiens， Mus musculus， Fugu rubripes， Drosophila melanogaster and Caenorhabditis elegans， were used as samples to study the circular matching features between the first intron and other introns. Firstly， the Smith-Waterman algorithm was used to obtain the optimal matched segments of the reverse complementary sequences between the first intron and other introns of each gene through local alignment， and then these fragments were classified into different types based on the 8-mer CG motif. After standardization of sequence length of all first intron sequences in the genome， the relative position distribution of the optimal matched segments and CG pairs in different types of fragments on the first intron were calculated， and their species differences were analyzed. The results showed that the relative position distribution of the optimal matched segments did not present species differences， while the relative position distribution of CG pairs in different types of fragments presented species differences. Moreover， the most probable positions of CG pairs in the 1CG type fragments exhibited a directed evolutionary trend with species evolution. The result suggested that the optimal matched segments might be some kinds of organized and functional elements， and the CG motif reflected important characteristics of these functional elements.

Graphical abstract

关键词

内含子 / 最佳匹配片段 / CG模体 / 相对位置分布

Key words

intron / optimal matched segments / CG motifs / relative position distribution

内含子作为真核基因组的主要组成部分，在真核生物中普遍存在^［1⁃2］，因此被研究者们广泛关注。研究表明，内含子可以通过剪接来提高mRNA的稳定性、促进mRNA的输出和翻译，进而提高基因的表达^［1］。并且，内含子序列作为一种调控元件，可以调节基因的表达^［2⁃4］。近年来的研究还发现一些内含子可以影响mRNA代谢的多个阶段，包括基因的初始转录、pre-mRNA的编辑以及mRNA的核输出、翻译和衰变^［5］。此外，内含子含有microRNA、snoRNA等多种非编码RNA^［6］，参与多种非编码RNA的功能活动^［4］。

基因表达调控过程通过内含子与内含子或内含子与外显子之间相互作用进行，而内含子与内含子，或内含子与外显子之间最基本的相互作用是碱基配对。例如，正确的密码子反密码子对的形成是保证翻译效率和保真度的关键，外显子环化或内含子环化是由互补的长侧链内含子碱基配对形成的^［7］。侯守庆^［8］研究发现，二斑叶螨sDscamβ基因5′可变外显子盒的选择是通过内含子竞争性碱基配对来调节的。核酸通过分子内部碱基配对的相互作用形成二级结构，二级结构的形成对RNA能够正确发挥功能有重要影响，对分析非编码RNA起到一定帮助。此外，许多研究表明内含子互补匹配片段不仅是环状RNA形成的原因，也是转录或转录后水平上基因表达复杂性和多样性的潜在因素^［9⁃11］。因此，研究内含子的环状匹配问题就显得尤为重要。近年来，第一内含子因为具有存在于转录附近的独特特征，以及表观遗传标记和核小体密度在第一内含子DNA序列上的独特沉积而受到越来越多的关注^{［12⁃15］}，并认识到第一内含子在几种调节基因表达的机制中起关键作用。第一内含子与其他内含子对应的反向互补序列之间的匹配特征可能蕴含了大量重要信息。

除此之外，真核生物基因组DNA序列k-mer非随机适应也有许多研究报道。研究主要关注k-mer字频非随机使用的生物学意义^［16］、k-mer分布的概率模型^［17］、稀有k-mer和富含k-mer的片段及其在基因组序列上的分布^{［18⁃20］}。例如，Bina等^{［21⁃22］}从k-mer使用入手寻找功能位点的调控片段，预测RNA功能片段，给出基因组组装方法。也有研究发现，不同物种基因组序列8-mer具有独立选择定律，即三种CG类（0CG、1CG和2CG）8-mer是各自独立进化的，极大比例DNA序列皆由这三种CG类模体组合而成，并且具有 物种普适性^{［23⁃25］}。基于上述研究，本文认为结合CG模体，研究内含子间环状匹配特征或许可以更好地揭示内含子间相互作用的生物学意义。

本文在前期研究的基础上，利用CG模体的定义将最佳匹配片段分类，研究最佳匹配片段位置分布规律，并统计检验结果的物种差异性，探究内含子序列之间相互作用的特性，进一步挖掘内含子序列的生物学意义。

1 材料与方法

1.1 材料选取

选取了核糖体蛋白基因数据库（robosomal protein gene database，RPG）中所有的线粒体核糖体蛋白基因序列作为研究样本，包括人类（Homo sapiens）、小鼠（Mus musculus）、红鳍东方鲀（Fugu rubripes）、黑腹果蝇（Drosophila melanogaster）和秀丽隐杆线虫（Caenorhabditis elegans）五个物种。选取线粒体核糖体蛋白基因的原因是它们作为一种管家基因在生物学研究中有很多独特的性质，参与了所有蛋白质翻译的关键过程，具有很好的进化 保守性，且内含子长度和数量在所有真核生物中差异不大。通过选择这些保守基因，可以获得内含子之间可靠的功能性的相互作用。关于五个物种的蛋白质基因具体情况见表1。

1.2 比对方法

内含子序列之间最基本的作用方式是碱基配对，配对率越高，相互作用越强。最佳匹配片段是每一个基因的第一内含子与其他内含子反向互补序列配对率最高的片段，该片段最能够反映它们之间的相互作用特点，因此相互作用可用最佳匹配片段进行表征。从样本基因序列中获得内含子序列，除第一内含子之外，将其他内含子转换为反向互补序列。将得到的其他内含子反向互补序列与第一内含子序列通过Smith-Waterman局域比对方法，得到它们之间的最佳匹配片段。

1.3 匹配频率

由于各物种基因组的全部第一内含子序列长度不同，将基因组内的全部第一内含子序列进行标准化处理，得到标准化的第一内含子各碱基的相对位置。具体计算公式为^［26］

n_{ij}=\left\{\begin{array}{l}[l\times N_{ij}/L_i],l\times N_{ij}/L_i\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{e}\mathrm{r},\\ [l\times N_{ij}/L_i]+\mathrm{1},l\times N_{ij}/L_i\mathrm{i}\mathrm{s}\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{g}\mathrm{e}\mathrm{r},\end{array}\right.（1）

其中：n_ij 为标准化后第i条第一内含子的第j个碱基位点；L_i 为第i条第一内含子的长度；l为第一内含子序列标准化的长度（l=100 bp）；N_ij 为标准化前第i条第一内含子的第j个碱基位点。

对每条第一内含子的每个位点定义赋值函数，计算公式为^［26］

（2）

其中：f_ij 为基因组第i条标准化后第一内含子序列中第j位点的赋值函数；n_ia 为最佳匹配片段起始碱基位点，n_ib 为终止碱基位点。对于每一个最佳匹配片段，将与其比对的第一内含子序列转换为包含0和1的数值序列，其中第一内含子与最佳匹配片段配对区内的每个位置赋值为1，配对区外的每个位置赋值为0。最佳匹配片段在每个位点的相对位置频率定义为

F_{\mathrm{r}\mathrm{o}j}=\sum \begin{array}{l}m\\ i=\mathrm{1}\end{array}{f}_{ij}/\sum \begin{array}{l}m\\ i=\mathrm{1}\end{array}\left(n_{ib}-n_{ia}+\mathrm{1}\right)，（3）

其中：F_roj 为标准化第一内含子第j个碱基上最佳匹配片段的相对位置频率；f_ij 为第j个碱基在标准化内含子上的分值。n_ia 和n_ib 为最佳匹配片段的起始和终止碱基位点，m为某一基因中最佳匹配片段的总数。

1.4 最佳匹配片段CG对位点的频率

在高等哺乳动物中，k>6时k⁃mer频谱呈现出稳定的三峰分布^［19］。Chor等^［27］发现不同生物基因组序列的k-mer（k取值在7~11之间）频谱分布呈现单峰或多峰现象。He^［28］认为多细胞动物转录因子结合位点一般为8个碱基长度，并且在高等生物基因组序列中8-mer具有很好的统计意义。另外，有研究发现基因组序列具有8-mer CG独立选择定律^{［23⁃25］}。因此，为进一步揭示最佳匹配片段的位置分布规律，本文利用8-mer CG模体的定义将最佳匹配片段分类。具体方法为，以8 bp为移动窗口，1 bp为移动步长，获得该序列上每个8-mer出现的个数即频次。然后，按照8-mer模体CG二核苷分类法将最佳匹配片段进行分类。将不包含CG对的最佳匹配片段称为0CG类片段，包含一个CG对的称为1CG类片段，8-mer模体集合包含两个或两个以上CG对的称为2CG类片段。

将标准化后的各类最佳匹配片段，运用公式（4）计算得到最佳匹配片段CG对位点的频率。

F_{\mathrm{r}\mathrm{C}\mathrm{G}i}=\frac{n_{\mathrm{m}i}}{{\displaystyle \sum _{\mathrm{1}}^{\mathrm{100}}}{n}_{\mathrm{m}i}}，（4）

其中：F_rCGi 为最佳匹配片段CG对在第一内含子序列第i个位点相对位置频率；n_mi 为最佳匹配片段在第一内含子序列第i个位点出现CG对的个数。

2 结果分析

2.1 最佳匹配片段相对位置分布

为研究最佳匹配片段在第一内含子序列上的相对位置分布，依据公式（1）-（3）计算出标准化第一内含子第j个碱基上最佳匹配片段的相对位置频率（F_roj ），以标准化后的序列位点为横坐标，以F_roj 为纵坐标作图，如图1。

从图1可以看出，对于各物种最佳匹配片段在标准化第一内含子序列的相对位置频率，峰值的个数和峰值的位置均不相同，数据分布也不符合正态分布。因此，为了检验物种间数据的差异性，采用R软件显著性水平0.05的非参数检验方法，结果如图2所示。

从图2可以看出，五个物种最佳匹配片段在第一内含子序列上的相对位置分布两两之间的显著性水平均大于0.05，说明任意两个物种的最佳匹配片段在位置分布方面无差异，即最佳匹配片段在第一内含子序列上的相对位置在物种进化过程中体现出保守性，并未发现物种间的差异性。因此，利用8-mer CG模体的定义将最佳匹配片段分类，通过分析各类片段的特性，尝试发现物种间的差异性。

2.2 最佳匹配片段CG对相对位置分布

统计最佳匹配片段内的CG对在第一内含子上的相对位置分布，并采用R软件进行显著性水平0.05的非参数检验，检验各物种之间的差异。但因为0CG类片段内不存在CG对，因此统计1CG类片段和2CG类片段两种类型。

分别统计出的每个物种所有的最佳匹配片段、1CG类片段及2CG类片段为研究样本，统计其CG对在第一内含子序列上的相对位置频率，并进行五个物种之间的差异性分析，结果如图3所示。

从图3（a）可以看出，人类、小鼠、红鳍东方鲀和秀丽隐杆线虫与黑腹果蝇的差异显著性水平均小于0.000 1，意味着人类、小鼠、红鳍东方鲀、秀丽隐杆线虫分别与黑腹果蝇在最佳匹配片段CG对的相对位置分布呈极显著差异。同时，人类、小鼠分别与秀丽隐杆线虫呈极显著差异，人类和红鳍东方鲀呈显著差异，其余物种两两之间未出现差异。1CG类片段相对位置分布的差异性分析表明（见图3（b）），人类、小鼠、红鳍东方鲀和秀丽隐杆线虫分别与黑腹果蝇呈极显著差异，人类和秀丽隐杆线虫具有差异性。2CG类片段的差异性分析表明（见图3（c）），人类和小鼠分别与秀丽隐杆线虫呈极显著差异，人类和红鳍东方鲀呈极显著差异，小鼠和红旗东方鲀呈显著差异，人类和小鼠、红鳍东方鲀和秀丽隐杆线虫差异不显著。由此可知，最佳匹配片段CG对在第一内含子序列上的相对位置分布具有物种差异性，而且物种间进化距离越大，物种差异越显著。

3 结论与讨论

通过最佳匹配片段在第一内含子序列相对位置频率分布的计算和物种间相对位置分布差异性的分析，发现任意两个物种之间未体现出差异性。该结果说明，每个物种的最佳匹配片段是有组织、有功能的序列。这些功能片段的位置分布表现为物种普适性规律，可能意味着最佳匹配片段位置的信息包含了片段功能的特征。即一些具有特殊位置的最佳匹配片段是具有特殊生物学功能的非编码RNA。另外，朱孝先等^［23］研究发现，内含子和mRNA序列的最佳匹配片段的最可几长度分布与siRNA和miRNA等非编码RNA的长度相同，这意味着最佳匹配片段可能是一种像siRNA和miRNA一样的非编码RNA片段，结果又一次证实了内含子之间的最佳匹配片段也可能是一种具有特殊生物学功能的非编码RNA片段，它们很可能参与基因的表达和调节过程。而且，通过它们之间的相互匹配关系，第一内含子中具有特殊特征的最佳匹配片段，可能通过与其他内含子或外显子的RNA匹配竞争参与调控基因表达。

综合分析最佳匹配片段CG对在第一内含子序列相对位置分布结果，发现不论对总体最佳匹配片段、1CG类片段还是2CG类片段，最佳匹配片段CG对在第一内含子序列上的相对位置分布具有物种差异性，而且结果表明，物种间进化距离越大，物种差异越显著。因此，最佳匹配片段CG对可以很好地体现物种进化过程中的差异性，说明最佳匹配片段CG对中或许包含着与生物进化相关的重要信息。

此外，对比各物种的1CG类最佳匹配片段CG对相对位置分布峰值所对应的位点（如图4），人类在69 bp处、小鼠在64 bp处、红鳍东方鲀在53 bp处、果蝇在30 bp处、秀丽隐杆线虫在25 bp处，分析这些位点的分布发现，随着物种进化程度升高，1CG类片段出现的最可几位置逐渐后移，验证了朱孝先等^［23］提出的8-mer频谱中1CG模体是定向进化的，且与物种的进化显著相关。

本文通过最佳匹配片段在第一内含子序列相对位置分布、最佳匹配片段CG对在第一内含子序列相对位置分布的物种差异性分析发现，最佳匹配片段没有体现出差异性，而最佳匹配片段CG对很好体现了物种在生物进化过程中的差异性。这些规律表明，一些内含子间的最佳匹配片段可能是具有物种普适性的非编码RNA功能片段，而不同进化水平的物种之间的内含子相互作用的复杂性存在很大差异，最佳匹配片段CG对的分布能够体现出其中的部分差异。研究结果为进一步研究非编码RNA的功能提供了一定的理论基础。

参考文献

[1] 王悦冰，郎志宏，黄大昉. 内含子对真核基因表达调控的影响［J］. 生物技术通报， 2008（4）：1-4.

[2] PALMITER R D， SANDGREN E P， AVARBOCK M R， et al. Heterologous introns can enhance expression of transgenes in mice［J］. Proceedings of the National Academy of Sciences of the United States of America， 1991， 88（2）： 478-482.

[3] LI N Q， YANG J， CUI L， et al. Expression of intronic miRNAs and their host gene Igf2 in a murine unilateral ureteral obstruction model［J］. Brazilian Journal of Medical and Biological Research， 2015， 48（6）： 486-492.

[4] ABOU A M， CELLI L， BELOTTI G， et al. GC-AG Introns features in long non-coding and protein-coding genes suggest their role in gene expression regulation［J］. Frontiers in Genetics， 2020 （11）： 488.

[5] BO S L， LI H， ZHANG Q， et al. Potential relations between post-spliced introns and mature mRNAs in the caenorhabditis elegans genome［J］. Journal of Theoretical Biology， 2019，467： 7-14.

[6] MATTICK J S， GAGEN M J. The evolution of controlled multitasked gene networks： The role of introns and other noncoding RNAs in the development of complex organisms［J］. Molecular Biology and Evolution， 2001， 18（9）： 1611-1630.

[7] HAN Z P， CHEN H F， GUO Z H， et al. Circular RNAs and their role in exosomes［J］. Frontiers in Oncology， 2022 （12）： 848341.

[8] HOU S Q， LI G， XU B B， et al. Trans-splicing facilitated by RNA pairing greatly expands sDscam isoform diversity but not homophilic binding specificity［J］. Science Advances， 2022， 8（27）： eabn9458.

[9] ZHANG Y， ZHANG X O， CHEN T， et al. Circular intronic long noncoding RNAs［J］. Molecular Cell， 2013， 51（6）： 792-806.

[10] ZHANG X O， WANG H B， ZHANG Y， et al. Complementary sequence-mediated exon circularization［J］. Cell， 2014， 159（1）： 134-147.

[11] JIAO S H， WU S， HUANG S， et al. Advances in the identification of circular RNAs and research into circRNAs in human diseases［J］. Frontiers in Genetics， 2021 （12）： 665233.

[12] FU L Z， CRAWFORD L， TONG A， et al. Sperm associated antigen 7 is activated by T3 during xenopus tropicalis metamorphosis via a thyroid hormone response element within the first intron［J］. Development Growth & Differentiation， 2022， 64（1）： 48-58.

[13] SINGH O P， MISHRA S， SHARMA G， et al. Evaluation of intron-1 of odorant-binding protein-1 of anopheles stephensi as a marker for the identification of biological forms or putative sibling specie［J］. PLoS One， 2022， 17（7）： e0270760.

[14] SPIJKER H M V， STACKPOLE E E， ALMEIDA S， et al. Ribosome profiling reveals novel regulation of C9ORF72 GGGGCC repeat-containing RNA translation［J］. RNA， 2022， 28（2）： 123-138.

[15] VOSSEBERG J， SCHINKEL M， GREMMEN S， et al. The spread of the first introns in proto-eukaryotic paralogs［J］. Communication Biology， 2022， 5（1）： 476.

[16] CSUROS M， NOE L， KUCHEROV G. Reconsidering the significance of genomic word frequencies［J］. Trends in Genetics， 2007， 23（11）： 543-546.

[17] TULLER T， CHOR B， NELSON N. Forbidden penta-peptides［J］. Protein Science， 2007， 16（10）： 2251-2259.

[18] SUBIRANA J A， MESSEGUER X. The most frequent short sequences in non-coding DNA［J］. Nucleic Acids Research， 2010， 38（4）： 1172-1181.

[19] HARIHARAN R， SIMON R， PILLAI M R， et al. Comparative analysis of DNA word abundances in four yeast genomes using a novel statistical background model［J］. PLoS One， 2013， 8（3）： e58038.

[20] YU H J. Segmented K-mer and its application on similarity analysis of mitochondrial genome sequences［J］. Gene， 2013， 518（2）： 419-424.

[21] BINA M， WYSS P， LAZARUS S A， et al. Discovering sequences with potential regulatory characteristics［J］. Genomics， 2009， 93（4）： 314-322.

[22] BINA M， WYSS P， REN W H， et al. Exploring the characteristics of sequence elements in proximal promoters of human genes［J］. Genomics， 2004， 84（6）： 929-940.

[23] 朱孝先，杨镇，段成妍，等. 基因组序列8-mer频次使用规律及与物种进化的关系［J］. 生物信息学， 2016， 14（4）： 195-202.

[24] YANG Z H， LI H， JIA Y， et al. Intrinsic laws of k-mer spectra of genome sequences and evolution mechanism of genomes［J］. BMC Evolutionary Biology， 2020， 20（1）： 157.

[25] JIA Y， LI H， WANG J F， et al. Spectrum structures and biological functions of 8-mers in the human genome［J］. Genomics， 2019， 111（3）： 483-491.

[26] ZHANG Q， LI H， ZHAO X Q， et al. The evolution mechanism of intron length［J］. Genomics， 2016， 108（2）： 47-55.

[27] CHOR B， HORN D， GOLDMAN N， et al. Genomic DNA k-mer spectra： Models and modalities ［J］. Genome Biology， 2009， 10（10）： R108.

[28] HE X M， TILLO D， VIERSTRA J， et al. Methylated cytosines mutate to transcription factor binding sites that drive tetrapod evolution［J］. Genome Biology and Evolution， 2015， 7（11）： 3155-3169.

基金资助

内蒙古自治区自然科学基金资助项目“mRNA对蛋白质折叠速率影响机制的研究”(2023MS03019)