Heritability of Cold and Heat Patterns: A Twin Study

Abstract

In traditional East Asian medicine, cold–heat patterns have been widely used in the diagnosis and treatment of patients suffering from various diseases. The present study aimed to estimate the heritability of cold–heat patterns. Trained interviewers administered a cold–heat pattern questionnaire to 1,753 twins (mean age = 19.1 ± 3.1 years) recruited throughout South Korea. Correlations for the cold pattern (CP) were 0.42 (95% CI [0.28, 0.54]) for monozygotic (MZ) males, 0.16 (95% CI [-0.08, 0.39]) for dizygotic (DZ) males, 0.40 (95% CI [0.30, 0.49]) for MZ females, 0.30 (95% CI [0.12, 0.45]) for DZ females, and 0.07 (95% CI [-0.11, 0.25]) for opposite-sex DZ twins. The corresponding twin correlations for the heat pattern (HP) were 0.38 (95% CI [0.24, 0.51]), -0.22 (95% CI [-0.43, 0.02]), 0.34 (95% CI [0.24, 0.43]), 0.21 (95% CI [0.03, 0.37]), and 0.08 (95% CI [-0.10, 0.26]), respectively. These patterns of twin correlations suggested significant genetic effects on the HP and the CP. Model-fitting analysis revealed that heritability estimates in both sexes were 40% (95% CI [38, 42]) for the CP and 33% (95% CI [25, 42]) for the HP, with the remaining variances attributable to unique environmental variances. These estimates did not vary significantly with age during adolescence and young adulthood.

In traditional East Asian medicine, the cold–heat pattern differentiation system has long been used as a part of clinical diagnosis (Jian, 2005). The cold–heat pattern differentiation system is useful in clinical practice because it improves diagnostic accuracy of the diseases by giving clinicians additional information about the nature and location of physiological imbalances of the patient (WHO Regional Office for the Western Pacific, 2007). The cold pattern (CP) is characterized by a preference for warm temperature, intolerance of cold, hypothermia, paleness, diarrhea, peripheral chills, and spasms, whereas the heat pattern (HP) is characterized by aversion to hot temperature, diaphoresis, rapid pulse, flushed face, vexation, constipation, thirst, and deep-colored urine (Jian, 2005). Incorporation of the cold–heat differentiation system in clinical diagnosis has been shown to improve treatment outcomes as well. Lu et al. (2009) divided 194 patients with rheumatoid arthritis into CP and the HP and found that as compared to the patients with HP, those diagnosed with the CP responded better to a 24-week biomedical therapy consisting of diclofenac, methotrexate, and sulfasalazine. In response to a nationwide survey, Korean Medicine doctors reported that the cold–heat differentiation system enhanced prognosis for 10 diseases, including menopausal disorders, chronic rhinitis, dyspepsia, hwa-byung, diarrhea, dysmenorrhea, headache, inflammation in the digestive tract, coldness in hands and feet, and atopic dermatitis (Bae et al., 2017).

