In traditional East Asian medicine, the cold–heat pattern differentiation system has long been used as a part of clinical diagnosis (Jian, Reference Jian2005). 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, Reference Jian2005). Incorporation of the cold–heat differentiation system in clinical diagnosis has been shown to improve treatment outcomes as well. Lu et al. (Reference Lu, Zha, Chang, He and Lu2009) 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., Reference Bae, Lee, Park, Yoon, Mun and Lee2017).

Recently, there is growing research interest in the biological basis of cold–heat patterns. For instance, van Wietmarschen et al. (Reference van Wietmarschen, Yuan, Lu, Gao, Wang, Xiao and van der Greef2009) 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., Reference Hsu, Yu, Lee, Lee and Yang2003), Li et al. (Reference Li, Zhang, Wu, Zhang, Li and Wang2007) 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., Reference Jenssen, Lagreid, Komorowski and Hovig2001; Li et al., Reference Li, Wu and Zhang2006). 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. (Reference Li, Zhang, Wu, Zhang, Li and Wang2007) 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. (Reference Cherkas, Williams, Carter, Howell, Black, Spector and MacGregor2007) 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. (Reference Hur, Chae, Chung, Kim, Jeong, Kim and Kim2012) 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., Reference Heath, Howells, Kirk, Madden, Bucholz, Nelson and Martin1989), 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., Reference Ooki, Yamada and Asaka1993). 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, Reference Hur and Kwon2005). 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, Reference Park2011), 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., Reference Yeo, Park, Bae, Jang and Lee2016) 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., Reference Yeo, Park, Bae, Jang and Lee2016). 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. (Reference Yeo, Park, Bae, Jang and Lee2016)

Biometrical Model-Fitting

A standard twin model incorporating age as a moderator (Neale and Cardon, Reference Neale and Cardon1992; Purcell, Reference Purcell2002) 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., Reference Neale, Boker and Xie2003) 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, Reference Akaike1987) 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., Reference Pham, Lee, Kim, Song, Kim and Leem2016; Yoshino et al., Reference Yoshino, Katayama, Munakata, Horiba, Yamaguchi, Imoto and Watanabe2013).

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).