Genetic Regulation of Pigment Epithelium-Derived Factor (PEDF): An Exome-Chip Association Analysis in Chinese Subjects With Type 2 Diabetes

Pigment epithelium-derived factor (PEDF) is a multifunctional glycoprotein that belongs to the serine protease inhibitor (serpin) superfamily (1). PEDF is widely distributed in multiple organs and tissues, such as kidney, eye, liver, lung, and adipose tissues (1,2). PEDF possesses diverse biological functions in different tissues, including antiangiogenesis, retina protection, anti-inflammation, antifibrosis, stem cell renewal, neurogenesis, and neuroprotection (1,3–5). Previous studies reported that circulating PEDF levels were significantly higher in patients with type 2 diabetes (T2D) than in subjects without diabetes (6). Serum PEDF in T2D patients was also shown to be elevated after treatment with metformin, an antihyperglycemia agent (7). On the other hand, it has been suggested that PEDF plays a protective role against diabetic microvascular damages (4,8). Dysregulation of PEDF expression was reported to be involved in the pathogenesis of microvascular complications of diabetes (4,8). Downregulation of PEDF was observed in ocular tissues of patients with diabetic retinopathy (DR) (9,10), and the role of PEDF in diabetic nephropathy (DN) was suggested by the finding of reduced PEDF expression in kidneys of diabetic animals (4). PEDF has also been shown to suppress the expression of angiogenic, fibrogenic, and proinflammatory factors (5,11,12), thereby contributing to the pathological changes in early DN.

Adipose tissue and liver have been considered the predominant sources of circulating PEDF (2). In contrast to the observations in tissues and organs, elevated circulating PEDF levels have been reported in patients with DR and DN in both type 1 diabetes and T2D (13–16). In subjects with T2D, we had also shown that increased serum PEDF levels independently predicted nephropathy progression, in particular the development of albuminuria among subjects with relatively well-preserved renal function at baseline (16). It has been suggested that the increased circulating level of PEDF might represent a compensatory systemic response to the reduced PEDF expression in the diseased tissues and organs (1). Given the multifunctional properties of PEDF, genetic polymorphisms at or near the PEDF gene have been reported to be associated with various diseases, such as age-related macular degeneration (17), coronary artery disease (CAD) (18), overall adiposity, and obesity-related insulin resistance (19). Furthermore, several studies have previously identified pathogenic mutations or common variants of the SERPINF1 (Serpin Family F Member 1) gene, which encodes the PEDF protein, as genetic determinants of circulating PEDF levels (19–25). However, to date, no genome-wide or exome-wide association studies on circulating PEDF levels have been published. Furthermore, no other study has evaluated the causal role of PEDF on the risk of diabetes complications using a Mendelian randomization approach. Therefore, we performed the current study to identify the genetic determinants of circulating PEDF levels with an exome-chip association analysis using a custom Illumina HumanExome BeadChip (Asian Exome-chip) in Chinese subjects with T2D. We then examined the relationship between circulating PEDF and diabetes complications, including DN, sight-threatening DR (STDR), and CAD, and explored the causal effect of PEDF on diabetes complications.

This study consisted of an exome-chip association study evaluating the genetic determinants of the circulating PEDF level in 5,385 Chinese subjects with T2D. The discovery stage involved 2,936 subjects who had been examined in our previous exome-chip association studies (26–28), followed by a replication study in 2,449 independent subjects with T2D who had not been genotyped with the exome chip. The cross-sectional associations of circulating PEDF levels and the PEDF-associated single nucleotide polymorphisms (SNPs) with diabetes complications were then examined. Mendelian randomization analyses were also performed to investigate the causal role of PEDF on diabetes complications in the Hong Kong West Diabetes Registry (HKWDR) cohort.

All subjects were recruited from the HKWDR, consisting of unrelated Chinese subjects with T2D regularly followed up at the medical specialist clinics of the Hong Kong West Cluster. The study protocol was approved by the institutional review board of The University of Hong Kong/Hospital Authority Hong Kong West Cluster. Written informed consent was obtained from all recruited subjects prior to any study-related procedures. Details of the study cohort; definitions of DN, STDR, and CAD; and measurement of serum PEDF levels are described in the Supplementary Data.

