基于兩樣本孟德爾隨機(jī)化探究甲狀腺功能減退癥與膝骨關(guān)節(jié)炎的因果關(guān)系

徐小東 王配婕 謝文鵬 李琦 傅永升 王衛(wèi)國

【摘 要】目的：通過兩樣本孟德爾隨機(jī)化方法研究甲狀腺功能減退癥與膝骨關(guān)節(jié)炎發(fā)病風(fēng)險(xiǎn)的因果關(guān)系。方法：獲取甲狀腺功能減退癥和膝骨關(guān)節(jié)炎的全基因組關(guān)聯(lián)研究（GWAS）數(shù)據(jù)并進(jìn)行匯總，兩者的數(shù)據(jù)均來自歐洲人群；將單核苷酸多態(tài)性（SNP）作為工具變量（IVs）并選擇敏感的SNPs進(jìn)行分析，通過逆方差加權(quán)法（IVW）、加權(quán)中位數(shù)法（WME）、MR-Egger回歸法進(jìn)行兩樣本孟德爾隨機(jī)化分析，將OR值作為標(biāo)準(zhǔn)評估甲狀腺功能減退癥與膝骨關(guān)節(jié)炎發(fā)病風(fēng)險(xiǎn)的因果關(guān)系，通過異質(zhì)性和多效性檢驗(yàn)對得到的結(jié)果進(jìn)行驗(yàn)證。結(jié)果：共納入27個(gè)SNPs位點(diǎn)作為IVs，IVW法的OR = 4.737，95%CI = ［1.038，20.524］，P = 0.038，表明甲狀腺功能減退癥和膝骨關(guān)節(jié)炎之間有因果關(guān)系。MR-Egger回歸方程截距項(xiàng)為0.001 5（P = 0.775），表明該結(jié)果不受基因多效性的影響。結(jié)論：甲狀腺功能減退癥與膝骨關(guān)節(jié)炎存在因果關(guān)系，且膝骨關(guān)節(jié)炎的發(fā)病風(fēng)險(xiǎn)隨著甲狀腺功能減退癥發(fā)病風(fēng)險(xiǎn)增加而增加。

【關(guān)鍵詞】 膝骨關(guān)節(jié)炎；甲狀腺功能減退癥；孟德爾隨機(jī)化研究；全基因組關(guān)聯(lián)分析；工具變量；單核苷酸多態(tài)性；因果關(guān)系

Exploring the Causal Relationship Between Hypothyroidism and Knee Osteoarthritis Based on a Two Sample

Mendelian Randomization Method

XU Xiao-dong，WANG Pei-jie，XIE Wen-peng，LI Qi，F(xiàn)U Yong-sheng，WANG Wei-guo

【ABSTRACT】Objective：To investigate the causal relationship between hypothyroidism and the risk of knee osteoarthritis using a two sample Mendelian randomization method.Methods：Genome Wide Association Study（GWAS）data on hypothyroidism and knee osteoarthritis were obtained and summarized，both of which were from the European population.Single Nucleotide Polymorphisms（SNPs）were used as instrumental variables（IVs）and sensitive SNPs were selected for analysis.Two sample Mendelian randomization analysis was performed using inverse variance weighting（IVW），weighted median（WME），and MR Egger regression.OR values were used as criteria to evaluate the causal relationship between hypothyroidism and the risk of knee osteoarthritis.The results were validated through heterogeneity and pleiotropy tests.Results：A total of

27 SNPs loci were included as IVs，with an OR = 4.737，95%CI =［1.038，20.524］and P = 0.038 for the IVW method，indicating a causal relationship between hypothyroidism and knee osteoarthritis.The intercept term of the MR-Egger regression equation was 0.001 5（P = 0.775），indicating that the result is not affected by gene pleiotropy.Conclusion：There is a causal relationship between hypothyroidism and knee osteoarthritis，and the risk of knee osteoarthritis increases as the risk of hypothyroidism increases.

【Keywords】 knee osteoarthritis;hypothyroidism;mendelian randomization research;whole genome association analysis;instrumental variables;single nucleotide polymorphism;causal relationship

膝骨關(guān)節(jié)炎（knee osteoarthritis，KOA）使患者生活質(zhì)量下降，社會活動參與能力下降。在我國，本病的發(fā)病率隨著年齡的增長而顯著增加，并且女性發(fā)病率大于男性［1］。在所有膝關(guān)節(jié)疼痛的人群中，KOA患病率接近70%，近25%的KOA患者會遭受嚴(yán)重膝關(guān)節(jié)疼痛的折磨［2-3］。

