Susceptible gene of stasis-stagnation constitution from genome-wide association study related to cardiovascular disturbance and possible regulated traditional Chinese medicine

Background This study identified susceptible loci related to the Yu-Zhi (YZ) constitution, which indicates stasis-stagnation, found in a genome-wide association study (GWAS) in patients with type 2 diabetes and possible regulated traditional Chinese medicine (TCM) using docking and molecular dynamics (MD) simulation. Methods Non-aboriginal Taiwanese with type 2 diabetes were recruited. Components of the YZ constitution were assessed by a self-reported questionnaire. Genome-wide SNP genotypes were obtained using the Illumina HumanHap550 platform. The world’s largest TCM database (http://tcm.cmu.edu.tw/) was employed to investigate potential compounds for PON2 interactions. Results The study involved 1,021 unrelated individuals with type 2 diabetes. Genotyping data were obtained from 947 of the 1,021 participants. The GWAS identified 22 susceptible single nucleotide polymorphisms on 13 regions of 11 chromosomes for the YZ constitution. Genotypic distribution showed that PON2 on chromosome 7 was most significantly associated with the risk of the YZ constitution. Docking and MD simulation indicated 13-hydroxy-(9E_11E)-octadecadienoic acid was the most stable TCM ligand. Conclusions Risk loci occurred in PON2, which has antioxidant properties that might protect against atherosclerosis and hyperglycemia, showing it is a susceptible gene for the YZ constitution and possible regulation by 13-hydroxy-(9E_11E)-octadecadienoic acid. Electronic supplementary material The online version of this article (doi:10.1186/s12906-015-0761-x) contains supplementary material, which is available to authorized users.


Background
Differences exist between traditional Chinese medicine (TCM) and conventional western medicine. These differences include not only the treatment approach (such as herbal medicine and acupuncture in TCM) but also the underlying theories. A principle component of TCM theory is the concept of constitution, which provides a method for classifying patients according to type.
Constitution demonstrates individual differences in structure and function, temperament, and environmental adaptability. Patients with different constitutions have different susceptibilities, development, and prognoses for certain diseases. According to Huang Di Nei Jing, a textbook of TCM internal medicine written approximately 2,000 years ago, a certain constitution is partially developed from congenital factors [1]. This view is similar to "personalized medicine," which highlights the genetic background for disease susceptibility. Genetic studies to investigate the congenital factors of constitution are increasing in the post-genome era. Chen et al. reported that allele frequencies of human leukocyte antigens including DPB1*0501 in the Yin-deficiency group (Frequency 51.6 % vs 35.6 %, relative risk 1.9), DRB1*09012 in the Phlegmwetness group (Frequency 23.4 % vs 12.8 %, relative risk 2.1), and DQB1*03032 in the Qi-deficiency (Frequency 22.2 % vs 8.1 %, relative risk 3.2) and Phlegm-wetness groups (Frequency 19.9 % vs 8.1 %, relative risk 2.8) differ significantly from those in the normal constitution [2]. Wang et al. conducted an expression array and identified 785 upregulated genes and 954 downregulated genes in the Yang-deficiency constitution, compared with those in normal individuals. The most significant enriched Gene Ontology Cluster of upregulated genes is "response to stress" which contained interleukin factors and their receptors. The most significant enriched Gene Ontology Cluster of downregulated genes is "nucleobase, nucleoside, nucleotide and nucleic acid metabolism" which contained thyroid hormone receptor signal pathway [3]. A study of polymorphisms further identified the biased distribution of single nucleotide polymorphisms (SNPs) in PPARD (peroxisome proliferator-activated receptors delta) rs2267669 and rs2076167 and APM1 (adipose most abundant gene transcript 1) rs7627128 and rs1063539 in the Yang-deficiency constitution; PPARD rs2076167 and APM1 rs266729 and rs7627128 in Phlegm-wetness constitution; and in PPARG (peroxisome proliferator-activated receptors gamma) Pro12Ala in the Yin-deficiency constitution [4]. Gene expression is influenced by environmental factors in the posttranscriptional process and candidate gene studies are limited to a certain viewing region. Cardiovascular disease is a major complication in patients with diabetes mellitus (DM), especially in those with type 2 diabetes, resulting in both comorbidity and mortality [5]. One meta-analysis reported that DM tends to double the risk of cardiovascular disease [6]. Screening for high cardiovascular risk in patients and providing more effective protection for these patients are important in clinical practice. The Yu-Zhi (YZ) constitution in TCM indicates stasis and stagnation, which expressed dull, lusterless skin color; dry, cracked, scaly or tough skin; dull purple lips or tongue; localized pain or numbness; knotted, intermittent, or uneven pulse. It is one of the body constitutions that tend to express blood stasis syndrome (BSS), a morbid state caused by blood circulation disturbance, included extravasated blood, blood circulating sluggishly, or blood congested in viscera, that may turn into pathogenic factors. BSS is usually considered a link to cardiovascular complications. A study of hospitalized patients with coronary artery disease (CAD) noted that BSS was the most common TCM pattern in over three-quarters of the patients [7]. The risk of BSS increases with the carotid intima-media thickness in patients with dyslipidemia [8]. According to TCM theory, BSS constitutes the main mechanism of cardiovascular diseases, including diabetic cardiovascular complications. The pathogenesis of BSS includes microcirculation disturbance, abnormal hemorheological factors, and hemodynamic changes. Some small samples molecular studies also detected differences in cell surface antigens and gene expression between those with BSS and healthy controls [9,10]. As mentioned earlier, congenital factors are considered to be a principal component of constitutional formation; possible genetic variations underlying the YZ constitution are of interest. In the present study, a genome-wide association study was conducted to identify susceptible loci related to the YZ constitution in patients with type 2 diabetes. Genes related to the susceptible loci are considered susceptible genes of the YZ constitution. Computer-aided drug design (CADD) has been widely used in studies investigating new treatments [11][12][13], and could help accelerate the development of leading drugs [14,15]. CADD could be employed to approaches in the design of drugs for anti-inflammation [16], anti-virus [17,18], pain regulation [19], weight loss [20,21], stroke therapy [22][23][24], and cancer therapy [25][26][27][28]. Hence, we employed a TCM database (http://tcm.cmu.edu.tw) [29] and natural compounds to conduct virtual screening for proteins of susceptible genes by molecular docking to find potential TCM or natural compounds. Then we performed molecular dynamics (MD) simulation to study the protein-ligand interactions and stabilized conformations for the top candidates.

Study participants
Our study population comprised adult patients with type 2 diabetes in Taiwan. Patients willing to participate in a genetic study in nonaboriginal Taiwanese were recruited from the outpatient clinic of China Medical University Hospital (Taiwan) between September 2006 and June 2007. Type 2 diabetes was diagnosed according to the criteria of the 1997 American Diabetes Association [30]. Patients with type 1 diabetes, gestational diabetes, or maturity-onset diabetes of the young were all excluded. This research was approved by the China Medical University Hospital Institutional Review Board, and all participants gave informed consent.

Data collection
Patient information (age, sex, age at diagnosis of diabetes, smoking history) was collected by a questionnaire. Systolic and diastolic blood pressures were obtained by averaging two measurements with a resting interval of at least 5 min. Patients with hypertension were defined as having a systolic pressure of more than 130 mmHg, and a diastolic pressure of more than 80 mmHg, or those who were receiving antihypertensive agents. The body height and weight of subjects (wearing light clothing and no shoes) were measured by experienced research staff. The body mass index was calculated by dividing the body weight (kg) by the square of the body height (m). Blood samples were drawn between 8:00 and 10:00 a.m. after the patients had fasted overnight, and separated serum was stored at −70°C until assayed. Fasting plasma glucose was detected by the hexokinase method. Serum total cholesterol, triglycerides, high-(HDL) and low-density lipoprotein (LDL) cholesterol, creatinine, and uric acid levels were measured by standard laboratory methods. Highsensitivity C-reactive protein was measured by immunoturbidimetry (Integra 700; Roche, Mannheim, Germany). Hemoglobin A1c was gauged by the high-performance liquid chromatography method (HLC-723G7; TOSOH Bioscience, Tokyo, Japan). The albumin-to-creatinine ratio (ACR) was obtained from a morning spot urine test and data were categorized as normoalbuminuria (ACR ≦ 30 mg/g), microalbuminuria (30 mg/g < ACR < 300 mg/g) or macroalbuminuria (ACR ≧ 300 mg/g). Biochemical analyses were performed at the Taipei Institution of Pathology.

Yu-Zhi constitution questionnaire
The YZ body constitution was assessed by a questionnaire, which was developed using a psychometrically sound method as reported previously [31], and comprises eight self-reported symptomatic items(Supp.) [32]. Each item was assessed on a 5-point Likert scale (never, occasionally, sometimes, often, and always). The YZ score was obtained by a summation of the scores of the eight items, and a higher score indicated stronger intensity of the YZ constitution.

Genotyping
Genomic DNA was extracted from peripheral blood mononuclear cells for a genome-wide association study. The procedures for genomic DNA extraction, whole genome genotyping, genotype calling and quality control have been described previously [33]. PUREGENE DNA isolation kit (Gentra Systems, Minneapolis, MN) was use to extract genomic DNA from peripheral blood mononuclear cells, and Illumina HumanHap550-Duo Bead-Chips was used to perform the whole genome genotyping in deCODE genetics (Reykjavı'k, Iceland). The standard procedure implemented in BeadStudio, with default parameters suggested by the platform manufacturer was used to perform Genotype calling. Genotyping validation was performed using the Sequenom iPLEX assay (Sequenom MassARRAY system; Sequenom, San Diego, CA, USA). By examining several summary statistics, quality control of the genotype data was performed. First, by calculating the ratio of loci having heterozygous calls on the X chromosome, sex of the patients was doublechecked. Second, total successful call rate and the minor allele frequency were also calculated for each SNP. The exclusion situations for SNPs were as follows: no polymorphism, a total call rate of < 95 %, or a minor allele frequency of < 5 % and a total call rate < 99 %. A total of 560,184 SNPs were genotyped, 38,700 SNPs were excluded due to quality control criteria, 12,723 SNPs were excluded due to Hardy-Weinberg equilibrium principle (P <0.0001) and 508,761 SNPs were used in final analysis.

Statistical analysis
Data from continuous variables were expressed as mean ± standard deviation and categorical variables were expressed as percentages. Participants were divided into high and low YZ score groups according to the median score (YZ score = 10). Differences between the high and low YZ score groups were compared using Student's ttest for continuous variables and the χ 2 test for categorical variables. Association analysis was carried out to compare allele frequency and genotype distribution between the high and low YZ score groups using 6 single point methods for each SNP: genotype, allele, trend (Cochran-Armitage test), additive, dominant, and recessive models using PLINK (PLINK 1.07, http://pngu.mgh. harvard.edu/~purcell/plink/contact.shtml#cite) and SAS (SAS Institute Inc., 100 SAS Campus Drive, Cary, NC 27513-2414, USA).
The associated SNPs were selected from those at least showing p-values < 10 −5 under the most significant test statistic obtained from any of the 6 statistical models. We then used a multivariate logistic regression model to determine the genotype odds ratios (ORs) and 95 % confidence intervals (CIs) of associated SNPs in the best model. For ORs, p-values < 0.05 were considered statistically significant.

Structure preparation and docking study
Because the PON2 structure is not available in the PDB database, the sequence of PON2 (UniProt ID: Q15165) was obtained from the UniProt database for 3D structure modeling. The sequence of PON2 was submitted to the I-TASSER server [34][35][36] (iterative threading assembly refinement algorithm) to generate a 3D structure of PON2. For 3D structure validation, we also employed a Ramachandran plot [37], profile-3D (Discovery Studio Client v2.5; Accelrys, San Diego, CA, USA.) and PONDR-FIT [38] (DisProt, http://www.disprot.org/) to validate the PON2 modelling structure. The residue His114 of PON2 was regarded as the binding site for screening TCM compounds based on protein-ligand interaction [39]. We used the LigandFit module of DS 2.5 to calculate docking poses of TCM compounds in the PON2 protein structure. Each  [40]. The minimization step in each docking pose was performed with 1000 steps of the Steepest Descent and Conjugate Gradient. The generated conformations of ligands were docked into the defined binding site of the PON2 modeling structure. All of the TCM compounds with docking poses had various scoring functions including the piecewise linear potential (−PLP1, −PLP2), potential of mean force (−PMF), and Dock Scores.

Molecular dynamics simulation
The MD simulation was carried out by the GROMACS 4.5.5 package [41] to simulate the dynamic condition of the protein-ligand complex from the docking results of PON2. The edge of the box between the protein complexes was set as 1.2 nm. We choice the charmm27 force field for the simulation system [42]. The proteinligand complex containing water molecules was placed by TIP3P model cubic cells. Non-bonded interactions included Coulomb terms and van der Waals (VDW). The particle mesh Ewald method [43,44] was used to define Coulomb interactions as electrostatic in this experiment, and the cut-off distance of VDW residues was set at 1.4 nm. The linear constraint solver algorithm was used to fix all bond lengths among all atoms of the protein-ligand complexes. Topology files and parameters of TCM compounds in PON2 complexes were generated from the SwissParam web server [45] for the GROMACS simulation. For ion setting of solvent concentrations, the Na + and Cl − ions were randomly replaced the water molecules in the solvent system with a concentration of 0.145 M. The minimization of energy was used to stabilize the protein-ligand complex by the steepest descent method with 5000 steps, followed by equilibration performed under position restraints with 1 ns for balance of the water molecules between PON2 complexes. The condition of equilibration was under constant temperature dynamics (NVT type) at a temperature of 310 K. In the final step, a production run was performed for 20,000 ps with constant pressure and temperature dynamics (NPT type). The temperature of the production run was set at 310 K. All frames of MD conformation were sampled every 20 ps for trajectory analysis under GROMACS4.5.5.

Results
This study enrolled 1,021 patients with type 2 diabetes who were 20 years old or older. Participants were divided into a high YZ score group and low YZ score group according to the median (numbers of high YZ score: low YZ score = 583: 438). The patient information and clinical characteristics of the high and low YZ score groups are summarized in Table 1. Compared with the low YZ score group, the high YZ score group was significantly older (61.3 ± 11.4 vs 59.7 ± 10.5 years; p = 0.020), included more women (53.5 vs 47.3 %; p = 0.048), had a longer duration of diabetes (11.9 ± 7.4 vs 10.8 ± 6.8 years; p = 0.017), and displayed poorer control of serum glucose (hemoglobin A1c 8.0 ± 1.5 vs 7.8 ± 1.4 %; p = 0.015).
Genotyping data were obtained from 947 (numbers of high YZ score: low YZ score = 539: 408) of the 1,021 participants using Illumina HumanHap550duov3 chips. Quantile-quantile plots for each model were shown that the distribution of observed p-values deviated from expected p-values in Fig. 1. Manhattan plots of p-values across all chromosomes for each model were shown in Fig. 2. SNPs in autosomal chromosomes with a p-value < 9.8 × 10 −8 were not detected in all the 6 statistical  models. Table 2 Table 3 shows the results of multiple logistic regression analysis of the genotypic distribution of susceptible SNPs in patients with high YZ scores among the high and low YZ score groups. The results showed that 20 SNPs in 11 regions of 10 chromosomes were significantly associated with high YZ scores in the best model, after controlling for age, sex, diabetes duration, and hemoglobin A1c. The risk genotypes were defined by homozygous risk alleles (higher allele frequency in the high YZ score group than Fig. 3 Ramachandran plot of the PON2 modelling structure. The modelling structure was built from the ITASSER server. There are 86.6 % and 5.1 % of residues of PON2 in the *favoured and disfavoured regions, respectively The PON2 protein was built from the I-TASSER server and we validated the simulated PON2 structural residues by Ramachandran plot (Fig. 3), with 86.6 % of residues of PON2 located in the favoured region, and only 5.1 % of residues in the disfavoured regions. We also used 3Dprofiling to observe the reliability of each residue. Most PON2 residues had validation scores with positive score values (Fig. 4). The score for the active residue, His144, displayed that it was reliable from the modelling structure, indicating that this key residue was not affected by the docking screen process. For protein disorder analysis (Fig. 5), most of the sequence (including residue His114) revealed that the disorder disposition was below 0.5. The prediction of disorder analysis illustrated that the PON2 protein is a folded structure, and may not affect TCM compounds during the docking process [46,47].
The docking analysis was based on -PLP1, −PLP2, −PMF, and Dock Scores to select the docking poses of TCM compounds from database screening. From the scoring analysis, we analysed the top ten TCM ligands with high -PMF scores as candidates (Table 4). Furthermore, we found the docking poses of the top three candidates, divaricatacid, 13-hydroxy-(9E_11E)-octadecadienoic acid, and 9-hydroxy-(10E)-octadecenoic acid, were very close to the key residue His114 (Fig. 6). The docking poses showing that the three candidates can interact with His114, and have the potential to activate PON2 for antioxidation. In a further study, we performed MD simulation to observe the stability of TCM candidates in the PON2 structures under dynamic condition.
For the trajectory analysis of the MD simulation, the root mean square deviation (RMSD) and gyration of    protein atoms were used to observe the stability of the protein structure for PON2 and protein-ligand complexes. The value of the protein RMSD was around 0.2 and 0.3 nm from 2 ns to 20 ns (Fig. 7a). The apo form of the PON2 and PON2 complexes with TCM candidates revealed stable fluctuation during all of the MD simulation time, and 20 ns of simulation time facilitated all simulation systems into constant conditions. The gyration of the PON2 structure was state in 2.05 from 8 ns to 20 ns (Fig. 7b). The apo form revealed a decreasing value of gyration, which was more compact than the PON2 complexes with TCM candidates, illustrating that TCM compounds not leave aware from the PON2 structure during all MD simulation times. For the ligand RMSD, we can found that 9-hydroxy-(10E)-octadecenoic acid displayed fluctuating RMSD values from 6 ns to 20 ns. Divaricatacid and 13-hydroxy-(9E_11E)-octadecadienoic acid revealed stable ligand RMSDs during the overall simulation time (Fig. 7c). Hence, 9-hydroxy-(10E)-octadecenoic acid may not be suitable to interact with the PON2 structure. The total energy of the PON2 complexes ( Fig. 8) with divaricatacid, 13-hydroxy-(9E_11E)-octadecadienoic acid, and 9-hydroxy-(10E)-octadecenoic acid was −8.74 × 10 5 during the initial simulation time (from 0 ns to 2 ns), but in the final step of the simulation time, the total energy tended to be stable at −8.78 × 10 5 which was similar to the apo form of PON2. The results for total energy show that all systems of PON2 are stable after 20 ns of MD simulation time.
We further calculated the root mean square fluctuation (RMSF) values for each residue of the PON2 protein structure (Fig. 9). Interestingly, we found the residues revealed high fluctuations from 50 to 100 on the apo form of PON2 (Fig. 9a). The high RMSF values indicate the residue was flexible at all simulation times. The PON2 structure with docked TCM compounds was more stable than the apo form, which suggests that TCM compounds can stabilize the protein structure in the protein-ligand complex type. Because of the flexible residues on the apo form of PON2 near the key residue His114, reducing the variation of these residues denotes that TCM compounds could tightly interact with the PON2 structure. We also calculated the solvent accessible hydrophobic and hydrophilic surface areas for the three TCM compounds. The results show that divaricatacid is suitable for hydrophobic solvents, because the hydrophilic areas of 13-hydroxy-(9E_11E)-octadecadienoic acid and 9-hydroxy-(10E)-octadecenoic acid were wider than that of divaricatacid (Fig. 10).
To analyze the migration of each TCM compound during the MD simulation time, we computed the mean square displacement (MSD) value to measure the variation of each ligand in the PON2 structure. 9-hydroxy-(10E)-octadecenoic acid displayed a significantly higher MSD value during the 20 ns simulation time than the other two candidates (Fig. 11a). In addition, we further calculated the distance between PON2 and the three candidates. The distance between 13-hydroxy-(9E_11E)octadecadienoic acid and the PON2 structure did not change much during the simulation time of 20 ns (Fig. 11b), but divaricatacid and 9-hydroxy-(10E)-octadecenoic acid gradually moved away from the PON2 structure.
We employed CAVER 3.0 software [48] to predict the ligand tunnels in the PON2 structure for the three TCM candidates (Fig. 12). The predicted tunnels are represented by red, blue, green, yellow, cyan and orange solid phases. The apo form of the PON2 structure revealed a broad space of tunnels (Fig. 12a), because there was no docked ligand in the PON2 binding site. A comparison of the MSD values showed that 13hydroxy-(9E_11E)-octadecadienoic acid was the most stable for all MD simulation times, and hence, the predicted tunnel reveals a narrow space in Fig. 12c, which is represented by the green solid phase. For the other two TCM candidates, divaricatacid and 9-hydroxy-(10E)octadecenoic acid, the predicted tunnels display spacious areas in Fig. 12b and 12c. In the 20 ns snapshot (Fig. 13), 13-hydroxy-(9E_11E)-octadecadienoic acid is still close to the key residue His114, which confirms the ligand RMSD, migration analysis and the measurement of the proteinligand distance (Fig. 13b). The final snapshot showed large distances between His114 and both divaricatacid and 9hydroxy-(10E)-octadecenoic acid ( Fig.13a and 13c). This illustrates that 13-hydroxy-(9E_11E)-octadecadienoic acid is the best potential TCM compound to interact with the PON2 structure.

Discussion
The results of this genome-wide association study identified 22 YZ constitutionally susceptible SNPs, representing 13 regions of 11 chromosomes. Genotypic distribution showed that high YZ scores were significantly associated with PON2 on chromosome 7, PIEZO2 on chromosome 18, ZNF665 on chromosome 19, FREM2 on chromosome 13, and unknown genes on chromosomes 1p, 2q, and 16p. PCDH10 on chromosome 4q, and unknown genes on chromosome 5q and chromosome 14q were significantly associated with lower risk of the YZ constitution after controlling for age, sex, diabetes duration, and hemoglobin A1c.
Without doubt, the YZ constitution is a consequence of complicated polygenic influences. Genome-wide association studies can provide an overview of whole genomes, and this is an appropriate method for examining the genetic factors of the YZ constitution. This technique had been adapted to explore the genetic base of Korean Sasang constitutional medicine [49].
Common manifestations of YZ include dull, lusterless skin color; dry, cracked, scaly or tough skin; dull purple lips or tongue; and localized pain or numbness. A patient with a YZ constitution tends to express BSS, which according to TCM theory, indicates a morbid state of blood The rs7493, rs2299263, and rs17166875 polymorphisms, located in PON2 on chromosome 7q, belong to one of the paraoxonase (PON) gene families, which encode enzymes participating in the hydrolysis of organophosphates. The PON gene cluster contains 3 adjacent gene members, PON1, PON2, and PON3. All 3 PON genes share high sequential homology and a similar β propeller protein structure [50] and are thought to have antiatherosclerotic properties. Thus the PON gene cluster has been considered a target in the treatment of atherosclerosis [51,52]. PON2 has been shown to prevent LDL oxidation, to reverse the oxidation of mildly oxidized LDL, and to inhibit oxidized LDL-induced monocyte chemotaxis [53]. It also increases cholesterol efflux [54] and decreases the size of atherosclerotic lesions [55].
PON2 is a ubiquitously expressed intracellular protein that is expressed in a wide range of tissues [56,53]. PON2 exhibits antioxidant functions at the cellular level, in addition to a host of intracellular antioxidative enzymes that act against oxidative stress. PON2 is localized in the inner mitochondrial membrane, associated with respiratory complex III, and binds with high affinity to coenzyme Q10. Decreased activity of mitochondrial electron transport chain (ETC.) complexes is implicated in the development of many inflammatory diseases, including atherosclerosis. PON2 protects ETC. complexes against oxidative stress by lowering reactive oxygen species. The intracellular antioxidative effect plays a role in antiatherosclerosis by avoiding endothelial dysfunction caused by mitochondria dysfunction [57,58]. A common polymorphism rs7493, also known as Ser311Cys, a missense SNP in PON2, has also been associated with the risk of CAD [59]. In addition, PON2 plays a role in hepatic insulin signalling. PON2-deficient mice display elevated hepatic oxidative stress, coupled with an exacerbated inflammatory response, because of PON2-deficient macrophages. PON2 deficiency is associated with inhibitory insulin-mediated phosphorylation of hepatic insulin receptor substrate-1. PON2 may enhance the influence of the macrophage-mediated inflammatory response in hepatic insulin sensitivity [60]. The PON2 G148 variant has been associated with elevated fasting plasma glucose in patients with type 2 diabetes [61]. The role of PON2 provides the genetic basis underlying the YZ constitution. Patients with a strongly YZ constitution may have PON2 polymorphism with a low protein function which tends to  decrease its antioxidative efficacy, resulting in cardiovascular disturbance and hyperglycemia. Thus, PON2 may be a candidate gene for the YZ constitution. Treatment using herbal medicines or natural compounds that could potentially regulate PON2 might be useful in protecting type 2 diabetes patients with a YZ constitution from cardiovascular complications.
From the docking results of TCM database screening, we chose three potential TCM candidates based on -PMF scores, divaricatacid, 13-hydroxy-(9E_11E)-octadecadienoic acid and 9-hydroxy-(10E)-octadecenoic acid. We further simulated the interaction between PON2 and TCM compounds under dynamic conditions for 20 ns. 13-hydroxy-(9E_11E)-octadecadienoic acid was more stable than the other two candidates for binding with the PON2 structure, which was still connected with active residue His114 after an MD simulation time of 20 ns. According to this result, 13-hydroxy-(9E_11E)-octadecadienoic acid should be a ligand with the ability to regulate PON2. 13-hydroxy-(9E_11E)-octadecadienoic acid is isolated from the seed of Coix lacryma-jobi L. Coix oil had been reported the efficacy to decrease adipose tissue and LDL concentrations and increase the total antioxidant capacity in hyperlipidemic rats [62]. 13-hydroxy-(9E_11E)octadecadienoic acid may play a role in antiatherosclerosis by avoiding endothelial dysfunction by regulating the antioxidant effect of PON2.
The polymorphisms of rs1133146, rs12971799, rs4801958, rs12460170, and rs4803055 are located in the ZNF665 of chromosome 19q, belonging to the Kruppel zinc finger family. Zinc fingers are the most abundant DNA-binding motifs in humans. The zinc finger protein families are mainly involved in recognizing DNA sequences, but are also able to bind RNA, DNA-RNA hybrids, and even proteins [63]. They work as transcription factors to interact with the control region and achieve gene expression. Kruppel type zinc finger genes are widely present in the human genome, and are usually involved in cell growth and differentiation. To date, specific details of the function of ZNF665 have not been documented. We speculated that the polymorphisms of ZNF665 might lead to poor gene expression because of poor DNA binding ability, which might disturb cell growth and differentiation; in turn, this disturbance might impede epithelial repair and lead to the dry, cracked, scaly, or tough skin that is characteristic of patients with the YZ constitution. In addition, poor cell growth and differentiation in endothelial progenitor cells might disturb epithelial repair and lead to endothelial dysfunction, which is thought to be a key event in the development of atherosclerosis [64,65].
The rs12865228 and rs4526895 polymorphisms are located in FREM2 on chromosome 13q. This gene encodes a membrane protein that belongs to the FRAS1 family. This extracellular matrix protein forms a ternary complex localized on the basement membrane, and plays a role in epidermal-dermal interactions during morphogenetic Fig. 13 The final snapshots of PON2 with three TCM compounds: (a) divaricatacid, (b)13-hydroxy-(9E_11E)-octadecadienoic acid, and (c) 9-hydroxy-(10E)-octadecenoic acid from the results of MD simulation processes [66]. The protein is thought to be necessary in maintaining the integrity of skin epithelium, vascular stability [67], and the differentiated state of renal epithelia [68]. The polymorphism of FREM2 in patients with high YZ scores might be implicated in skin changes and microcirculation disturbances, leading to dry, cracked, scaly, tough, and bruised skin, and a dull and lusterless face.
Our study found that SNPs located in PIEZO2 on chromosome 18p were also associated with high YZ scores among patients with type 2 diabetes. PIZEO2 is a large transmembrane protein with 24 to 36 predicted transmembrane domains, and is a component of the mechanosensitive channel. This channel is required for the rapid adaptation of mechanically activated currents in dorsal root ganglia [69]. Mechanical stimuli drive many physiological processes, including touch and pain sensation, hearing, and blood pressure regulation [70]. Dysfunction of PIZEO2 might be the source of the abnormal sensations reported by patients with the YZ constitution, including numbness, tightness, tingling pain, and a dull sensation.
PCDH10 on chromosome 4q was associated with a lower risk of the YZ constitution. PCDH10 belongs to the protocadherin gene family, a subfamily of the cadherin superfamily. PCDH10 is a putative tumor suppressor gene [71], and is also known to guide the development of axons [72]. Furthermore, several SNPs located in the intergenic area on chromosomes 1p, 2q, 5q, 14q, and 16p require further investigation to clarify their relationship with the YZ constitution.

Conclusions
The findings of this study contribute to an understanding of the genetic susceptibility of patients with type 2 diabetes to the YZ constitution. Risk loci occurred in PON2 that encoded intracellular proteins with antioxidant properties, which normally protect against atherosclerosis and hyperglycemia. Disturbance of this genetic function might constitute one of the mechanisms of cardiovascular disturbance induced by the YZ constitution. Docking and molecular dynamic simulation showed that 13-hydroxy-(9E_11E)-octadecadienoic acid is a stable ligand of PON2 that may have the ability to regulate the antioxidant effects of PON2. Other related genes included ZNF665, FREM2, PIZEO2, PCDH10 and several SNPs located in genes of unknown function.