All subjects were genotyped using a specially designed Asian exome chip at the Centre for Genomic Sciences of The University of Hong Kong. Details of the Asian exome chip have been described previously (26–28). Sample-level quality control was carried out with regard to sex mismatch, duplication, biological relatedness, and possible sample contamination. A principal component (PC) analysis was performed using a panel of over 20,000 independent common SNPs (minor allele frequency [MAF] >0.05) for the detection of possible existence of non-Chinese samples, and outliers were removed. For SNP-level quality control, variants with MAF <0.1%, or >2% missingness, significantly deviated from Hardy-Weinberg equilibrium (HWE) with P_HWE < 1 × 10⁻⁵, or originally designed for quality control purpose were excluded from the analysis. After all quality control measures, a total of 2,936 subjects and 76,951 polymorphic SNPs (MAF ≥0.1%) remained in the exome-chip association analysis. In order to account for the between-SNP linkage disequilibrium (LD), the P value–informed LD-based clumping approach with the “–clump” command implemented in PLINK version 1.9 (29) was performed. Each index SNP had the strongest association P value within each clumped region. The other variants within the same clumped region were in LD (r² ≥0.2) with the index SNP and within ±500 kb from it.

In the replication stage, all eight index SNPs that achieved a P_discovery < 5 × 10⁻⁵ were selected for replication in another 2,449 subjects with T2D. Seven of the selected SNPs were genotyped by the MassARRAY Sequenome iPLEX Gold platform (San Diego, CA) at the Centre for Genomic Sciences of The University of Hong Kong. Genotyping of SERPINF1 rs1136287 was conducted with the TaqMan predesigned genotyping assay (Assay ID: C___1841779_20) according to the manufacturer’s instructions. SNPs that showed a low genotyping call rate of less than 90% or significantly deviated from HWE (P_HWE < 0.006 [= 0.05/8]) were excluded from further analysis. TUBGCP6 rs76062207 was excluded from further analysis due to low genotyping call rate. Consequently, seven SNPs were included in the replication analysis. The average genotyping call rate and concordance rate of these SNPs were 99.17% and 99.63%, respectively.

All statistical analyses in the discovery and replication stages were performed with PLINK version 1.9 (29). Serum PEDF level was natural-logarithmically transformed before analysis. In the discovery stage, single-variant association analysis was carried out on the standardized residuals of PEDF level generated from the multiple linear regression analysis, under the additive genetic model, with adjustment for age, sex, and the first two PCs (model 1). The first two PCs showed P values <0.05 in the Tracy-Widom test and were therefore included in the adjustment model to control for population stratification. The test statistics, as visualized in a quantile-quantile plot (Supplementary Fig. 1), appeared well calibrated. In order to assess for adiposity-independent associations, BMI was included as an additional covariate in the multiple linear regression analysis (model 2). As the previous study demonstrated that the use of metformin may affect circulating PEDF levels (7), the use of metformin was also included as a covariate in the multiple logistic regression analysis (model 3). Genome- and exome-wide significance were defined as P < 5 × 10⁻⁸ and P < 6.53 × 10⁻⁷ (= 0.05/76,951), respectively. In the replication stage, associations of the index SNPs and circulating PEDF levels were assessed by multiple linear regression analyses with adjustment for age and sex (model 1), with BMI as an additional covariate (model 2). A Bonferroni-corrected P value < 0.007 (= 0.05/7) was used as the threshold for successful replication. Meta-analysis of the association results of the discovery and replication stages was conducted using GWAMA (Genome-Wide Association Meta-Analysis) software (30). The inverse variance fixed-effect method was used to meta-analyze the summary statistics of the two stages. The heterogeneity of effect was assessed using Cochran Q test and I² index. The cross-sectional associations between the PEDF-associated SNPs with diabetes complications were examined by multiple logistic regression analyses, with adjustment for traditional risk factors including age, sex, duration of diabetes, glycated hemoglobin (HbA_1c), and presence of hypertension.