甲狀腺功能減退癥是比較常見的甲狀腺功能障礙疾病，女性比男性患病率更高［4］。有研究揭示甲狀腺功能障礙疾病影響某些骨科疾病的發(fā)

生［5-8］，抗甲狀腺過氧化物酶抗體指標(biāo)較高的患者膝關(guān)節(jié)軟骨鈣化程度更高［9］，有更大的概率患KOA。除KOA外，甲狀腺功能異常對骨質(zhì)疏松癥等其他退行性骨關(guān)節(jié)疾病存在影響［10-12］。

全基因組關(guān)聯(lián)分析（GWAS）將人類基因組中的單核苷酸多態(tài)性（SNP）作為分子遺傳標(biāo)記物，進(jìn)行全基因組水平上的對照分析或相關(guān)性分析，通過比較發(fā)現(xiàn)影響復(fù)雜性狀的基因變異。孟德爾隨機(jī)化（MR）［13］使用GWAS數(shù)據(jù)分析不同暴露和結(jié)果之間的因果關(guān)系，所以基因和疾病之間的關(guān)聯(lián)不會受到環(huán)境等常見混雜因素的干擾，使效應(yīng)估計(jì)值更接近真實(shí)情況。MR分析有效克服了傳統(tǒng)隨機(jī)對照研究的局限性，是分析甲狀腺功能減退癥與KOA之間因果關(guān)系的可行方法。本研究進(jìn)行兩樣本MR分析，使用GWAS探索甲狀腺功能減退癥與KOA之間的因果關(guān)系。

1 資料與方法

1.1 研究設(shè)計(jì) 本研究使用兩樣本MR分析甲狀腺功能減退癥與KOA的因果關(guān)系，將甲狀腺功能減退癥作為暴露因素，利用與甲狀腺功能減退癥相關(guān)的SNPs作為工具變量（IVs），將KOA作為結(jié)局變量，分別進(jìn)行異質(zhì)性檢驗(yàn)、多效性分析、留一法分析，以確保結(jié)果的可靠性。MR分析主要基于3個(gè)假設(shè)（見圖1）［14］：IVs與暴露之間有很強(qiáng)的相關(guān)性，所選的IVs與混雜因素相互獨(dú)立，IVs只能通過暴露影響結(jié)局（無水平多效性）。在本研究中，即甲狀腺功能減退癥相關(guān)的SNPs與KOA之間相關(guān)性強(qiáng)，與混雜因素相互獨(dú)立，且只能通過甲狀腺功能減退癥影響KOA的發(fā)生。

1.2 資料來源 通過IEU Open GWAS項(xiàng)目（https：//gwas.mrcieu.ac.uk/）搜索獲得甲狀腺功能減退癥和KOA的GWAS數(shù)據(jù)［15］。甲狀腺功能減退癥的數(shù)據(jù)來源于英國布里斯托爾大學(xué)的MRC綜合流行病學(xué)研究所的IEU GWAS數(shù)據(jù)庫，共包含463 010個(gè)樣本量，有9 851 867個(gè)SNPs；KOA數(shù)據(jù)由Tachmazidou I為作者收集，包含403 124個(gè)樣本量，有29 999 696個(gè)SNPs。以上數(shù)據(jù)均來自歐洲人群，見表1。

1.3 工具變量 為了選擇滿足MR 3個(gè)假設(shè)的SNPs，參照千人全基因組，執(zhí)行如下操作。篩選與甲狀腺功能減退癥差異有統(tǒng)計(jì)學(xué)意義（P < 0.05）的SNPs，為了排除SNPs之間的連鎖不平衡（LD），對基因組數(shù)據(jù)進(jìn)行處理，將R2閾值定為0.001，千堿基對設(shè)置為10 000 kb。為了消除工具變量之間的水平多效性，使用工具［16-17］（http：//www.phosanner.medschl.cam.ac.uk/）查找受SNPs強(qiáng)影響的表型，并且剔除與結(jié)局相關(guān)的SNPs。得到與甲狀腺功能減退癥明顯相關(guān)的有效SNPs作為IVs。

1.4 MR分析 本研究中的雙樣本MR數(shù)據(jù)分析基于R軟件中的TwoSampleMR程序包。首先使用逆方差加權(quán)［18］（IVW）作為主要方法初步評估甲狀腺功能減退癥與KOA之間的因果關(guān)系，這個(gè)方法的可靠性取決于IVs的多效性。此外，還選擇了MR-Egger回歸和加權(quán)中位數(shù)（WME）作為補(bǔ)充方法［19-20］，如果所有納入的SNP都符合有效的IVs假設(shè)，那么IVW可以被認(rèn)為是最可靠的結(jié)果［21］。所有結(jié)果均以O(shè)R值和95%置信區(qū)間（95%CI）表