Recently, there is growing research interest in the biological basis of cold–heat patterns. For instance, van Wietmarschen et al. (2009) classified 33 patients with rheumatoid arthritis into CP and HP and found significantly different gene expression and metabolite profiles for regulation of apoptosis between the two subtypes. Given that patients with cold–heat patterns were shown to have abnormal functions in the neuro-endocrine-immune (NEI) system (Hsu et al., 2003), Li et al. (2007) explored whether and which NEI-related genes and chemical messengers (CM: hormones, cytokines, and neurotransmitters) are related to the CP and the HP network by using a co-occurrence literature mining approach (Jenssen et al., 2001; Li et al., 2006). This approach assumes that when two biological entities are co-cited in the same text unit, there should be a potential biological relationship between them. Li et al. (2007) extracted over 21,000 abstracts containing NEI, NEI-related CM, and cold–heat patterns-related keywords from the PubMed database and NEI-related genes from the Human Genome Organization (HUGO) glossary. Using the Neato program in the Graphviz software (AT&T; http://www.research.att.com/sw/tools/graphviz/), Li et al. visualized NEI-related gene/CM networks and the CP and HP networks, computed connectivities among the elements, and evaluated connectivities statistically. The results showed that hypothalamus-pituitary-adrenal (HPA) and hypothalamus-pituitary-thyroid (HPT) hormones were significantly related to the CP network, while the cytokines of the immune systems were significantly related to the HP network. Similarly, while HPA and HPT hormone-related genes were found to be predominant in the CP network, immune-related genes were mainly present in the HP network.

While several studies have been undertaken to determine the underlying biological pathways for cold–heat patterns, to date their heritability has been very rarely explored. In a large population-based sample of UK female twins (aged 30–60 years), Cherkas et al. (2007) estimated heritability of Raynaud's phenomenon (RP), a phenotype similar to cold hypersensitivity in the hands and feet. Cherkas et al. reported heritability of RP to be 55%, suggesting that it has a genetic underpinning. In a sample of South Korean twins aged between 12 and 24 years, Hur et al. (2012) found heritability of 64% in males and females for a single question about cold hypersensitivity in the hands and feet. Although these studies employed large population-based twin samples, they estimated heritability either on the basis of a single question or RP, and it is therefore difficult to apply the results to cold–heat patterns that encompass symptoms associated with a broad range of diseases and related risk factors.

In the present study, we report the heritability of cold–heat patterns in adolescent and young adult twins using a questionnaire that includes multiple questions measuring symptoms related to cold–heat patterns. We also examined whether or not heritability differs across sexes and whether it changes with age during adolescence and young adulthood.

Materials and Methods

Sample

The sample included 1,753 South Korean twins unselected for health status. The mean age of the twins was 19.1 ± 3.1 years (range: 12–29 years). Twins under 20 years of age were recruited mostly from schools throughout South Korea, while those aged 20 years or older were recruited from Facebook, twin clubs on the internet, and colleges throughout South Korea.

Table 1 presents demographic characteristics of the sample by zygosity. As with most volunteer twin samples (Heath et al., 1989), the present sample has an over-representation of females (62%). An over-representation of females was also due in part to the fact that some of the male twins were in the military service, as army service is mandatory for young adult males in South Korea. Zygosity of the twins was assessed using a three-item zygosity questionnaire. When compared to DNA analysis, this approach has been shown to achieve over 90% accuracy (Ooki et al., 1993). The number of monozygotic (MZ) twins was much greater than that of dizygotic (DZ) twins in the present sample, which likely reflected the low DZ twin birth rates in the South Korean population for the birth cohorts in the present study rather than sampling bias (Hur and Kwon, 2005). The mean levels of BMI and height of our sample were very close to those of adolescents and young adults in South Korea (KOSIS, 2017; Park, 2011), suggesting that our sample is fairly representative of South Korean adolescents and young adults in terms of physical characteristics.

TABLE 1 Means and Standard Deviations of the Age, Height, and BMI of the Sample by Zygosity

MZM = monozygotic male, DZM = dizygotic male, MZF = monozygotic female, DZF = dizygotic female, OSDZ = opposite-sex dizygotic twins.

Measure

Trained interviewers administered a questionnaire about cold–heat patterns (Yeo et al., 2016) to twins via a telephone interview. The cold–heat patterns questionnaire used in the present study consisted of eight items for typical symptoms of the CP and seven items for typical symptoms of the HP (Table 2). This questionnaire was developed for the general population using factor analysis and clinicians’ diagnostic consensus. The items in the questionnaire have been demonstrated to have adequate psychometric properties (Yeo et al., 2016). Twins were instructed to respond to five-point Likert scale (1 = not at all true, 5 = very true) for each item. Cronbach alpha reliabilities in the present sample were acceptable: They were 0.66 for the CP and 0.70 for the HP.

TABLE 2 Items of the Cold–Heat Patterns Questionnaire Developed by Yeo et al. (2016)

Biometrical Model-Fitting

A standard twin model incorporating age as a moderator (Neale and Cardon, 1992; Purcell, 2002) was applied to the raw data for the CP and the HP, respectively. The full model included additive genetic factors (A) (i.e., the sum of the average effect of all alleles that influence a phenotype), shared environmental factors (C) (i.e., those environmental factors shared between the two members of a twin pair reared in the same home), and unique environmental factors (E), (i.e., those environmental factors unique to each member of a twin pair and measurement error). On the basis of the degree of genetic relatedness, the correlations for A were set at 1.0 for MZ and 0.50 for DZ twins. The correlation for C was set at 1.0 for both types of twins because these twins were raised together in the same home. The correlation for E was set at zero because by definition, this factor was uncorrelated between the two members of a twin pair. The model also includes the standard paths, a, c, and e, indicating the magnitude of moderation effects related to A, C, and E, respectively, which are allowed to vary as a function of sex, and a moderator, age (M). In the model, the total phenotypic variance (V_p) and the covariance for MZ and DZ twins for each sex can be expressed as follows:

$$\begin{eqnarray*} {{\rm{V}}_{\rm{p}}} &=& {({\rm{a}} + {{\rm{\beta }}_{\rm{\alpha }}}M)^2} + {({\rm{c}} + {{\rm{\beta }}_{\rm{c}}}M)^2} + {({\rm{e}} + {{\rm{\beta }}_{\rm{e}}}M)^2}\\ {\rm{CO}}{{\rm{V}}_{{\rm{MZ}}}} &=& {({\rm{a}} + {{\rm{\beta }}_{\rm{\alpha }}}M)^2} + {({\rm{c}} + {{\rm{\beta }}_{\rm{c}}}M)^2}\\ {\rm{CO}}{{\rm{V}}_{{\rm{DZ}}}} &=& 0.5\;{({\rm{a}} + {{\rm{\beta }}_{\rm{\alpha }}}M)^2} + {({\rm{c}} + {{\rm{\beta }}_{\rm{c}}}M)^2} \end{eqnarray*}$$

In the equations above, β_α, β_c, and β_e represent the magnitude of the moderating effects of age associated with A, C, and E, respectively. Thus, if β_α is significantly different from zero, then the influence of additive genetic factors varies as a function of age. The same applies to β_c and β_e.

Model-fitting analyses were performed using the maximum likelihood raw data option in M_x (Neale et al., 2003) that calculates twice the negative log-likelihood (-2 LL) of the data. To determine the best-fitting, most parsimonious model, parameters in the full model are constrained sequentially. The likelihood ratio test (LRT) and Akaike's information criterion (AIC = -2 LL - 2df; Akaike, 1987) were used to evaluate alternative models.

Results

Descriptive Statistics and Twin Correlations

Means, standard deviations, and correlations with age for the CP and the HP are presented in Table 3. For CP, females showed a significant positive correlation with age, while males did not, suggesting that CP increases with age during adolescence and young adulthood in females but not in males. The opposite pattern was observed in the HP: The correlation with age was significant and positive only in males. Note, however, that the magnitudes of correlations with age were relatively modest in both patterns and in both sexes.

TABLE 3 Correlations With Age and BMI, Means, and Standard Deviations for the Cold Pattern (CP) and Heat Pattern (HP) for Males and Females

M = males, F = females.

**Correlation significantly different from zero at p < .01.

††Means and SDs are significantly different between males and females at p < .01.

^a Age was adjusted.

Whereas the HP did not show any significant sex difference in means or standard deviations, CP was significantly higher in females than in males in both means and standard deviations. Age-adjusted correlations of BMI with the CP and the HP were -0.23 (p < .01) for males and -0.19 (p < .01) and 0.16 (p < .01), respectively, for females. These sex and age differences in the CP and the HP and correlations with BMI were consistent with the results from prior studies (Pham et al., 2016; Yoshino et al., 2013).

As a preliminary step, we computed twin correlations across zygosity groups to predict the results from model-fitting analysis. As can be seen in Table 4, for the CP and HP, MZ twin correlations were greater than DZ twin correlations in both sex groups, suggesting the presence of genetic influences. DZ twin correlations were greater than half the MZ twin correlations in females but not in males, suggesting a hint of shared environmental influences in females. These observations of the twin correlations were formally tested using model-fitting analysis.

TABLE 4 Maximum Likelihood Twin Correlations and Their 95% Confidence Interval for the Cold Pattern (CP) and Heat Pattern (HP) by Zygosity

95% CI are in square brackets. MZM = monozygotic male twins, DZM = dizygotic male twins, MZF = monozygotic female twins, DZF = dizygotic female twins, OSDZ= opposite-sex dizygotic twins.

Model-Fitting Analysis

Tables 5 and 6 show the results of the model-fitting analysis for the CP and HP, respectively. From the full model for the CP, we first removed age moderator effects associated with A, C, and E for males and females altogether (β_am, β_af, β_cm, β_cf, β_em, β_ef), which resulted in a non-significant change in chi-square (Model 1), suggesting that genetic and environmental influences do not significantly change with age. Next, in Model 2, we further removed shared environmental factors for each sex (c_m, c_f,) from Model 1, whereas in Model 3, we removed additive genetic factors (a_m, a_f). The model fit was acceptable in Model 2 but not in Model 3, suggesting that shared environmental influences are not significant, whereas additive genetic influences are critically important. Finally, we equated additive genetic and unique environmental factors across males and females (a_m = a_f,, e_m = e_f), which produced a non-significant chi-square difference, suggesting that the magnitudes of additive genetic and unique environmental factors are the same for both sexes. From these model comparisons using LRT, Model 4 was chosen as the best fit. In agreement with the results from LRT, AIC suggested that Model 4 was the best fit because it showed the lowest value.

TABLE 5 Biometrical Model-Fitting Results for the Cold Pattern (CP)

LL = Log likelihood. The best-fitting model is indicated in bold type. a = additive genetic path, c = shared environmental path, e = unique environmental path. β_a = the effect of age moderator associated with ‘a’, β_c = the effect of age moderator associated with ‘c’, β_e = the effect of age moderator associated with ‘e’. Subscripts ‘m’ and ‘f’ denote males and females, respectively.

TABLE 6 Biometrical Model-Fitting Results for the Heat Pattern (HP)

LL = Log likelihood. The best-fitting model is indicated in bold type. a = additive genetic path, c = shared environmental path, e = unique environmental path. β_a = the effect of age moderator associated with ‘a’, β_c = the effect of age moderator associated with ‘c’, β_e = the effect of age moderator associated with ‘e’. Subscripts ‘m’ and ‘f’ denote males and females, respectively.

For the HP in Table 6, elimination of all age moderators (β_am, β_af, β_cm, β_cf, β_em, β_ef) from the full model yielded a borderline significant change in chi-square (Model 1). As with CP, additive genetic effects were significant, but shared environmental effects were not (Model 2 vs. Model 3). Constraining additive genetic and unique environmental effects to be equal between males and females (a_m = a_f,, e_m = e_f) resulted in a non-significant change in Model 4. Thus, we concluded that Model 4 was the best-fitting one for the HP. Model 4 had the lowest AIC for the HP.

In the best-fitting models, heritability estimates were 40% (95% CI [38, 42]) for the CP and 33% (95% CI [25, 42]) for the HP in both males and females. Unique environmental influences including measurement error were 60% (95% CI [58, 62]) for the CP and 67% (95% CI [58, 75]) for the HP in both sexes. These estimates did not change significantly with age during the age period we studied (i.e., 12 years to 29 years).

Discussion

Our study clearly demonstrated that cold–heat patterns have a genetic etiology. These results can facilitate future molecular genetic studies to identify genes associated with cold–heat patterns as well as genes for various diseases related to the two patterns. Genetic influences found in the present study may reflect genes for body temperature regulation and/or genetic susceptibility to various diseases associated with cold–heat patterns.

Interestingly, our heritability estimate for the CP was lower than heritability for the cold hypersensitivity in the hands and feet (64%; 95% CI [55, 72]) reported in Hur et al. (2012) or heritability estimates for RP (55%; 95% CI [40, 70]) or the cold sensitivity (53%; 95% CI [42, 64)] found in Cherkas et al. (2007), with little overlapping of confidence intervals. The discrepancies in findings among the three studies may be due to the differences in measures, samples, or age of the twins. However, the discrepant results also suggest that the CP, the cold hypersensitivity in the hands, and feet, and RP may be overlapping but different phenotypes. Multivariate genetic analysis of the three phenotypes may elucidate the etiological relationships of the three phenotypes.

Our analysis also showed that environmental factors unique to each member of twin pairs were substantial, pointing to the importance of detecting environmental factors that influence cold–heat patterns. Candidate environments that may influence cold–heat patterns include lifestyle factors such as exercise, dietary factors, nutrient intakes, psychological stress, and viruses. Our twin model was based on the assumption that there is no interaction between genetic and unique environmental factors. However, there is increasing recognition in the field of medicine that the effects of genes for complex diseases and disease-related phenotypes may vary across different levels of environmental exposures (Riz et al., 2017). Given that cold–heat patterns are complex, multifactorial phenotypes, efforts should be made in future to unravel how genetic and environmental factors act and interact in the development of cold–heat patterns.

There were a few limitations in the present study that need to be addressed. First, most of our twins had passed puberty. The estimates of heritability and environmental influences on the CP and the HP in prepubescent children or adults after menopause may be different from what we found, given that hormones are known to be involved in molecular pathways of the CP (Ma et al., 2010). To more fully delineate how genetic and environmental influences change with age, future research should broaden the age range of the subjects to include young children, middle-aged, and old-aged adult twins. Second, we did not find significant sex differences in genetic and environmental influences on the CP or the HP. However, there was some indication of shared environmental influences on CP only in females. These shared environmental influences on CP in females were observed in the cold hypersensitivity in the hands and feet as well (Hur et al., 2012). To draw definite conclusions on sex differences in genetic and environmental influences on cold–heat patterns, replication of our findings with larger statistical power is necessary. Finally, in the present study, a dimensional approach was used to measure the CP and the HP, assuming that there is an underlying continuum for symptoms of cold–heat patterns. Caution is needed when our results are generalized to the clinical diagnosis of the cold and heat syndromes where classifications are made by clinicians’ subjective diagnosis or by using cut-off scores.