To enhance the power for detecting low-frequency variants, a gene-based analysis was conducted using three tests implemented in RVTESTS (31) to evaluate their aggregate effect in each gene. These included 1) the unweighted combined multivariate and collapsing (CMC)-Wald burden test (32), which collapsed and combined rare variants and then performed Wald test; 2) the sequence kernel association test (SKAT) (33); and 3) the variable threshold test (34). Only missense and loss-of-function (stop-gain and splicing) variants as predicted to be damaging by KGGSeq (35) were grouped into the gene sets. Two MAF thresholds (MAF <1% and MAF <5%) were used for the CMC-Wald and SKAT methods, whereas only a MAF threshold of <5% was used for the variable threshold method. Only genes with at least five copies of rare alleles were considered (n = 10,039 for MAF <1%, and n = 11,157 for MAF <5%). The gene-based significance thresholds were defined as 0.05/10,039 = 4.98 × 10⁻⁶ for MAF <1% and 0.05/11,157 = 4.40 × 10⁻⁶ for MAF <5%.

As all subjects from the discovery and replication stages were recruited from the HKWDR cohort and measurement of circulating PEDF was conducted using the same assay kit in the same laboratory, data from the two stages were combined for analyses. Binary logistic regression analyses were used to examine for the associations of circulating PEDF with DN, STDR, and CAD. Multiple logistic regression analyses with adjustment for traditional risk factors were used to investigate the independent association of circulating PEDF.

The IBM SPSS Statistics 25 and R version 3.4.3 (available at www.r.project.org) were used for the Mendelian randomization analyses. The age- and sex-standardized PEDF level was first calculated. SERPINF1 rs1136287, SMYD4 (SET And MYND Domain Containing 4) rs7224496, and SERPINF2 (Serpin Family F Member 2) rs2070863, which showed genome-wide significant associations with circulating PEDF in the current study, were selected as the instrumental variables. The inverse-variance weighted and weighted median methods implemented in the R package “MendelianRandomization” (version 0.3.0) (36), which allow for analyses with correlated variants, were used to calculate the Mendelian randomization estimates using summarized data. The MR-Egger method was then used to test for pleiotropy (36). Details of the Mendelian randomization analysis are described in the Supplementary Data.

In the discovery stage, which involved 2,936 Chinese subjects with T2D (Table 1), 76,951 polymorphic SNPs with MAF ≥0.1% were examined for their associations with circulating PEDF level in the single-variant association analysis (Fig. 1). Of these SNPs, 48.93% altered protein composition and 19.53% were Asian-specific variants with MAF 0.1%–5%. Supplementary Table 1 shows the results of the single-variant association analysis. After LD-clumping, eight index SNPs were found to be significantly associated with circulating PEDF level at P < 5 × 10⁻⁵, including three missense variants that achieved genome-wide significance (P < 5 × 10⁻⁸) after adjustment for age, sex, and two PCs. The most significant association was observed at a missense variant of SERPINF1, rs1136287 (p.Met72Thr; P_discovery = 3.93 × 10⁻³⁵; β [SE] −0.32 [0.03]). SMYD4 rs7224496 (p.Arg131Ile; P_discovery = 1.07 × 10⁻¹⁸; β [SE] 0.25 [0.03]) and SERPINF2 rs2070863 (p.Arg33Trp; P_discovery = 3.73 × 10⁻⁸; β [SE] −0.18 [0.03]) were also found to be strongly associated with circulating PEDF level at genome-wide significance (Table 2). Similar results were obtained when BMI (model 2) and the use of metformin (model 3) were further included in the multiple logistic regression analyses (Supplementary Table 1). The overall proportion of variance in circulating PEDF level explained by these three SNPs was 0.059. After conditioning on the top SNP rs1136287 at SERPINF1, the associations of SMYD4 rs7224496 (P_conditioned = 2.05 × 10⁻⁴) and SERPINF2 rs2070863 (P_conditioned = 5.64 × 10⁻⁴) with circulating PEDF levels were attenuated but remained highly suggestive of the presence of independent effects. No other variants reached exome-wide significance in the conditional analysis.