示，暴露與結(jié)果表型之間的相關(guān)性被認(rèn)為差異有統(tǒng)計(jì)學(xué)意義（P < 0.05）。統(tǒng)計(jì)強(qiáng)度F值設(shè)定為 > 10，計(jì)算公式為F = β2exposure = SE2exposure［22］。

1.5 敏感性分析 MR結(jié)果需要利用敏感性分析檢測結(jié)論是否穩(wěn)健。由于兩樣本MR的暴露和結(jié)局分別來自不同的樣本，可能存在人群不同等異質(zhì)性，因此需要通過異質(zhì)性檢驗(yàn)說明。異質(zhì)性檢驗(yàn)?zāi)康氖菣z驗(yàn)每一個(gè)IVs之間有無差異，不同IVs之間的差異越大代表這些IVs的異質(zhì)性就大，P > 0.05表示不存在異質(zhì)性；多效性檢驗(yàn)［23］主要用于檢驗(yàn)納入研究的IVs之間是否存在水平多效性，常用MR-Egger法的截距項(xiàng)表示，如果該值無限接近于0，說明選擇的IVs之間不存在水平多效性；逐個(gè)剔除檢驗(yàn)法［24］是指每次分析剔除一個(gè)IVs，計(jì)算剩余IVs的MR結(jié)果，觀察效應(yīng)值是否發(fā)生較大改變，從而估計(jì)單個(gè)SNP對因果關(guān)聯(lián)的影響。如果剔除某個(gè)IV后，其他IVs估計(jì)出來的MR結(jié)果和總結(jié)果差異很大，說明MR結(jié)果對該IVs敏感。

2 結(jié) 果

2.1 工具變量 經(jīng)過篩選之后，排除有連鎖不平衡的IVs，被選為與甲狀腺功能減退有效IVs的SNPs為27個(gè)，見表2。所有用作IVs的SNP的F值均> 10（30.418～290.893），表明本研究所包括的IVs均有效，受到弱工具變量偏倚影響較小。

2.2 甲狀腺功能減退癥與KOA的MR分析 使用IVW法、MR-Egger法和WME法評估甲狀腺功能減退癥和KOA風(fēng)險(xiǎn)之間是否存在因果關(guān)系。根據(jù)IVW結(jié)果，研究發(fā)現(xiàn)甲狀腺功能減退癥與KOA呈正相關(guān)，甲狀腺功能減退癥可增加KOA的風(fēng)險(xiǎn)（OR = 4.737，95%CI = ［1.038，20.524］，P = 0.038），見表3、圖2；MR-Egger法和加權(quán)中位數(shù)法比較，差異無統(tǒng)計(jì)學(xué)意義（P > 0.05）。但3種方法得到的效應(yīng)方向均為甲狀腺功能減退癥可增加KOA的發(fā)病風(fēng)險(xiǎn)，見圖3。

2.3 敏感性分析 異質(zhì)性檢驗(yàn)分別采用IVW法的q檢驗(yàn)，P = 0.828 ﹥ 0.05，表示SNPs之間不存在異質(zhì)性。多效性分析MR-Egger的egger_intercept截距項(xiàng)為0.001 5，無限接近于0，P = 0.775 ﹥ 0.05，表示不存在多效性。漏斗圖顯示納入的所有SNPs基本對稱，代表因果關(guān)聯(lián)的點(diǎn)對稱分布，說明IVs之間差異較小，表明因果關(guān)聯(lián)受到潛在偏倚影響的可能性較小，見圖4。敏感性分析采用留一法逐一排除各SNPs后，未發(fā)現(xiàn)對因果關(guān)聯(lián)估計(jì)值影響較大的SNPs，見圖5。

3 討 論

KOA主要特征是膝關(guān)節(jié)軟骨丟失，軟骨下骨改變，滑膜炎癥以及半月板變性［25］，影像學(xué)表現(xiàn)為關(guān)節(jié)間隙變窄，邊緣骨贅形成。本病可導(dǎo)致膝關(guān)節(jié)疼痛、晨僵、腫脹、僵硬、畸形。

