Chao Wang
Yong Zhu1†
Zhi Liu
Zhe Ruan
Shushan Zhao- 1Department of Orthopaedics, Xiangya Hospital, Central South University, Changsha, China
- 2National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
Backgrounds: Epidemiological studies have repeatedly investigated the association between obesity related anthropometric indicators and body compositions and osteoarthritis (OA). However, the results have remained inconsistent. This work aimed to investigate the genetic correlation and causal associations of obesity related anthropometric indicators and body compositions with knee and hip OA.
Methods: Single-nucleotide polymorphisms associated with the exposures were searched from the recent genome-wide association studies (GWAS) to obtain full statistics. Summary-level results of knee and hip OA were from the UK Biobank and arcOGEN. First, linkage disequilibrium score regression (LD score regression) was applied to detect the genetic correlation (rg). We further performed a series of sensitivity analyses as validation of primary mendelian randomization (MR) results and the specific evidence of potential causal effects was defined.
Results: We found that genetic components in OA had significant correlation with obesity related traits, except waist-to-hip ratio. In the univariable MR analysis, with the exception of waist-to-hip ratio, obesity related anthropometric indicators were causally associated with increased risks of knee and hip OA. For obesity related body compositions, higher fat-free mass in arm, leg, and whole body increased the risk of knee OA but only fat-free mass in leg showed a significant association with hip OA. Meanwhile trunk fat mass and trunk fat percentage, were associated with knee but not with hip OA. Higher fat mass, and fat percentage in arm, leg, and whole body increased the risk of both knee and hip OA. After adjusting for BMI, the multivariable MR showed maintained results in knee OA. However, in hip OA, only fat mass and fat-free mass in arm, leg, trunk and whole body were significantly associated with the risk of hip OA.
Conclusion: The present study suggests genetic evidence for certain causal associations of obesity related anthropometric indicators and body compositions with knee and hip OA, which may provide important insights for the prevention and treatment on OA.
1 Introduction
Osteoarthritis (OA) is a common disabling disease. With the comprehensive impacts of the aging and the increasing of obesity in the global population, an estimated 250 million people around the world are affected by OA, which brings a huge and increasing health burden to the society. Clinically, osteoarthritis most often occurs in the knee, followed by the hand and hip (1).
Genome-wide association study (GWAS) is a research method used to identify genomic variants associated with the risk of a disease or a specific trait. Recent GWAS have identified some statistically independent risk variants, which explain a small proportion of the phenotypic variance and are mainly related to osteoarthritis affecting the knee and hip joints (2). The heritability of osteoarthritis is about 50% and a total of 21 loci, spanning hip and knee OA, have been identified in previous genetic studies (3).
Obesity and osteoarthritis are two interrelated health care problems. It is estimated that by 2030, the increasing weight among population will result in nearly 1.3 billion adults being considered overweight and 573 million classified as obese (4, 5). Epidemiological studies showed that obesity was one of the main risk factors of knee OA (6, 7), and more contradictory results were found in the association of hip OA (8–11). At present, a large number of epidemiological studies explored the relationship between body mass index(BMI), a surrogate marker for adiposity, and OA, but there is relatively little evidence on the association between body compositions (fat mass and fat-free mass) and the risk of OA. What’s more, BMI is not an exclusive measure of adiposity though, but rather reflects the composite of fat, muscle and bone mass (6). Meanwhile, observational studies are commonly susceptible to confounding factors, reverse causality, limited follow-up time and small samples, which may lead to erroneous conclusions.
Although Randomized Controlled Trials (RCTs) are gold standard for causality inference, they are often prohibitively expensive or ethically restrictive (12). Recently, many robust causal inference methods have emerged, which could overcome the limitations of observational research. As one of the most popular genetic correlation analysis methods (13), linkage disequilibrium score regression (LDSC) can evaluate the genetic correlation (rg) between obesity and osteoarthritis, though it is hard to determine the causal relationship between them. Mendelian randomization (MR) analysis, which uses genetic variations as instrument variables (IVs), can be used to investigate the causal relationship between exposure and outcome. Since inheritance of genetic variants is random at conception and is accordingly unaffected by environmental factors or disease process, MR is less susceptible to bias from confounding factors and reverse causality (14). In our study, we used the genome-wide significant single nucleotide polymorphisms (SNPs) retrieved from the large-scale genome-wide association studies as IVs, aiming to investigate the genetic correlation and causal associations of obesity related anthropometric indicators and body compositions with knee and hip OA through LD score regression and two-sample MR analysis. Anthropometric indicators of obesity include: weight, BMI, obesity class (I II III), waist circumference (WC), hip circumference (HC) and waist-to-hip ratio (WHR). Body compositions includes: arm fat mass (AFM), arm fat-free mass (AFFM), arm fat percentage (AFP), leg fat mass (LFM), leg fat-free mass (LFFM), leg fat percentage (LFP), trunk fat mass (TFM), trunk fat-free mass (TFFM), trunk fat percentage (TFP), whole body fat mass (WBFM), whole body fat-free mass (WBFFM) and body fat percentage (BFP).
2 Materials and methods
2.1 Study design
In this study, first, we used the cross-trait LDSC method as implemented in LDHub (15) (https://github.com/bulik/ldsc) to estimate the genetic correlation between anthropometric indicators, body compositions and knee (and/or hip) osteoarthritis. Then we used univariable two-sample MR analysis to evaluate the causal relationship between exposures and outcomes. The exposure variables were divided into anthropometric indicators and body compositions of obesity, while knee OA, hip OA and knee or hip OA were classified as outcome. In order to obtain impartial results, three hypotheses of MR must be met: (a) the genetic IVs should be highly correlated with exposure; (b) there is no confounding of the instrument-outcome association; (c) the genetic IVs affected the outcome only through exposure but not through other biological pathways (16). Therefore, the genome-wide significance SNPs (p ≤5×10-8) were used as instrumental variables, and then a series of sensitivity analyses were carried out, including four MR methods (MR-Egger, weighted median, inverse-variance weighted and weighted mode), as well as heterogeneity and horizontal pleiotropy tests (such as Cochran’s Q test, MR-Egger intercept test and MR-PRESSO test) to verify our results. Bidirectional mendelian randomization analysis and MR Steiger test were performed to estimate the potential reverse causal impact of obesity on knee and hip OA. For the purpose of assessing whether the causal impact of other obesity-related characteristics on the outcomes were mediated by BMI, we adjusted each anthropometric indicator and body composition with BMI for multivariable MR analysis by multivariable inverse-variance weighted method (MV-IVW) (17). Finally, we calculated the proportion of BMI’s effect on knee and hip OA mediated by other obesity-related traits with a mediation analysis (18). The whole workflow of MR analysis were shown in Figure 1.
Figure 1 The whole workflow of MR analysis. MR, Mendelian randomization; GWAS, genome-wide association studies; SNPs, single nucleotide polymorphisms; MAFs, minor allele frequencies; IVW, inverse-variance weighted.
2.2 Data sources
All analyses were based on publicly shared databases and no additional ethical approvals were required. Participants in the genome-wide association study were of European descent on all variables.
Anthropometric indicators of obesity were obtained from the summary statistics of the Genetic Investigation of Anthropometric Traits (GIANT) consortium. A meta-analysis of GWAS using data up to 263407 individuals of European ancestry from 51 studies participating in the GIANT consortium, subdivided indicators into obesity class I (n = 98697), class II (n = 72546) and class III (n = 50364) (19). Summary-level GWAS data for waist to hip ratio (n = 224452), waist circumference (n = 134593), and hip circumference (n = 127997) were acquired from a meta-analysis of 224459 individuals of European ancestry from the GIANT consortium (20). Summary-level data for weight (n = 73137) and BMI (n = 51852) were obtained from previous large-scale GWAS, including male and female participants from European population (21).
Genetic associations with the body compositions of obesity were obtained from summary statistics from the MRC-IEU and Neale Lab. Data of arm fat mass (n = 454757), arm fat-free mass (n = 454753), leg fat mass (n = 454846), leg fat-free mass (n = 454835), leg fat percentage (n = 454826) and body fat percentage (n = 454633) obtained from males and females of European ancestry from the MRC-IEU consortium. While statistics of arm fat percentage (n = 331249), trunk fat mass (n = 331093), trunk fat-free mass (n = 331030), trunk fat percentage (n = 331113), whole body fat mass (n = 330762) and whole-body fat-free mass (n = 331291) were acquired from the Neale Lab consortium (http://www.nealelab.is/blog/2017/7/19/rapid-gwas-of-thousands-of-phenotypes-for-337000-samples-in-the-uk-biobank).
The aggregated data on osteoarthritis was from a large meta-analysis of GWAS with 455221 individuals contained both males and females from the UK Biobank and arcOGEN (2). 24955 patients with knee OA, 15704 patients with hip OA, 39427 patients with knee or hip OA and 378169 control participants from the European population were included in the study.
Some of our exposure data among body compositions of obesity were also from the UK Biobank. Sample overlap may be biased in the direction of the observational association between the risk factor and outcome (22). However, a recent study (23) had shown that for one single large dataset, such as the UK Biobank, main two-sample MR methods can be safely used, except MR-Egger. In addition, we calculated the value of I2 GX (24) of these overlapped sample (Supplementary Table 1). Large I2 GX value (more than 90%) showing a strong genetic variant–exposure associations. In the case of very high variability in instrument strength of different variants, the bias from MR-Egger was greatly reduced (24). Meanwhile, we used a recently developed MR method, MRlap, that is robust to bias introduced by sample overlap, winner’s curse and weak instruments. MRlap revealed almost identical results compared to our primary analysis (Supplementary Table 1), suggesting sample overlap may not substantially bias estimates (25). Finally, we used MR-PRESSO (Mendelian Randomization Pleiotropy RESidual Sum and Outlier) method to detect outliers and MV-IVW for multivariable MR analysis, since they are variations based on the IVW method, it could be safely used (23, 26).
2.3 Genetic instruments selection
We screened the single nucleotide polymorphisms (SNPs) according to following procedures: (1) a genome-wide significance threshold was adopted (p ≤5×10-8); (2) SNPs with minor allele frequencies (MAFs) less than 0.05 were excluded; (3) SNPs with mismatched alleles (i.e., C/T and C/A) were excluded; (4) palindromic variants with equivocal strands (i.e., C/G or A/T) were excluded; (5) In Plink, lead SNPs with parameters (window size= 10000 kb, r2< 0.001) were clumped from the genetic variation selected above.
2.4 Statistical analysis
We calculated the F statistic of SNPs to measure the intensity of IVs. IVs with F statistics < 10 were excluded.
The inverse-variance weighted method (IVW) was used as the main analysis, which is considered as a method to provide the most accurate estimation when all IVs are effective. MR-Egger regression, weighted median and weighted model were used for sensitivity analysis. Cochran’s Q-test was used to test the heterogeneity of IVW estimates. The pleiotropy was detected by MR-Egger intercept test (p for intercept < 0.05). We also performed MR-PRESSO method, including global, outlier and distortion tests, to detect and correct pleiotropy and potential outliers (27). A recently developed MR method was performed to correct the bias introduced by sample overlap (25). Meanwhile, the direction of causality is inferred based on the Steiger test (28).
We defined the evidence of potential causal effects when the following criteria were met: after Bonferroni correction, the MR results of IVW passed the multiple comparisons adjusted p value <0.0025 (0.05/20); The remaining MR methods showed similar size and direction to IVW (p <0.05 in weighted medium and weighted mode, as well as MR-Egger showed a similar effect size); There was no evidence of heterogeneity and level pleiotropy after exclusion of potential outliers. Funnel plots and scatter plots of MR analysis were used to visually assess level pleiotropy and heterogeneity (Supplementary Figure 1-120). MR and other analysis were performed in R (version 4.0.3) using R software packages “two-sample MR”, “Mendelian Randomization”, “MRlap” and “MR-PRESSO”. p values were bilateral, and p value < 0.05 is considered to be suggestively significant for the remaining methods (except IVW).
3 Results
3.1 Genetic correlation and genetic instruments
The characteristics of the SNPs for the exposures and outcomes are showed in Tables 1, 2. First, the GWAS summary data of obesity related parameters and OA were used for genetic correlation analysis, and we found that osteoarthritis genetic components had significant correlation with obesity related traits (rg >0, p ≤0.05), except WHR (Figure 2; Table 3). Then, as genetic instruments we selected related SNPs. Summary statistics of instrumental SNPs as genetic instrumental variables were presented in Table 1 and Supplementary Table 2. The average values of F statistic are shown in Table 1. F statistic < 10 was considered as weak IV and was excluded.
Table 1 Overview of genome-wide association studies used in the analyses for the exposures and characteristics of the Single Nucleotide Polymorphisms (SNPs) considered as instrumental variables.
Table 2 Overview of genome-wide association studies used in the analyses for the outcomes.
Figure 2 Heat map of genetic correlations between knee and hip osteoarthritis and 20 obesity traits grouped in 2 categories from GWAS consortia.
Table 3 Genetic correlation between obesity related anthropometric indicators and body compositions and knee and hip arthritis.
3.2 Univariable mendelian randomization
3.2.1 Causal effects of anthropometric indicators and body compositions of adiposity on knee OA
As shown in Figure 3, in the univariable MR of knee OA, the genetically predicted BMI was significantly correlated with knee OA. For the 1-SD increase in BMI level, the odds ratios were 1.09 (p = 9.09 × 10-8) for knee OA with the IVW method. The remaining MR methods showed consistent size and direction to IVW. Among other anthropometric indicators of obesity, genetic prediction of weight (OR: 1.66; p value: 2.35×10-7), obesity class I (OR: 1.17; p value: 3.57×10-8), obesity class II (OR: 1.15; p value: 9.62×10-15) and obesity class III (OR: 1.11; p value: 3.45×10-6), WC (OR: 1.53; p value: 6.01×10-9) and HC (OR: 1.57; p value: 1.54×10-9) were associated with the risk of knee OA, except WHR (OR: 0.87; p value: 0.39).
Figure 3 Effects of obesity-related exposures on knee OA. Results from univariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on knee OA with IVW, MR-Egger, weighted median, and weighted mode methods.
MRlap were performed since the sample overlap. Supplementary Table 1 showed the preliminary analysis results of genetic prediction levels of body compositions of adiposity (increased 1-SD) and knee OA risk. In the MRlap analysis, genetically raised AFM (OR: 2.17; p value: 2.81×10-65), AFFM (OR: 1.82; p value: 6.91×10-16), AFP (OR: 2.33; p value: 3.77×10-26), LFM (OR: 2.54; p value: 3.03×10-54), LFFM (OR: 1.87; p value: 8.41×10-18), LFP (OR: 3.38; p value: 2.52×10-49), TFM (OR: 1.75; p value: 1.02×10-19), TFP (OR: 1.92; p value: 2.78×10-21), WBFM (OR: 1.92; p value: 1.37×10-34), WBFFM (OR: 1.67; p value: 9.08×10-13) and BFP (OR: 2.39; p value: 1.09×10-37) were associated with increased knee OA risk.
3.2.2 Causal effects of anthropometric indicators and body compositions of adiposity on hip OA
In Figure 4, the genetically predicted BMI was predicted to increase risk of hip OA (OR: 1.12; p value: 7.79×10-9). Consistent results were also observed for other obesity-related measures, such as waist circumference (OR: 1.54; p value: 3.90×10-7) and hip circumference (OR: 1.62; p value: 1.00×10-4).
Figure 4 Effects of obesity-related exposures on hip OA. Results from univariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on hip OA with IVW, MR-Egger, weighted median, and weighted mode methods.
For body compositions of adiposity, AFM (OR: 1.66; p value: 8.39×10-24), AFP (OR: 1.71; p value: 5.51×10-11), LFM (OR: 1.74; p value: 3.62×10-19), LFFM (OR: 1.68; p value: 3.39×10-11), LFP (OR: 1.88; p value: 5.43×10-10), WBFM (OR: 1.57; p value: 1.48×10-14), BFP (OR: 1.58; p value: 7.98×10-8) all showed significant associations with risk of hip OA in MRlap analysis. Higher levels of these factors were associated with increased risk of hip OA.
3.2.3 Causal effects of anthropometric indicators and body compositions of adiposity on knee or hip OA
In Figure 5, the IVW method determined that BMI is associated with knee or hip OA risk (OR 1.10; p = 2.37 × 10-14), suggesting that each 1-SD increment in body mass index was predicted to increase risk of knee or hip OA by 1.10. Except for TFP and WHR, the other exposures showed significant correlation with knee or hip OA risk.
Figure 5 Effects of obesity-related exposures on knee or hip OA. Results from univariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on knee or hip OA with IVW, MR-Egger, weighted median, and weighted mode methods.
3.2.4 Sensitivity analysis
As a further test of the robustness of our findings across MR methodologies, we performed sensitivity analysis and the results were shown in Figures 3–5 and Supplementary Table 2. MRlap revealed almost identical results compared to our primary analysis (Supplementary Table 1). The results from MR sensitivity analysis methods were consistent with those of IVW (p < 0.05 in weighted median and weighted mode, and MR-Egger showed a similar effect size). Although the MR-Egger method in hip circumference in knee OA (OR: 1.68; p value: 6.72×10-2), weight in knee OA (OR: 2.20; p value: 7.44×10-2), hip circumference in hip OA (OR: 2.24; p value: 8.86×10-2) and obesity class 2 in hip OA (OR: 1.25; p value: 5.54×10-2) was not statistically significant, this method is mainly used to detect possible pleiotropic effects in our study.
On knee OA, possible heterogeneity and pleiotropic effects were detected (Supplementary Table 3). Cochran’s Q-test showed that the MR results of AFM, AFFM, AFP, LFM, LFFM, LFP, WBFM, WBFFM, BFP, TFP, TFM, TFFM, weight and WHR were heterogeneous (p <0.05). Except AFFM, WBFFM and TFFM, there was no horizontal pleiotropy as assessed by the MR-Egger intercept test (p intercept > 0.05). Outliers were detected by the MR-PRESSO method, but no significant differences between the casual estimates before and after correction for outliers were found, and p values were >0.05 for distortion tests. On hip OA, these association were maintained after Bonferroni correction for multiple testing (p < 0.05/20 = 2.5 × 10−3). MR-PRESSO identified outliers and provided revised estimates. After removing the abnormal SNPs, the recalculated MR estimates were similar to the results above (Supplementary Table 3). On knee or hip OA, results remained consistent after MR-PRESSO correction for outliers (Supplementary Table 3). Funnel plots and scatter plots of MR analysis were used to visually assess level pleiotropy and heterogeneity (Supplementary Figures 1-120).
3.2.5 Bidirectional mendelian randomization
In the other direction, with genetic liability for knee and hip OA as exposures, we performed MR to explore the causal effect of knee and hip OA on obesity-related traits. After the multiple comparisons adjusted p value (0.05/20), across IVW method, we found no evidence of causal relationships of knee and hip OA with obesity-related traits (Supplementary Table 4). Besides, no evidence of reverse causality across the analyses in the MR Steiger test (Supplementary Table 4).
3.2.6 The Venn diagram
The Venn diagram of the univariable MR was shown in Figure 6 and Supplementary Table 5. The results showed that TFP was an independent risk factor for knee OA. BMI, weight, obesity class (I II III), WC, HC, AFM, AFP, LFM, LFFM, LFP, WBFM and BFP were common risk factors for knee and hip OA.
Figure 6 The Venn diagram of the univariable MR on knee and hip OA. See also Table S5.
3.3 Multivariable mendelian randomization and mediation analysis
Due to the high relevance of the included obesity-related traits with BMI, we further used multivariable MR analysis to adjust the effects of these traits for BMI (Figures 7–9). The conditional F-statistics were shown in Supplementary Table 6. All the F-statistics were >10, implicating absence of weak instrument bias. Multivariable MR demonstrated that other obesity-related traits except WHR causally influenced knee OA even after controlling for BMI, supporting the results of our univariate MR models. In hip OA, after adjusting for BMI, only AFM, AFFM, LFM, LFFM, TFM, TFFM, WBFM, WBFFM were significantly associated with the risk of hip OA. Finally, we further calculated the proportion of BMI’s effect on the outcomes mediated by other obesity-related traits by using mediation analysis. The results showed that the effect of obesity-related traits was partially mediated by BMI (Supplementary Table 7).
Figure 7 Effects of obesity-related exposures on knee OA adjusted for BMI. Results from multivariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on knee OA after adjusting for the effect of BMI, represented by red line.
Figure 8 Effects of obesity-related exposures on hip OA adjusted for BMI. Results from multivariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on hip OA after adjusting for the effect of BMI, represented by red line.
Figure 9 Effects of obesity-related exposures on knee or hip OA adjusted for BMI. Results from multivariable MR analyses showing the effects of genetic predicted obesity related anthropometric indicators and body compositions on knee or hip OA after adjusting for the effect of BMI, represented by red line.
4 Discussion
At present, many previous studies have attributed obesity (defined by BMI) to the pathogenesis of osteoarthritis (29). However, BMI is not an exclusive measure of adiposity though, but rather reflect the composite of fat, muscle and bone mass. Using the recent large-scale GWAS dataset to date, our study performed LD score regression and two-sample MR analysis to determine the genetic correlation and causal associations of obesity related anthropometric indicators and body compositions with knee and hip OA. We investigate whether osteoarthritis genetic components are shared with obesity related traits, and then significant genetic correlation are found. Besides, as shown in Figure 1 and Table 3, we find all the genetic correlations between obesity and knee osteoarthritis are higher than with hip osteoarthritis, such as BMI in knee OA (rg =0.36, p=1.11×10-20) and hip OA (rg=0.21, p=1.08×10-06).
It is well known that osteoarthritis is one of the most common clinical diseases. A number of reviews and clinical studies on the epidemiology of OA have been conducted in the past decades, which had a similar conclusion with our results that obesity has been considered as a risk factor for knee OA (29). However, the relationship between obesity defined by body mass index (BMI) and hip OA remains controversial in many cross-sectional and longitudinal studies. Some previous studies have shown that the association between obesity and hip OA is weak (1, 9, 10). However, recently a prospective cohort study has found that there is an independent association between weight gain and the diagnosis of hip OA (5), which is consistent with our result. In the subsequent studies, WC, HC and WHR were studied separately from BMI.
In addition, conflicting results have emerged regarding the relationship between the body compositions of adiposity (fat mass/muscle mass) and osteoarthritis. Where between fat mass or fat percentage and knee OA or knee cartilage, some studies reported a negative association (30, 31), other studies did not find an association (32). For fat-free mass (FFM), it accounts for a large part of muscle mass and has been shown to be related to both muscle mass and muscle strength. Therefore, FFM is a valuable representative of muscle mass (33). A previous population-based cohort study showed that low fat free mass may interact with knee OA, resulting in an additional reduction in health-related quality of life (HRQOL) in patients with knee OA (34). A cross-sectional study using the data from Korea National Health and Nutrition Examination Survey (KNHANES) showed that the low percentage of lower limb muscle mass was more correlated with knee OA in women (p < 0.05) (35). Another study found that the higher the total extensor muscle cross-sectional area (CSA) and vastus medialis CSA, the greater the loss of patellofemoral cartilage over time (36). In a longitudinal study registered in Sweden, the higher the knee extensor strength of adolescent men, the greater the risk of knee OA in middle age was found (37). While in the cross-sectional study of hip joint muscle strength and joints in subjects with hip OA, better self-reported physical function was found in the greater isometric muscle strength of hip joint and thigh muscle group (38). However, from the results of Venn and MR in our study, that AFFM, TFFM and WBFFM did not meet the definition of potential causal effects, while interestingly, LFFM was found to be significantly associated with both OA of knee and hip.
The difference between our results of MR analysis and epidemiological studies may be due to: (1) Small sample size of the epidemiological study; (2) Confounding factors linked to the chosen exposure and outcome in epidemiological study; (3) Some epidemiological studies only discussed FFM, but not subdivided it into AFFM, TFFM, WBFFM and LFFM; (4) There may be a reverse causal relationship in epidemiological study, for example, low muscle mass may also be a consequence of knee OA, since OA related knee pain leads to muscle abandonment.
4.1 Mechanisms underlying genetic causal links
Some interesting results are also obtained from this MR analysis. The MR result of BMI is consistent with most previous studies, that is, obesity defined by BMI is a risk factor for both knee and hip OA. Regarding the body compositions of adiposity, the results of MR and Venn have showed that AFM, AFP, WBFM, BFP, LFM and LFP are all risk factors for knee and hip OA. Due to a strong correlation (39), AFP and BFP can be used to measure the overall fat content, so as to replicate the results of BMI. TFP can be used to study the effect of visceral fat content on OA, while LFP represents the fat content of legs. These results suggest that the fat mass predicted by genes is highly correlated with the occurrence of OA. Obese adipose tissue seems to create a chronic inflammatory environment. Specific biological mechanisms to explore the link between obesity and OA risk are centered along metabolic and endocrine consequences of body fat accumulation. In basic research, previous studies have attributed the effect of adipose tissue to the adipokines it produces, such as leptin, resistin and adiponectin. Studies on non-weight-bearing arthritis in hands or fingers may support this hypothesis (40). Meanwhile, many studies have proved that adipose tissue can produce inflammatory cytokines such as IL-1, IL-6, and TNF-a. Elevated blood IL-6 levels were thought to be associated with decreased physical function and frailty (41), while higher levels of IL-6 were found to be associated with an increased risk of knee OA progression (42). The serum levels of IL-6 and TNF-a were correlated with articular space stenosis and cartilage loss in the elderly (43). Adipokines and inflammatory cytokines secreted by adipose tissue participate in cartilage degradation, synovitis and bone erosion, leading to OA (44), which is consistent with our study results. However, in this paper, TFP is only related to knee OA, but not hip OA. Whether this discrepancy is related to the different internal structures of knee and hip joint deserves further study. What’s more, as proxy for muscle mass, LFFM is found to be significantly correlated with knee and hip OA, and the correlation still existed after adjusting BMI by multivariable MR analysis. It indicates that genetic predicted LFFM is a risk factor for both knee and hip OA. Recently, skeletal muscle has been identified as an endocrine organ. We believe that some cytokines and inflammatory factors produced by muscle tissue may also lead to OA in knee and hip joints, such as IL-6 (45, 46).
As for the anthropometry indicators of obesity, we find an interesting phenomenon, that is, in univariable MR analysis, with the increase of obesity class (I II III), the odds ratio of both hip OA and knee OA gradually decreased (still greater than 1). In general, abnormal biomechanical loading is considered to be the main mechanism of obesity induced osteoarthritis. According to this result, we consider that, with the gradual increase of BMI, there may be some factors that weaken the effect of biomechanical loading. Of course, this does not exclude chances. Further research could explore this phenomenon. In addition, after the adjustment for BMI with multivariable MR, WC, HC, weight, obesity class (I II III) still have statistical significance in knee OA, while these causal effects disappeared in hip OA. It indicates that the other anthropometric indicators of obesity in our study, except WHR, are strongly correlated with knee OA independent from BMI.
Although the exact pathophysiological mechanism of obesity leading to articular cartilage degeneration is still unclear, mechanical factors may play a partial role in the pathogenesis of OA. With obesity, being overweight will increase the joint load, resulting in a detrimental effect on the weight-bearing joints (47). Maquet et al. (48) found that two to five times a person’s body weight is transmitted across the knee joint. Reilly et al. (49) showed that the reaction force of the entire patellofemoral joint increases from three times the weight when climbing or descending stairs to seven to eight times the weight when squatting. Gait analysis showed that each pound of weight loss could reduce the load on the knee joint by four times (50). Articular cartilage lacks regenerative ability when it bears long-term abnormal mechanical loads, which is prone to degenerative lesions to osteoarthritis. A review of mechanical loading on osteoarthritis showed that mechanical load activated multiple inflammatory pathways, such as IL-1β, TNF-α, inducing chondrocyte apoptosis, synovial inflammation and subchondral bone dysfunction, finally leading to OA (51). In conclusion, obese patients have a higher risk of osteoarthritis related outcomes due to higher mechanical forces exerted on weight-bearing joints under excessive mechanical loading conditions.
4.2 Strengths and limitations
In previous study, Boer et al. (52) found that there was a significant correlation between the genetic components of osteoarthritis and anthropometric traits (BMI, obesity, weight and fat mass) through LD score regression, which was consistent with our research results, but it did not study the relationship between osteoarthritis and other anthropometric traits such as waist circumference, hip circumference or waist-to-hip ratio, nor the genetic correlation between OA and fat mass or fat-free mass in different parts. Previous MR researches have suggested potential associations between BMI and OA (2, 53, 54). With regard to the relationship between other obesity related measurement indicators and OA, recent MR studies (53) indicated that some anthropometric traits such as waist circumference or hip circumference may display modest, but significant, causal associations with OA. However, few studies have adjusted for the effect of BMI through multivariable MR analysis. What is more, no known empirical research has focused on exploring relationships between OA and fat mass, fat free mass or fat mass percentage in different parts. Therefore, we have conducted a more comprehensive investigation.
This study has several strengths. First, this is a more comprehensive study used large-scale GWAS datasets to assess the genetic effect of both anthropometry indicators and body compositions of obesity on knee or hip OA through LD score regression and two-sample MR analysis. Second, given the advantages of MR analysis, it is less likely to be interfered by confounding bias and reverse causality compared with epidemiological studies. Third, we defined strict evidence of potential causal effects, robust methods enabled us to reinforce findings by considering different patterns of pleiotropy and the multivariable method allowed us to distinguish the direct effects of exposures.
However, there are certain limitations in this study. A part of sample overlaps may be biased in the direction of the observational association between risk factors and outcomes. However, a recent study had shown that for one single large dataset, such as the UK Biobank, main two-sample MR methods can be safely used, except MR-Egger, because of its bias (23). Therefore, we define a more rigorous potential causal effects to reduce the possible bias. Meanwhile, we have performed a recently developed MR method, MRlap, to correct the bias introduced by sample overlap. MRlap revealed almost identical results compared to our primary analysis, suggesting sample overlap may not substantially bias estimates. However, as MR-IVW estimates, the corrected results will also be biased if the genetic confounder, through which some of the IVs affect the exposure of interest, is existence. In our study, to avoid the genetic confounder, also called exchangeability or independence, we select the IV with a strict criterion as mentioned in the method section. What’s more, this analytical derivation hinges on a genetic architecture of the exposure, which could reduce the efficiency of the bias correction. However, the effect of the bias will decrease as the sample size increases, and for very large sample size, such as UKBB, the bias is expected to be less pronounced. In addition, since the data in the database includes both male and female, this study does not research the sex difference of obesity on OA. As we all recognized, there are differences in fat, muscle and bone mass between men and women, which may lead to changes in the results. It is also worth studying separately in the future. Finally, the population is restricted to individuals of European ancestry in this study, which limits the transferability to other ethnicities worldwide.
5 Conclusion
Osteoarthritis remains the leading cause of global morbidity and medical costs. Various individual and joint level risk factors are associated with the development and progression of the disease. In this paper, the results of LD score regression have showed the genetic correlation between obesity and osteoarthritis. Mendelian randomization analysis is used to minimize confounding and to avoid reverse causality, and further confirm the effect of genetic predicted anthropometry and body compositions of obesity on knee or hip OA, which may help to better understand the role of BMI, WC, HC, fat and muscle in OA, and finally prevent the occurrence of OA.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
All analyses were based on publicly shared databases and no additional ethical approvals were required.
Author contributions
CW and YZ contributed equally to this work. CW, YZ, and SZ designed the study. SZ and ZL performed the analyses. CW, YZ, ZL, ZR, and SZ analyzed data. CW, YZ, HL, and SZ interpreted the data, wrote the manuscript, and revised the manuscript. All authors approve of the final version of the manuscript.
Funding
This work was financially supported by the National Natural Science Foundation of China (Nos. 81902222, and Nos. 82172399) and the Natural Science Foundation of Hunan Province (2020JJ4928).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.1011896/full#supplementary-material.
Supplementary Table 1 | Results of I2 GX of the overlapped sample and results of MRlap.
Supplementary Table 2 | Extra information about the selected related SNPs.
Supplementary Table 3 | Results of the Cochran’s Q, MR-Egger-intercept and MR-PRESSO tests for detecting horizontal and directional pleiotropy and between SNP-heterogeneity in the univariable Mendelian randomization setting and results of the MR analysis after SNPs screened by PhenoScanner V2.
Supplementary Table 4 | Results of the bidirectional mendelian randomization and MR steiger test.
Supplementary Table 5 | The results of Venn diagram of the univariable MR.
Supplementary Table 6 | Results of the conditional F-statistics in multivariable MR.
Supplementary Table 7 | Results of mediation analysis.
Abbreviations
WC, waist circumference; HC, hip circumference; WHR, waist-to-hip ratio; BMI, body mass index; AFM, arm fat mass; AFFM, arm fat-free mass; AFP, arm fat percentage; LFM, leg fat mass; LFFM, leg fat-free mass; LFP, leg fat percentage; WBFM, whole body fat mass; WBFFM, whole body fat-free mass; BFP, body fat percentage; TFM, trunk fat mass; TFFM, trunk fat-free mass; TFP, trunk fat percentage; OA, osteoarthritis; MR, Mendelian randomization; GWAS, genome-wide association studies; SNPs, single nucleotide polymorphisms; GIANT, Genetic Investigation of Anthropometric Traits; LDSC, linkage disequilibrium score regression; IVs, instrument variables.
References
1. Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet (2019) 393(10182):1745–59. doi: 10.1016/S0140-6736(19)30417-9
2. Tachmazidou I, Hatzikotoulas K, Southam L, Esparza-Gordillo J, Haberland V, Zheng J, et al. Identification of new therapeutic targets for osteoarthritis through genome-wide analyses of UK biobank data. Nat Genet (2019) 51(2):230–6. doi: 10.1038/s41588-018-0327-1
3. Cibrián Uhalte E, Wilkinson JM, Southam L, Zeggini E. Pathways to understanding the genomic aetiology of osteoarthritis. Hum Mol Genet (2017) 26(R2):R193–r201. doi: 10.1093/hmg/ddx302
4. Kelly T, Yang W, Chen CS, Reynolds K, He J. Global burden of obesity in 2005 and projections to 2030. Int J Obes (2005) (2008) 32(9):1431–7. doi: 10.1038/ijo.2008.102
5. Reyes C, Leyland KM, Peat G, Cooper C, Arden NK, Prieto-Alhambra D. Association between overweight and obesity and risk of clinically diagnosed knee, hip, and hand osteoarthritis: A population-based cohort study. Arthritis Rheumatol (Hoboken NJ) (2016) 68(8):1869–75. doi: 10.1002/art.39707
6. Misra D, Fielding RA, Felson DT, Niu J, Brown C, Nevitt M, et al. Risk of knee osteoarthritis with obesity, sarcopenic obesity, and sarcopenia. Arthritis Rheumatol (Hoboken NJ) (2019) 71(2):232–7. doi: 10.1002/art.40692
7. Kulkarni K, Karssiens T, Kumar V, Pandit H. Obesity and osteoarthritis. Maturitas (2016) 89:22–8. doi: 10.1016/j.maturitas.2016.04.006
8. Jiang L, Rong J, Wang Y, Hu F, Bao C, Li X, et al. The relationship between body mass index and hip osteoarthritis: A systematic review and meta-analysis. Joint Bone spine (2011) 78(2):150–5. doi: 10.1016/j.jbspin.2010.04.011
9. Mork PJ, Holtermann A, Nilsen TI. Effect of body mass index and physical exercise on risk of knee and hip osteoarthritis: Longitudinal data from the Norwegian HUNT study. J Epidemiol Community Health (2012) 66(8):678–83. doi: 10.1136/jech-2011-200834
10. Reijman M, Pols HA, Bergink AP, Hazes JM, Belo JN, Lievense AM, et al. Body mass index associated with onset and progression of osteoarthritis of the knee but not of the hip: The Rotterdam study. Ann rheumatic dis (2007) 66(2):158–62. doi: 10.1136/ard.2006.053538
11. Grotle M, Hagen KB, Natvig B, Dahl FA, Kvien TK. Obesity and osteoarthritis in knee, hip and/or hand: An epidemiological study in the general population with 10 years follow-up. BMC musculoskeletal Disord (2008) 9:132. doi: 10.1186/1471-2474-9-132
12. Ma B, Li C, Pan J, Zhang S, Dong H, Wu Y, et al. Causal associations of anthropometric measurements with fracture risk and bone mineral density: A mendelian randomization study. J Bone mineral Res Off J Am Soc Bone Mineral Res (2021) 36(7):1281–7. doi: 10.1002/jbmr.4296
13. Luo P, Cheng S, Zhang F, Feng R, Xu K, Jing W, et al. A large-scale genetic correlation scan between rheumatoid arthritis and human plasma protein. Bone Joint Res (2022) 11(2):134–42. doi: 10.1302/2046-3758.112.BJR-2021-0270.R1
14. Nitsch D, Molokhia M, Smeeth L, DeStavola BL, Whittaker JC, Leon DA. Limits to causal inference based on mendelian randomization: A comparison with randomized controlled trials. Am J Epidemiol (2006) 163(5):397–403. doi: 10.1093/aje/kwj062
15. Bulik-Sullivan BK, Loh PR, Finucane HK, Ripke S, Yang J, Patterson N, et al. LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat Genet (2015) 47(3):291–5. doi: 10.1038/ng.3211
16. Davies NM, Holmes MV, Davey Smith G. Reading mendelian randomisation studies: A guide, glossary, and checklist for clinicians. BMJ (2018) 362:k601. doi: 10.1136/bmj.k601
17. Burgess S, Thompson SG. Multivariable mendelian randomization: The use of pleiotropic genetic variants to estimate causal effects. Am J Epidemiol (2015) 181(4):251–60. doi: 10.1093/aje/kwu283
18. Burgess S, Thompson DJ, Rees JMB, Day FR, Perry JR, Ong KK. Dissecting causal pathways using mendelian randomization with summarized genetic data: Application to age at menarche and risk of breast cancer. Genetics (2017) 207(2):481–7. doi: 10.1534/genetics.117.300191
19. Berndt SI, Gustafsson S, Mägi R, Ganna A, Wheeler E, Feitosa MF, et al. Genome-wide meta-analysis identifies 11 new loci for anthropometric traits and provides insights into genetic architecture. Nat Genet (2013) 45(5):501–12. doi: 10.1038/ng.2606
20. Shungin D, Winkler TW, Croteau-Chonka DC, Ferreira T, Locke AE, Mägi R, et al. New genetic loci link adipose and insulin biology to body fat distribution. Nature (2015) 518(7538):187–96. doi: 10.1038/nature14132
21. Randall JC, Winkler TW, Kutalik Z, Berndt SI, Jackson AU, Monda KL, et al. Sex-stratified genome-wide association studies including 270,000 individuals show sexual dimorphism in genetic loci for anthropometric traits. PloS Genet (2013) 9(6):e1003500. doi: 10.1371/journal.pgen.1003500
22. Burgess S, Davies NM, Thompson SG. Bias due to participant overlap in two-sample mendelian randomization. Genet Epidemiol (2016) 40(7):597–608. doi: 10.1002/gepi.21998
23. Minelli C, Del Greco MF, van der Plaat DA, Bowden J, Sheehan NA, Thompson J. The use of two-sample methods for mendelian randomization analyses on single large datasets. Int J Epidemiol (2021) 50(5):1651–9. doi: 10.1093/ije/dyab084
24. Bowden J, Del Greco MF, Minelli C, Davey Smith G, Sheehan NA, Thompson JR. Assessing the suitability of summary data for two-sample mendelian randomization analyses using MR-egger regression: The role of the I2 statistic. Int J Epidemiol (2016) 45(6):1961–74. doi: 10.1093/ije/dyw220
25. Mounier N, Kutalik Z. Correction for sample overlap, winner’s curse and weak instrument bias in two-sample mendelian randomization. bioRxiv (2021) 2021:03.26.437168. doi: 10.1101/2021.03.26.437168
26. Burgess S, Davey Smith G, Davies NM, Dudbridge F, Gill D, Glymour MM, et al. Guidelines for performing mendelian randomization investigations. Wellcome Open Res (2019) 4:186. doi: 10.12688/wellcomeopenres.15555.1
27. Freuer D, Linseisen J, O'Mara TA, Leitzmann M, Baurecht H, Baumeister SE, et al. Body fat distribution and risk of breast, endometrial, and ovarian cancer: A two-sample mendelian randomization study. Cancers (2021) 13(20)5053:1–12. doi: 10.3390/cancers13205053
28. Hemani G, Tilling K, Davey Smith G. Orienting the causal relationship between imprecisely measured traits using GWAS summary data. PloS Genet (2017) 13(11):e1007081. doi: 10.1097/BOR.0000000000000479
29. Vina ER, Kwoh CK. Epidemiology of osteoarthritis: Literature update. Curr Opin Rheumatol (2018) 30(2):160–7. doi: 10.1097/BOR.0000000000000479
30. Berry PA, Wluka AE, Davies-Tuck ML, Wang Y, Strauss BJ, Dixon JB, et al. The relationship between body composition and structural changes at the knee. Rheumatology (2010) 49(12):2362–9. doi: 10.1093/rheumatology/keq255
31. Teichtahl AJ, Wang Y, Wluka AE, Cicuttini FM. Obesity and knee osteoarthritis: New insights provided by body composition studies. Obes (Silver Spring Md) (2008) 16(2):232–40. doi: 10.1038/oby.2007.30
32. Abbate LM, Stevens J, Schwartz TA, Renner JB, Helmick CG, Jordan JM. Anthropometric measures, body composition, body fat distribution, and knee osteoarthritis in women. Obes (Silver Spring Md) (2006) 14(7):1274–81. doi: 10.1038/oby.2006.145
33. Lafortuna CL, Maffiuletti NA, Agosti F, Sartorio A. Gender variations of body composition, muscle strength and power output in morbid obesity. Int J Obes (2005) (2005) 29(7):833–41. doi: 10.1038/sj.ijo.0802955
34. Visser AW, de Mutsert R, Bloem JL, Reijnierse M, Kazato H, le Cessie S, et al. Do knee osteoarthritis and fat-free mass interact in their impact on health-related quality of life in men? Results from a population-based cohort. Arthritis Care Res (2015) 67(7):981–8. doi: 10.1002/acr.22550
35. Suh DH, Han KD, Hong JY, Park JH, Bae JH, Moon YW, et al. Body composition is more closely related to the development of knee osteoarthritis in women than men: A cross-sectional study using the fifth Korea national health and nutrition examination survey (KNHANES V-1, 2). Osteoarthritis cartilage (2016) 24(4):605–11. doi: 10.1016/j.joca.2015.10.011
36. Goldman LH, Tang K, Facchetti L, Heilmeier U, Joseph GB, Nevitt MC, et al. Role of thigh muscle cross-sectional area and strength in progression of knee cartilage degeneration over 48 months - data from the osteoarthritis initiative. Osteoarthritis cartilage (2016) 24(12):2082–91. doi: 10.1016/j.joca.2016.07.004
37. Turkiewicz A, Timpka S, Thorlund JB, Ageberg E, Englund M. Knee extensor strength and body weight in adolescent men and the risk of knee osteoarthritis by middle age. Ann rheumatic dis (2017) 76(10):1657–61. doi: 10.1136/annrheumdis-2016-210888
38. Hall M, Wrigley TV, Kasza J, Dobson F, Pua YH, Metcalf BR, et al. Cross-sectional association between muscle strength and self-reported physical function in 195 hip osteoarthritis patients. Semin Arthritis rheumatism (2017) 46(4):387–94. doi: 10.1016/j.semarthrit.2016.08.004
39. Rask-Andersen M, Karlsson T, Ek WE, Johansson Å. Genome-wide association study of body fat distribution identifies adiposity loci and sex-specific genetic effects. Nat Commun (2019) 10(1):339. doi: 10.1038/s41467-018-08000-4
40. Duclos M. Osteoarthritis, obesity and type 2 diabetes: The weight of waist circumference. Ann Phys Rehabil Med (2016) 59(3):157–60. doi: 10.1016/j.rehab.2016.04.002
41. Beavers KM, Ambrosius WT, Nicklas BJ, Rejeski WJ. Independent and combined effects of physical activity and weight loss on inflammatory biomarkers in overweight and obese older adults. J Am Geriatrics Soc (2013) 61(7):1089–94. doi: 10.1111/jgs.12321
42. Livshits G, Zhai G, Hart DJ, Kato BS, Wang H, Williams FM, et al. Interleukin-6 is a significant predictor of radiographic knee osteoarthritis: The chingford study. Arthritis rheumatism (2009) 60(7):2037–45. doi: 10.1002/art.24598
43. Stannus O, Jones G, Cicuttini F, Parameswaran V, Quinn S, Burgess J, et al. Circulating levels of IL-6 and TNF-α are associated with knee radiographic osteoarthritis and knee cartilage loss in older adults. Osteoarthritis cartilage (2010) 18(11):1441–7. doi: 10.1016/j.joca.2010.08.016
44. Chang J, Liao Z, Lu M, Meng T, Han W, Ding C. Systemic and local adipose tissue in knee osteoarthritis. Osteoarthritis cartilage (2018) 26(7):864–71. doi: 10.1016/j.joca.2018.03.004
45. Karstoft K, Pedersen BK. Skeletal muscle as a gene regulatory endocrine organ. Curr Opin Clin Nutr Metab Care (2016) 19(4):270–5. doi: 10.1097/MCO.0000000000000283
46. Giudice J, Taylor JM. Muscle as a paracrine and endocrine organ. Curr Opin Pharmacol (2017) 34:49–55. doi: 10.1016/j.coph.2017.05.005
47. Sridhar MS, Jarrett CD, Xerogeanes JW, Labib SA. Obesity and symptomatic osteoarthritis of the knee. J Bone Joint Surg Br vol (2012) 94(4):433–40. doi: 10.1302/0301-620X.94B4.27648
48. Maquet PG, Pelzer GA. Evolution of the maximum stress in osteo-arthritis of the knee. J biomechanics (1977) 10(2):107–17. doi: 10.1016/0021-9290(77)90074-4
49. Reilly DT, Martens M. Experimental analysis of the quadriceps muscle force and patello-femoral joint reaction force for various activities. Acta orthopaedica Scandinavica (1972) 43(2):126–37. doi: 10.3109/17453677208991251
50. Messier SP, Gutekunst DJ, Davis C, DeVita P. Weight loss reduces knee-joint loads in overweight and obese older adults with knee osteoarthritis. Arthritis rheumatism (2005) 52(7):2026–32. doi: 10.1002/art.21139
51. Fang T, Zhou X, Jin M, Nie J, Li X. Molecular mechanisms of mechanical load-induced osteoarthritis. Int orthopaedics (2021) 45(5):1125–36. doi: 10.1007/s00264-021-04938-1
52. Boer CG, Hatzikotoulas K, Southam L, Stefánsdóttir L, Zhang Y, Coutinho de Almeida R, et al. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations. Cell (2021) 184(18):4784–818.e17. doi: 10.1016/j.cell.2021.07.038
53. Zengini E, Hatzikotoulas K, Tachmazidou I, Steinberg J, Hartwig FP, Southam L, et al. Genome-wide analyses using UK biobank data provide insights into the genetic architecture of osteoarthritis. Nat Genet (2018) 50(4):549–58. doi: 10.1038/s41588-018-0079-y
Keywords: obesity, anthropometric measurement, body compositions, knee arthritis, hip arthritis, genetic correlation, Mendelian randomization, causal association
Citation: Wang C, Zhu Y, Liu Z, Long H, Ruan Z and Zhao S (2022) Causal associations of obesity related anthropometric indicators and body compositions with knee and hip arthritis: A large-scale genetic correlation study. Front. Endocrinol. 13:1011896. doi: 10.3389/fendo.2022.1011896
Received: 04 August 2022; Accepted: 20 September 2022;
Published: 29 September 2022.
Edited by:
Ting Zheng, Hospital for Special Surgery, United StatesReviewed by:
Xiaoguang Xu, The University of Manchester, United KingdomCarolina Muniz Carvalho, Universidade Federal de São Paulo, Brazil
Copyright © 2022 Wang, Zhu, Liu, Long, Ruan and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shushan Zhao, zhaoshuiquan@126.com
†These authors have contributed equally to this work