Among the seven SNPs that passed quality control in the replication stage, which involved another 2,449 subjects with T2D (Table 1), four SNPs showed consistent direction of effect as in the discovery stage (Table 2). The associations of SERPINF1 rs1136287 (P_replication = 9.32 × 10⁻³²; β [SE] −0.33 [0.03]) and SMYD4 rs7224496 (P_replication = 5.07 × 10⁻¹⁰; β [SE] 0.19 [0.03]) were successfully replicated. SERPINF2 rs2070863 showed marginal association with circulating PEDF level in the replication analysis (P_replication = 0.021; β [SE] −0.09 [0.04]) (Table 2). Similar results were obtained when further adjusted for BMI. In the meta-analysis of the two stages, these SNPs showed even stronger associations with circulating PEDF level: SERPINF1 (rs1136287; P_combined = 2.06 × 10⁻⁵⁷; β [SE] −0.33 [0.02]; I² = 0, P_{heterogeneity} = 0.810), SMYD4 (rs7224496; P_combined = 7.56 × 10⁻²⁵; β [SE] 0.22 [0.02]; I² = 0.57, P_{heterogeneity} = 0.128), and SERPINF2 (rs2070863; P_combined = 8.22 × 10⁻¹⁰; β [SE] −0.15 [0.02]; I² = 0.73, P_{heterogeneity} = 0.059) (Table 2).

To leverage the power for detecting novel PEDF-associated genes, a gene-based association test was conducted to assess the aggregate effect of low-frequency or rare variants across genes, as the single-variant association test tends to show limited power to detect the association of individual low-frequency or rare variants compared with common variants. The ITGAV gene (P = 3.82 × 10⁻⁵; β [SE] −0.46 [0.14]) exhibited the strongest effect among others in the gene-level test (SKAT <1%). A missense variant of this gene, rs61757099 (p.Asn769Ser), has also shown a suggestive association in the single-variant association test (Table 2). Supplementary Table 2 shows the genes demonstrating suggestive associations (P < 5 × 10⁻⁴) with circulating PEDF level in the gene-based association test, including SERPINF1, which showed the strongest association with circulating PEDF level in the single-variant association test. However, none of these genes reached the exome-wide significance level in the gene-based analysis.

We further investigated the associations between circulating PEDF and both microvascular and macrovascular complications, including DN, STDR, and CAD, in the HKWDR cohort. Circulating PEDF levels showed significant associations with DN and STDR (both P_unadjusted < 0.001) but not with CAD (P_unadjusted = 0.240). Elevated PEDF levels were found to be independently associated with both DN (P_adjusted < 0.001; odds ratio [OR] [95% CI] 6.66 [4.70–9.42]) (Table 3) and STDR (P_adjusted < 0.001; OR [95% CI] 2.06 [1.52–2.78]) (Table 4) even after adjustment for the traditional risk factors.

In view of the potential causal role of PEDF in DN and STDR, we also examined the associations between the PEDF-associated SNPs and these diseases. The SMYD4 rs7224496 showed a trend for association with STDR in the discovery stage cohort (P = 0.061; OR [95% CI] 0.85 [0.72–1.01]) after adjustment for traditional risk factors, but this was not replicated in the replication stage cohort. None of the other SNPs showed significant association with DN or STDR (Table 5).

We then proceeded to examine the causal role of PEDF on the risk of DN and STDR using the three PEDF-associated SNPs (SERPINF1 rs1136287, SMYD4 rs7224496, SERPINF2 rs2070863) as instrumental variables in the HKWDR cohort (n = 5,385). All three SNPs were strongly associated with the age- and sex-standardized circulating PEDF levels, and SERPINF1 rs1136287 showed the highest F-statistic and R² value among all SNPs (Supplementary Table 3). However, there was no association between these three SNPs and the risk of DN or STDR (Supplementary Table 4). All three SNPs did not show significant association with the potential confounders after correction for multiple testing (Supplementary Table 5). Mendelian randomization estimates of the associations between an SD increase in genetically determined age- and sex-standardized PEDF and risk of DN or STDR are presented in Table 6. Based on the Mendelian randomization estimates calculated using the weighted median method, there was suggestive evidence that the genetically determined age- and sex-standardized PEDF was associated with a reduced risk of STDR (P = 0.085; OR [95% CI] 0.80 [0.62–1.03]). However, there was no evidence supporting the causal role on the risk of DN. No substantial horizontal pleiotropy effects have been observed, as indicated by the insignificant P values of intercepts calculated by the MR-Egger method (Supplementary Table 6).

The current study reported the first exome-chip association analysis on circulating PEDF level in subjects with T2D. To our knowledge, no previous studies have reported an association of genetic variants with circulating PEDF level at genome-wide significance. We have demonstrated the independent associations of circulating PEDF levels with DN and STDR and undertaken a Mendelian randomization study to investigate the causal role of PEDF on the risk of DN and STDR.