研究發(fā)現(xiàn)，甲狀腺功能減退癥與各種運(yùn)動系統(tǒng)疾病存在某些聯(lián)系，甲狀腺癌可引起局部性肌張力障礙，在甲狀腺功能減退癥患者服用安定類藥物后，藥物可導(dǎo)致運(yùn)動障礙［26-27］。有學(xué)者發(fā)現(xiàn)，甲狀腺激素的異常信號轉(zhuǎn)導(dǎo)可以增加骨關(guān)節(jié)炎和其他退行性骨關(guān)節(jié)疾病的發(fā)病風(fēng)險(xiǎn)［11-12］。研究表明，膝關(guān)節(jié)周圍軟骨內(nèi)甲狀腺激素受體表達(dá)情況與KOA有較大的聯(lián)系［28］。在增生的骨贅中，尤其是軟骨內(nèi)的成骨部位廣泛分布著甲狀腺α1、變異型甲狀腺α2、甲狀腺β1和甲狀腺β2的mRNA基因。在骨改建部位，超過90%的成骨細(xì)胞表達(dá)甲狀腺α1基因，約75%的細(xì)胞表達(dá)變異型甲狀腺α2基因和甲狀腺β1基因。甲狀腺受體基因的表達(dá)在關(guān)節(jié)骨贅和破骨細(xì)胞中的分布與其受體蛋白的表達(dá)相似；但甲狀腺受體基因在破骨細(xì)胞和軟骨細(xì)胞中高表達(dá)，與其受體蛋白的表達(dá)相反。在KOA患者的膝關(guān)節(jié)軟骨組織中，脫碘酶的存在可以調(diào)節(jié)局部甲狀腺激素水平，該化合物可以降低關(guān)節(jié)軟骨周圍甲狀腺激素的濃度［29］。

大量實(shí)驗(yàn)研究證明，骨骼對甲狀腺激素敏感，甲狀腺激素在骨和關(guān)節(jié)軟骨中的生物學(xué)作用也越來越受到人們的關(guān)注。但是，甲狀腺功能減退癥與KOA之間正向因果關(guān)系背后的具體生物學(xué)機(jī)制仍不清楚。既往研究發(fā)現(xiàn)，甲狀腺激素是骨代謝的重要調(diào)節(jié)劑，對骨重塑至關(guān)重要，甲狀腺激素分泌不足或過量都可能導(dǎo)致骨代謝紊亂［11，30］。

本文采用雙樣本MR研究甲狀腺功能減退癥與KOA之間的因果關(guān)系，證明甲狀腺功能減退癥與KOA之間存在正向因果關(guān)系。甲狀腺功能減退癥會增加患KOA的風(fēng)險(xiǎn)。本研究為預(yù)測甲狀腺功能減退癥患者KOA的發(fā)生、早期預(yù)防和延緩病情進(jìn)展提供了建議，但是，研究結(jié)果還需要通過基礎(chǔ)實(shí)驗(yàn)進(jìn)一步證實(shí)。

參考文獻(xiàn)

［1］ TANG X，WANG S，ZHAN S，et al.The prevalence of symptomatic knee osteoarthritis in China：results from the China health and retirement longitudinal study［J］.Arthritis Rheumatol，2016，68（3）：648-653.

［2］ MICKLE AM，DOMENICO LH，TANNER JJ，et al.

Elucidating factors contributing to disparities in pain-related experiences among adults with or at risk for knee osteoarthritis［J］.Front Pain Res，2023，22（4）：1058476-1058483.

［3］ BARBOUR KE，BORING M，HELMICK CG，et al.Prevalence of severe joint pain among adults with doctor-diagnosed arthritis-united states，2002—2014［J］.MMWR Morb Mortal Wkly Rep，2016，65（39）：1052-1056.

［4］ SCHNEIDER SA，TSCHAIDSE L，REISCH N.Thyroid disorders and movement disorders-a systematic

review［J］.Mov Disord Clin Pract，2023，10（3）：

360-368.

［5］ EK?I M?，ORHUN ?，DEMIR YN，et al.Are serum thyroid hormone，parathormone，calcium，and vitamin D levels associated with lumbar spine degeneration？

A cross-sectional observational clinical study［J］.Eur Spine J，2023，32（5）：1561-1574.

［6］ DENG Y，HAN Y，GAO S，et al.The Physiological functions and polymorphisms of typeⅡ deiodinase［J］.Endocrinol Metab，2023，38（2）：190-202.

［7］ SZULC M，?WI?TKOWSKA-STODULSKA R，PAW?OWSKA E，et al.Vitamin D3 metabolism and its role in temporomandibular joint osteoarthritis and autoimmune thyroid diseases［J］.Int J Mol Sci，2023，24（4）：4080-4100.

［8］ NOH HM，PARK YS，LEE J，et al.A cross-sectional study to examine the correlation between serum TSH levels and the osteoporosis of the lumbar spine in healthy women with normal thyroid function［J］.Osteoporosis Int，2015，26（3）：997-1003.