The strongest association with circulating PEDF level was observed at SERPINF1 rs1136287 (p.Met72Thr) in the single-variant association analysis, reaching genome-wide significance (P_combined = 2.06 × 10⁻⁵⁷). Notably, in a previous study, this variant was shown to be associated with a lower level of plasma PEDF in healthy control subjects (n = 72; P < 0.01) but not in patients with colorectal cancer (n = 80) (25). Furthermore, carriers of SERPINF1 pathogenic mutations have shown extremely low levels of circulating PEDF compared with noncarriers (20–24). Thus, using different approaches, the current study and previous studies (19–25) have demonstrated that the SERPINF1 gene is the major genetic determinant of circulating PEDF levels. However, our observed associations of the SMYD4 and SERPINF2 variants with circulating PEDF have provided some novel insights into the genetic regulation of PEDF, although these variants only showed suggestive associations with circulating PEDF levels in the conditional analysis. Further confirmation of these newly identified associations is required. The r² between the index SNPs of the SERPINF1, SMYD4, and SERPINF2 gene regions were <0.2 in our population (Supplementary Fig. 2). SERPINF1 is flanked by SERPINF2 located ∼6.5 kb upstream on the same DNA strand and SMYD4 located ∼1.5 kb downstream on the reverse DNA strand. As discussed in the following sections, SERPINF2 and SMYD4 may also be involved in the regulation of PEDF via different signaling pathways.

SERPINF2 rs2070863 (p.Arg33Trp) was found to be significantly associated with decreased circulating PEDF levels in the current study. SERPINF2 encodes the α-2-antiplasmin (37). α-2-Antiplasmin specifically induces TGFβ1 production during the process of fibrosis (38). On the other hand, TGFβ1 has been reported to show direct positive regulation (39,40) and indirect negative regulation on the PEDF expression (41), depending on the cell environment and cell type used in the studies. It is possible that SERPINF2 rs2070863 leads to a decreased expression of α-2-antiplasmin and hence a reduction in TGFβ1 production, thereby resulting in reduced TGFβ1-induced PEDF expression.

The missense variant of SMYD4, rs7224496 (p.Arg131Ile), was associated with elevated circulating PEDF levels in the current study. This variant has been predicted to be functional by multiple prediction tools (SIFT score = 0.003; PROVEAN score = −3.02; PolypPhen2 HDIV score = 0.993; Polyphen2 HVAR score = 0.961). SMYD4 is a member of the SMYD gene family that comprises both an MYND domain and a SET domain and often carries the lysine methyltransferase activity for histone modification (42). Histone lysine methylation and DNA methylation are highly interrelated, and these two systems rely on each other mechanistically for normal chromatin function (43,44). One previous study showed that PEDF promoter methylation status negatively correlated with the PEDF mRNA expression level in locally advanced rectal carcinoma tissues (45). Therefore, SMYD4 may influence PEDF expression via its potential effect on epigenetic modification. However, SMYD4 has been proposed to function as a tumor suppressor through inhibiting platelet-derived growth factor (PDGF) receptor α (PDGFRα) expression (46). Binding of PDGF-A, -B and -C chains to PDGFRα induces its dimerization and autophosphorylation (47), and PDGFRα signaling has been shown to be critical for the maximum induction of c-Jun N-terminal kinases (JNK)-1 by PDGF-BB (48). Furthermore, reduced PDGFRα activation in response to PDGF-BB resulted in decreased downstream p38 signaling (49). Interestingly, PEDF expression has also been shown to be downregulated by PDGF-BB, with JNK1/2 and p38 signaling being critical for this repression (41). Taken together, we postulate that SMYD4 rs7224496 may lead to an increased SMYD4 expression, which results in the reduced expression of PDGFRα and hence decreased PDGF-BB activation of JNK-1 and p38, leading to a reduction in PDGF-BB–mediated PEDF repression and thus elevated PEDF expression and circulating PEDF levels (Supplementary Fig. 3). Nonetheless, the roles of both α-2-antiplasmin (SERPINF2) and SMYD4 in the regulation of PEDF expression remains to be elucidated in further functional analyses.

Despite the strongly significant associations of the SERPINF1, SMYD4, and SERPINF2 SNPs with circulating PEDF levels, we were unable to detect significant association of these SNPs with diabetes complications. Given the potential involvement of PEDF in diabetes complications, whether these PEDF-associated variants may be involved in the development of diabetes complications warrants further investigation in other independent cohorts. Similarly, further replications are also required to validate our findings on the gene-based association analysis.