［9］ TAGOE CE，WANG W，WANG S，et al.Association of anti-thyroid antibodies with radiographic knee osteoarthritis and chondrocalcinosis：a NHANES Ⅲ study［J］.Ther Adv Musculoskelet Dis，2021，13（4）：199-208.

［10］ WILLIAMS GR.Thyroid hormone actions in cartilage and bone［J］.Eur Thyroid J，2013，2（1）：3-13.

［11］ BASSETT JHD，WILLIAMS GR.Role of thyroid hormones in skeletal development and bone mainten-

ance［J］.Endocr Rev，2016，37（2）：135-187.

［12］ NEHLSV.Osteoarthropathies and myopathies associated with disorders of the thyroid endocrine system［J］.Dtsch Med Wochenschr，2018，143（16）：1174-1180.

［13］ BIRNEY E.Mendelian randomization［J］.Cold Spring Harb Perspect Med，2022，12（4）：a041302-a041315.

［14］ DAVIES NM，HOLMES MV，DAVEY SMITH G.Reading Mendelian randomisation studies：a guide，

glossary，and checklist for clinicians［J］.BMJ，2018（362）：k601-k608.

［15］ TACHMAZIDOU I，HATZIKOTOULAS K，SOUTHAM L，et al.Identification of new therapeutic targets for osteoarthritis through genome-wide analyses of UK Biobank data［J］.Nat Genet，2019，51（2）：230-236.

［16］ KAMAT MA，BLACKSHAW JA，YOUNG R，et al.PhenoScanner V2：an expanded tool for searching human genotype-phenotype associations［J］.Bioinformatics，2019，35（22）：4851-4853.

［17］ STALEY JR，BLACKSHAW J，KAMAT MA，et al.PhenoScanner：a database of human genotype-phenotype associations［J］.Bioinformatics，2016，32（20）：3207-3209.

［18］ LAWLOR DA，HARBORD RM，STERNE JAC，et al.

Mendelian randomization：using genes as instruments for making causal inferences in epidemiology［J］.Stat Med，2008，27（8）：1133-1163.

［19］ BOWDEN J，DAVEY SMITH G，HAYCOCK PC，et al.

Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator［J］.Genet Epidemiol，2016，40（4）：304-314.

［20］ BOWDEN J，DEL GRECO MF，MINELLI C，et al.Assessing the suitability of summary data for two-sample Mendelian randomization analyses using MR-Egger regression：the role of the I2 statistic［J］.Int J Epidemiol，2016，45（6）：1961-1974.

［21］ HARTWIG FP，DAVIES NM，HEMANI G，et al.

Two-sample Mendelian randomization：avoiding the down sides of a powerful，widely applicable but potentially fallible technique［J］.Int J Epidemiol，2016，45（6）：1717-1726.

［22］ BURGESS S，THOMPSON SG，CRP CHD GENETICS COLLABORATION.Avoiding bias from weak instruments in Mendelian randomization studies［J］.Int J Epidemiol，2011，40（3）：755-764.

［23］ BOWDEN J，DAVEY SMITH G，BURGESS S.Mendelian randomization with invalid instruments：effect estimation and bias detection through Egger regression［J］.Int J Epidemiol，2015，44（2）：512-525.

［24］ GRONAU QF，WAGENMAKERS EJ.Limitations of bayesian leave-one-out cross-validation for model selection［J］.Comput Brain Behav，2019，2（1）：1-11.

［25］ MADRY H，KON E，CONDELLO V，et al.Early osteoarthritis of the knee［J］.Knee Surg Sports Traumatol Arthrosc，2016，24（6）：1753-1762.

［26］ JAMORA RDG，ESPIRITU AI，SANTIAGO JJ，et al.

Blepharospasm as the presenting feature of papillary thyroid cancer and parathyroid adenoma［J］.J Clin Neurosci，2020，72（1）：460-463.

［27］ SANDYK R.Tardive dyskinesia induced by sulpiride in a patient with hypothyroidism［J］.Clin Neuropharmacol，1986，9（1）：100-101.

［28］ ABU EO，HORNER A，TETI A，et al.The localization of thyroid hormone receptor mRNAs in human bone［J］.Thyroid，2000，10（4）：287-293.

［29］ P?RINGS AS，LOWIN T，DUFNER B，et al.A thyroid hormone network exists in synovial fibroblasts of rheumatoid arthritis and osteoarthritis patients［J］.Sci Rep，2019，9（1）：13235-13248.

［30］ TUCHENDLER D，BOLANOWSKI M.The influence of thyroid dysfunction on bone metabolism［J］.Thyroid Res，2014，7（1）：12-16.

收稿日期：2023-09-05；修回日期：2023-10-10