We demonstrated the independent associations of circulating PEDF levels with DN and STDR, but not with CAD, in the HKWDR cohort, in line with the findings in previous publications (13–16). Our results have provided further support for the potential use of circulating PEDF level as a biomarker for diabetic microvascular complications, such as DN and STDR. Our finding was in agreement with the notion that the increased circulating level of PEDF in disease status might represent a compensatory systemic response (1), thereby providing a protective effect against further damage. In the Mendelian randomization study, our results are not supportive of a causal role of PEDF on the risk of DN. However, we found suggestive evidence for a protective role of PEDF on STDR. Noteworthy, our Mendelian randomization analysis was obviously lacking in statistical power. Our sample sizes were only able to achieve ∼25% and ∼9% power, respectively, for detecting the causal effect of PEDF on the risk of STDR and DN, assuming a R² of 0.059 at a significance level of 0.05. Therefore, a definite conclusion on the role of PEDF on the risk of DN and STDR could not be drawn from the current study. Identification of more PEDF-associated variants to provide stronger instrumental variables will serve to improve the study power. Mendelian randomization studies with sufficient power in other populations to investigate further the causal role of PEDF on the risk of diabetes complications are warranted.

The current study has some limitations. First, this study lacks an external validation. Further replication study in independent cohorts in other populations would serve to validate our findings. Second, despite the use of the Asian exome chip, which had allowed us to examine more functional coding variants with a lower MAF compared with the conventional genotyping arrays, our study, with its relatively small sample size, has limited power to detect the association of rare or low-frequency variants. In view of this, we conducted the gene-based association analysis to assess the combined effect of these rare or low-frequency alleles. A whole-genome imputation using publicly available reference panels, such as the 1000 Genome Project, would increase the total evaluable SNPs in studies using array chips. However, the exome chip used in this study has been designed with enrichment of rare exonic variants identified primarily from the European population. This feature implies that a large proportion of variants are likely to be unique to the exome chip design and are therefore not present in the publicly available reference panels of non-European populations, such as Asians, and may in turn lead to reduced imputation accuracy. Future meta-analysis of genome-wide or exome-wide association studies would be useful to identify additional genetic variants influencing the circulating PEDF level. Moreover, the current study lacks functional analyses to validate the proposed mechanisms whereby the two newly identified genetic variants may regulate PEDF expression, which remain to be confirmed in future functional studies. Finally, we have limited power to investigate the cause–effect relationship between PEDF and the risk of diabetes complications. Further Mendelian randomization studies in larger cohorts may provide more insights into the causal role of PEDF on the risk of diabetes complications.

In summary, we have conducted an exome-chip association study and identified, for the first time, three missense variants of SERPINF1, SMYD4, and SERPINF2 showing genome-wide significant associations with circulating PEDF level in subjects with T2D. Our findings have shed light on the genetic regulation of circulating PEDF level and provided further support for its potential application as a biomarker for DN and STDR. The causal role of PEDF on diabetes complications deserves further investigations.

Acknowledgments. The authors thank Rachel Wong (Department of Medicine, The University of Hong Kong, Hong Kong, China) for her technical assistance in the measurement of serum PEDF level and all the study subjects and clinical and research staffs of the HKWDR for their contribution to this research study.

Funding. This work was supported by the Hong Kong Research Grants Council Theme-based Research Scheme (T12-705/11).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. C.Y.Y.C. and C.-H.L. wrote the first draft of the manuscript. C.Y.Y.C. analyzed and interpreted the data. C.Y.Y.C., C.-H.L., A.X., K.-W.A., C.H.Y.F., K.K.K.N., K.H.M.K., W.-S.C., Y.-C.W., M.M.A.Y., J.H., K.C.B.T., and T.-H.L. were involved in the sample collection and selection and phenotype data preparation. C.S.T. and P.-C.S. provided useful comments to data analysis. H.-F.T., P-.C.S., and K.S.L.L. conceived the study, undertook project leadership, and contributed equally to the supervision of this work. All authors contributed to the drafting and critical revision of the manuscript and approved the final version of the manuscript. H.-F.T., P-.C.S., and K.S.L.L. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.