Yasukochi et al.
Journal of Physiological Anthropology(2023) 42:2
https://doi.org/10.1186/s40101-023-00319-2

ORIGINAL ARTICLE

Journal of
Physiological Anthropology

Open Access

Cold‑induced vasodilation response
in a Japanese cohort: insights from cold‑water
immersion and genome‑wide association
studies
Yoshiki Yasukochi1* , Toshihiro Sera2, Taiki Kohno3, Yusuke Nakashima3, Musashi Uesugi4 and Susumu Kudo2*

Abstract
Background Cold-induced vasodilation (CIVD) occurs after blood vessels in the skin are constricted due to local
cold exposure. Although many CIVD studies have been conducted, the underlying molecular mechanisms are yet
to be clarified. Therefore, we explored genetic variants associated with CIVD response using the largest-scale dataset reported to date in a CIVD study involving wavelet analysis; thus, the findings improve our understanding of themolecular mechanisms that regulate the CIVD response.
Methods We performed wavelet analysis of three skin blood flow signals [endothelial nitric oxide (eNO)-independent, eNO-dependent, and neurogenic activities] during finger cold-water immersion at 5 °C in 94 Japanese young
adults. Additionally, we conducted genome-wide association studies of CIVD using saliva samples collected from theparticipants.
Results We found that the mean wavelet amplitudes of eNO-independent and neurogenic activities significantly
increased and decreased prior to CIVD, respectively. Our results also implied that as many as ~ 10% of the Japanese
subjects did not show an apparent CIVD response. Our genome-wide association studies of CIVD using ~ 4,040,000
imputed data found no apparent CIVD-related genetic variants; however, we identified 10 genetic variants, including
2 functional genes (COL4A2 and PRLR) that are associated with notable blunted eNO-independent and neurogenic
activity responses in individuals without CIVD response during local cold exposure.
Conclusions Our findings indicate that individuals without CIVD response differentiated by genotypes with COL4A2
and PRLR genetic variants exhibited notable blunted eNO-independent and neurogenic activity responses duringlocal cold exposure.
Keywords Cold-induced vasodilation, Genome-wide association study, Japonica array, Single-nucleotide
polymorphism, Wavelet analysis*Correspondence:
Yoshiki Yasukochi
yasukocy@hirakata.kmu.ac.jp
Susumu Kudo
kudo.susumu.237@m.kyushu-u.ac.jp
Full list of author information is available at the end of the article© The Author(s) 2023. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​
mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

BackgroundPhysiological thermoregulatory systems that act against
extreme cold are essential for the survival of organisms. For example, constriction of skin vessels occurs in
response to local cold exposure to minimize heat loss inthe human body, resulting in attenuated skin blood flow
(SkBF). After several minutes of cold exposure, peripheral blood flow increases due to subsequent vasodilation,
known as “cold-induced vasodilation” (CIVD) or “hunting reaction” [1]. CIVD is thought to play an important role in biological defense against local cold-related
injuries, such as frostbite of the fingers, although theunderlying mechanisms of CIVD remain unclear [2].
Understanding these mechanisms could help prevent
cold injuries. Interestingly, the CIVD response reportedly
varies among ethnic groups; however, it is difficult to distinguish between the effects of genetic adaptation to cold
environments (i.e., genetically fixed resistance to suchenvironments in a population or ethnic group) and physiological acclimatization caused by cold stress (i.e., plastic
physiological changes in an individual due to chronic or
repeated cold stress) [2]. Although the timing of CIVD
responses estimated using SkBF in the middle finger canbe reproducible [3], the heterogeneity of CIVD responses
among fingers and toes has been reported previously[4]. Such complexity increases the difficultly of elucidating the underlying mechanism of CIVD. In previous
studies, wavelet analysis of SkBF measured using a laser
Doppler velocimeter has been applied to examine CIVD
responses because the wavelet transform can divide themeasured perfusion signals into endothelial, neurogenic,
myogenic, respiratory, and cardiac oscillatory components [5]. Such analysis allows researchers to estimate theactivities of these components via a noninvasive monitoring technique.
Several hypotheses on the physiological mechanismof CIVD have been proposed. It has been suggested
that CIVD is caused by neurogenic activity, such as the
cessation of adrenergic neurotransmission, the releaseof dilating substances [e.g., nitric oxide (NO)] in theendothelium of blood vessels, and myogenic activity, suchas the relaxation of the smooth muscle cells [2, 6–8]. Our
previous study revealed that both neurogenic activity andendothelial NO (eNO)-independent activity may affectthe finger SkBF response during immersion in cold waterat 5 °C [9]. Indeed, it is possible that eNO-independent
activity derived from endothelium-derived hyperpolarizing factor (EDHF) may affect the vasodilation of small
vessels [10, 11]. Therefore, it is important to determine
which molecules involved in eNO-independent and neurogenic activities primarily affect the CIVD response.
To elucidate the underlying physiological mechanismof CIVD, a genotype–phenotype association study is

Page 2 of 16

required because physiological functions are systematically regulated by interactions among various molecules.
Although many CIVD studies have been conducted, theassociation between genetic variants and CIVD has yet
to be revealed. This might be because substantial experimental effort is required to collect wavelet analysis data
on SkBF; indeed, to the best of our knowledge, wavelet
analysis of > 40 individuals has yet to be performed [12].
Here, we conducted such analysis using 94 study subjects.
To quantify the magnitude of CIVD, wavelet analysis of
SkBF signals was implemented to measure the magnitudeof three possible vascular regulatory factors: eNO-independent, eNO-dependent, and neurogenic activities.
As mentioned, whether genetic polymorphism is
attributable to the interethnic variability of the CIVD
response is debatable [2]. To address this issue, we must
explore the association between genetic polymorphisms
and CIVD-related phenotypes to clarify the contribution of genetic variants to CIVD, which will provide crucial insights into the interethnic variability of the CIVD
response. Therefore, we aimed to identify genetic variants related to CIVD and conducted genome-wide association studies (GWASs) on CIVD in 94 young Japanese
adults. Our study is the first in which a finger cold-water
immersion experiment and GWAS were performed, andwe present the largest dataset in the field of CIVD todate. Moreover, our results begin to elucidate the underlying mechanism of CIVD by revealing the molecularfunctions related to the CIVD response.

Methods
Study subjects

In total, 94 healthy Japanese volunteers (93 males and1 female; mean age ± standard error: 22.3 ± 0.1 years;
height: 171.8 ± 0.6 cm; weight: 64.0 ± 1.1 kg) wererecruited from Kyushu University (Fukuoka, Japan). The
participants were nonsmokers and were not taking any
medication. We refer to this cohort as the “JPQ cohort.”
Since wavelet amplitudes of SkBF signals and anthropometric data such as body mass index and weight in the
one female were not deviated from the corresponding
measurements in males (Additional file 1: Table S1), we
used the genotype and measurement data of the female
for further analyses.
Experimental cold‑water immersion of fingers

The attenuated SkBF of the 94 study subjects was measured using laser Doppler flowmetry (FLO-C1, OMEGAWAVE). A probe for SkBF (ML type, OMEGAWAVE) anda thermistor for sensing temperature (TSD202F, BIOPAC
Systems, Inc., USA) were attached to the ventral aspect ofthe distal phalanx of the middle finger at room temperature (25 °C). Before the experiment began, the subjects

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

had rested in the sitting position for at least 30 min until
all parameters to be examined were stable, and they further rested quietly for 10 min to collect baseline data.
Subsequently, the middle finger of the participants was
immersed up to the second phalanx for 30 min in a custom-made circulating cold water bath (5 °C). We used the60-s moving average filter of the measurements for SkBF
and finger skin temperature to avoid small fluctuations.
Finger SkBF was sampled at 200 Hz using the MP150
program (BIOPAC Systems, Inc.). Oscillation data were
downsampled to 5 Hz using AcqKnowledge software(BIOPAC Systems, Inc.), and a Morlet mother wavelet was applied to conduct wavelet analysis of the SkBF
signals using the Biomass program (ELMEC Inc.). The
wavelet amplitudes were presented as a ratio of thetotal power between 0.005 and 0.150 Hz covering thefrequency bands related to eNO-independent (0.005–
0.010 Hz), eNO-dependent (0.01–0.02 Hz), and neurogenic (0.02–0.05 Hz) activities [5, 8, 12, 13]. We defined
a CIVD event according to the definition of Sera et al. [9],which was as follows: “A CIVD event was determined as aminimum increase in finger temperature after the bloodflow reached the minimum first, due to local cooling, andthen increased.” In other words, followed by the reactionthat the blood flow reached the minimum first prior toreascending it, the finger temperature also reached theminimum and then increased. This event was determined as a CIVD. The authors also defined three phases
resulting from marked changes in the SkBF during finger cold exposure: vasoconstriction (phase 1), prior to

Page 3 of 16

CIVD (phase 2), and CIVD (phase 3). These were definedaccording to the timing of the increased rate of SkBF due
to the onset of the CIVD [9]. Figure 1 shows representative finger SkBF and temperature in response to finger
cooling. Phases 1, 2, and 3 represent the periods from the
beginning of cold exposure to minimum blood flow, fromthe minimum blood flow to the onset of CIVD, and fromthe onset of CIVD to maximum blood flow, respectively.
We investigated the wavelet amplitudes at these phases
as well as at the baseline. Detailed experimental and datacollection methods were described previously by Sera
et al. [9].
Genome‑wide genotyping and imputation

An Oragene DNA Kit (DNA Genotek, Ottawa, Canada)
was used to collect ≥ 2 mL of saliva from each subject.
From these samples, genomic DNA extraction and genotyping were conducted using an Axiom Japonica array(Toshiba, Tokyo, Japan) [14] by Cell Innovator Co., Ltd(Fukuoka, Japan), resulting in ~ 660,000 genetic variants among 94 study subjects, with a < 3% missing rateof the genotype data. Subsequently, we filtered out thegenetic variants described below using PLINK ver. 1.90
[15]. After checking gender matching between phenotypic records and heterozygosity of genetic variants on
chromosome X and cryptic relatedness between DNA
samples (proportion identity with a descent threshold of ≥ 0.1875), we performed quality control (QC) of
the genotype data and discarded the following variants:
variants with a call rate of < 0.97, variants for which theFig. 1 Representative blood flow and temperature in response to finger cooling for individuals with CIVD response at room temperature (25 °C).
This figure is comparable to Fig. 1 in our previous study [9]. White and black dots indicate blood flow and temperature, respectively. The finger was
exposed to cold water at 5 °C for 30 min following 10 min of rest (Baseline). The first minimum and maximum temperatures are indicated by blackarrows and arrowheads, respectively, suggesting the presence of a CIVD event. The first maximum and minimum blood flows are indicated by whitearrowheads and arrows, respectively. The onset time of CIVD (CIVD onset_sbf ) is indicated by small arrows. Phase 1 (vasoconstriction; from the onsetof cooling to minimum blood flow), phase 2 (prior to CIVD; from minimum blood flow to CIVD onset_sbf ), and phase 3 (CIVD; from CIVD onset_sbf
to maximum blood flow) are indicated in parentheses as 1–3, respectively. AU, arbitrary unitsYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

genotype distribution significantly (p < 0.0001) deviatedfrom the Hardy–Weinberg equilibrium, variants with aminor allele frequency (MAF) of < 0.005, and variants onsex chromosomes and the mitochondrial genome. Consequently, 622,446 genetic variants were used for genotype imputation.
After prephasing the genotype data that passed QC
via SHAPEIT v2.r904 [16, 17], untyped genotype data
were imputed with the 1000 Genomes Project (1KGP)
reference panel (phase 3) using IMPUTE2 ver. 2.3.2 [18](Ne = 20,000; chunk size = 5 Mb). Genetic variants with
low imputation quality (info score of < 0.5) were discarded when the imputed output data were converted
into PLINK format data using GTOOL ver. 0.7.5 (https://​www. ​ w ell. ​ ox. ​ a c. ​ u k/ ​ ~ cfree ​ m an/ ​ s oftw ​ a re/ ​ g was/​ g tool.​html). We also conducted QC for the imputed genotype
data in the same manner, although a 1% MAF threshold
was set for the imputed data.
To assess population stratification on a genome-wide
scale, we performed principal component analysis (PCA)
using the Japonica array genotype data that passed QC
in the JPQ cohort after removing genetic variants with a
MAF of < 0.05 as well as datasets of five East Asian populations [CDX (Chinese Dai in Xishuangbanna, China),
CHB (Han Chinese in Beijing, China), CHS (Southern Han Chinese), JPT (Japanese in Tokyo, Japan), andKHV (Kinh in Ho Chi Minh City, Vietnam)] retrieved
from 1KGP [19] phase 3 reference panels (NCBI Build
GRCh37; http://​bochet.​gcc.​biost​at.​washi​ngton.​edu/​beagle/​1000_​Genom​es_​phase3_​v5a). The datasets of 1KGP
populations were generated from variant call format(VCF) data using the PLINK program after the VCF data
were converted to the binary VCF data using SAMtools/
BCFtools ver. 1.9 [20]. To detect possible population outliers, PCA was also conducted using the genotype data of
the JPQ cohort only.
Construction of a network of CIVD‑related categorical data

We focused on two SkBF signals, eNO-independent andneurogenic activities, as these factors can largely affectthe finger SkBF response during local cold exposure [9].
Based on the results in our finger cold-water immersion
experiment, we categorized the physiological responsesof the study subjects in the following manner. First, thesubjects were divided into two groups with either thepresence or absence of a CIVD response. Second, thesubjects were assigned to two additional groups in which
they showed either higher or lower amplitudes of each
SkBF signal (i.e., eNO-independent and neurogenic
activities) at phase 2 compared with the respective mean
wavelet amplitudes in individuals who did not show anapparent CIVD response. Consequently, we generated
a binary dataset of three categorical variables: (i) the

Page 4 of 16

presence or absence of a CIVD response and higherthan-average or lower-than-average wavelet amplitudesof (ii) eNO-independent and (iii) neurogenic activities.
Using the dataset, we constructed a network of the three
CIVD-related categorical variables using the medianjoining method [21] via NETWORK 10.1.0.0 (https://​www.​fluxus-​engin​eering.​com; Fluxus Technology Ltd.,Colchester, England).
Genome‑wide association studies of CIVD

We tested the significance of the association between
the three aforementioned CIVD-related categorical datatypes and single-nucleotide polymorphisms (SNPs) using
logistic regression analysis and Fisher’s exact test with1,000,000 max(T) permutations in the additive model
via the PLINK program. In addition, two-way analysis ofvariance (ANOVA) was conducted using CIVD response(presence or absence) and SNP genotypes as independent variables as well as the measurements of SkBF signals
related to eNO-independent and neurogenic activities atphase 2, at which stage the amplitudes were remarkably
changed [9], as dependent variables. To conduct two-way
ANOVA, PLINK format data were converted into thesuitable format for the analysis using the Bioconductor
R package snpStats ver. 1.20.0 [22] via R ver. 3.4.2 [23].
Furthermore, the genotype data for all subjects wereconverted into numeric data with dominant and additive genetic inheritance models. Two-way ANOVA was
conducted using R ver. 3.5.3 [24] in RStudio ver. 1.2.5042[25].
SNPs affecting CIVD are expected to be shared
between GWASs of phenotypes derived from SkBF signals. Moreover, SkBF is likely controlled by various
genetic variants (i.e., polygenic inheritance or epistasis)
because quantitative traits related to vascular function
and diseases, such as peripheral artery disease and hypertension, can be affected by two or more genes [26–29];
thus, the effect of each SNP on CIVD may not be large.
Therefore, we identified candidate CIVD-related genetic
variants as follows: SNPs associated with the two CIVDrelated SkBF signals or three CIVD-related categorical
data (see “Construction of a network of CIVD-related
categorical data”) were obtained with a p value threshold of < 1 × ­10−5, which was considered a borderline significant level. We then searched SNPs for which the pvalue was < 5 × ­10−8 at a genome-wide level in GWASs
for either eNO-independent or neurogenic activities. Of
the SNPs with p values at this significant level, we identified SNPs associated with the other activity at a borderline significant level as candidate CIVD-related genetic
variants.
Using a generalized estimating equation (GEE)model [30, 31] via the R package “geepack” [32], weYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

also examined the association between repeated measurements of eNO-independent or neurogenic activities across all phases and combined phenotypes of SNP
genotypes with CIVD response (presence or absence).
Because the wavelet amplitudes were continuous data, aGaussian distribution for the family with an identity link
was applied to assess the association in the GEE method.
The working correlation matrix (correlation of repeatedmeasurements within an individual) was set to a firstorder autoregressive model, and the wave argument wasused to specify the ordering of repeated measurements
within individuals.
Manhattan plots of p values were produced using theR package “CMplot” (https://​github.​com/​YinLi​Lin/R-​
CMplot). Linkage disequilibrium (LD) between pairs of
SNPs in the JPQ cohort was estimated using Haploview
version 4.2 [33], whereas the LD in JPT from the 1KGP
database was assessed using LDproxy in the LDlink webbased tool (https://​ldlink.​nci.​nih.​gov/) [34]. Informationon the association of genetic variants with vascular function or disease was searched in the Genome-Wide Repository of Associations Between SNPs and Phenotypes

Page 5 of 16(GRASP, https://​grasp.​nhlbi.​nih.​gov/​Overv​iew.​aspx) [35]and GWAS Catalog (https://​www.​ebi.​ac.​uk/​gwas/​home)[36] databases. To predict the effects of SNPs on proteinfunction, Combined Annotation Dependent Depletion(CADD) scores (https://​cadd.​gs.​washi​ngton.​edu/) [37]were used. The relevant biological pathways of genes
containing focal SNPs were estimated using the Kyoto
Encyclopedia of Genes and Genomes (KEGG) database(https://​www.​genome.​jp/​kegg/).

Results
Physiological responses to cold‑water immersion

In total, 94 healthy Japanese volunteers (93 men and1 woman), referred to here as the “JPQ cohort,” wererecruited from Kyushu University. The characteristics of the 94 subjects in this cohort are summarized inTable 1. In the finger cold-water immersion experiment(see “Methods”), wavelet analyses of SkBF signals wereconducted to quantify the activities of the following
vascular regulatory factors: eNO-independent, eNOdependent, and neurogenic activities. Interestingly, noapparent CIVD response was observed in 10 subjects

Table 1 Characteristics of study subjects with or without CIVD response
Characteristic

NonCIVDCIVD

No. of subjects

1084Age, years

22.80 ± 0.44

22.24 ± 0.15

Height, cm
Weight, kg
Body mass index, kg/m2
Body temperature, °C
Minimum finger temperature during local cooling, °C
Mean finger temperature between 15 and 35 min, °C
Time for the first minimum blood flow, ­mina
Neurogenic ­activityb (Baseline)
Neurogenic ­activityb (Phase 1)c
Neurogenic ­activityb (Phase 2)c
Neurogenic ­activityb (Phase 3)c
eNO-dependent activity b (Baseline)
eNO-dependent activity b (Phase 1) c
eNO-dependent activity b (Phase 2) c
eNO-dependent activity b (Phase 3) c
eNO-independent activity b (Baseline)
eNO-independent activity b (Phase 1) c
eNO-independent activity b (Phase 2) c
eNO-independent activity b (Phase 3) c
Data are means ± standard error of the mean

Abbreviations: CIVD Cold-induced vasodilation, eNO endothelial nitric oxideaTime corresponding to the white arrow in Figs. 1 and 3bWavelet amplitude estimated from finger skin blood flow (arbitrary units)cPhases 1, 2, and 3 represent “vasoconstriction,” “prior to CIVD,” and “CIVD,” respectively [9]171.30 ± 2.30
62.22 ± 2.73

21.24 ± 0.93

36.61 ± 0.13

12.9 ± 2.0

16.6 ± 2.3

171.83 ± 0.60
64.25 ± 1.20

21.73 ± 0.38

36.45 ± 0.04
10.3 ± 0.4

13.6 ± 0.4

11.9 ± 0.3

11.7 ± 0.1

0.24 ± 0.05

0.19 ± 0.01

0.20 ± 0.04

0.14 ± 0.01

0.35 ± 0.03

0.14 ± 0.05
0.25 ± 0.03

0.31 ± 0.01

0.06 ± 0.01
0.25 ± 0.01

0.27 ± 0.04

0.26 ± 0.01

0.20 ± 0.04

0.17 ± 0.01

0.20 ± 0.04
0.24 ± 0.04

0.14 ± 0.01
0.30 ± 0.01

0.38 ± 0.07

0.48 ± 0.02

0.45 ± 0.08

0.58 ± 0.02

0.58 ± 0.10

0.76 ± 0.02

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2(~ 10%), whereas the remaining 84 subjects exhibited atleast one CIVD response. In the 84 subjects with CIVD
responses, the patterns of wavelet amplitudes across thefour phases [baseline, vasoconstriction (phase 1), prior
to CIVD (phase 2), and CIVD (phase 3); Fig. 2a] were inagreement with those reported in our previous study [9].
The amplitudes of eNO-independent activity increased
until phase 2, after which they decreased; conversely,
the amplitudes of neurogenic activity decreased until
phase 2, after which they increased. The differences in
the wavelet amplitudes of eNO-independent and neurogenic activities between the baseline and phases 1–3
were statistically significant (p < 0.05, Dunnett’s test).
For eNO-dependent activity, the amplitudes between
the baseline and phase 1 did not differ significantly(p = 0.734), whereas the amplitudes between the baseline and phases 2 and 3 differed significantly (p < 0.0001).
The representative patterns of finger SkBF and temperature in response to finger cooling for individuals with

Page 6 of 16

CIVD response and without apparent CIVD response areshown in Figs. 1 and 3, respectively.
In the 10 subjects who did not exhibit a CIVD
response, there were no significant differences in eNOdependent activity between the baseline and phases 1–3(p = 0.6674–0.9149), but eNO-independent and neurogenic activities between the baseline and phase 2 or 3
(p ≤ 0.0018) differed significantly (Fig. 2b). We postulatethat eNO-dependent activity is not a major factor in the
CIVD response because the mean values of the amplitude at phase 2 were lower than those at the baseline(Fig. 2). In addition, the decreased activity of eNO, which
is a vasodilator, at phase 2 is likely related to reduced
vasodilation, which would indicate that eNO-dependent
activity does not primarily influence CIVD. The wavelet amplitudes of each activity in all study subjects wereshown in Fig. 2c. Overall, wavelet analysis indicated that
eNO-independent and neurogenic activities potentiallycontribute to the CIVD response.

Fig. 2 Wavelet amplitudes of endothelial NO-independent, endothelial NO-dependent, and neurogenic activities in four phases, including thebaseline phase, in the study subjects. These subjects showed the presence (a) or absence (b) of a CIVD response (a n = 84; b n = 10). All study
subjects (n = 94) are shown in c. Phases 1, 2, and 3 represent “vasoconstriction,” “prior to CIVD,” and “CIVD,” respectively [9]. Data are means ± standard
error of the mean. * p < 0.05 and ** p < 0.01 (compared with the baseline using Dunnett’s test). a.u., arbitrary unitsYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Page 7 of 16

Fig. 3 Representative blood flow and temperature in response to finger cooling for individuals without CIVD response at room temperature (25 °C).
White and black dots indicate blood flow and temperature, respectively. The finger was exposed to cold water at 5 °C for 30 min following 10 min ofrest (Baseline). The first maximum and minimum blood flows are indicated by white arrowheads and arrows, respectively. AU., arbitrary units

Genotype data in the JPQ cohort

In the JPQ cohort, ~ 660,000 genetic variants among
the 94 subjects were genotyped using an Axiom Japonica array. Following QC of the genotype data, 622,446
genetic variants were generated. Population stratification was assessed on a genome-wide scale using the
genotype datasets of the JPQ cohort as well as five East
Asian populations (CDX, CHB, CHS, JPT, and KHV)
from the 1KGP database after removing genetic variants
with a MAF of < 0.05 (Fig. 4a). A PCA plot showed thatthe eigenvectors of the JPQ cohort and JPT (consisting ofJapanese individuals) formed a distinct cluster from thoseof the other East Asian populations. The JPQ genotype
data were also used in isolation to examine population
outliers in the PCA, and such outliers were not detectedin this analysis (Fig. 4b); thus, we used all JPQ cohort
genotype data for genotype imputation. After imputation and QC of the genotype data (see “Methods” for thedetailed methods of both), 4,661,825 SNPs among the 94
Japanese subjects were used for GWASs of CIVD-related
traits.
Association of SNPs with CIVD‑related traits in the JPQ
cohortIn GWASs, we focused on eNO-independent and neurogenic activities because these factors likely affect the
CIVD response. However, in GWASs in which logistic regression analysis and Fisher’s exact test were conducted using the three CIVD-related categorical data,
there was a lack of SNPs associated with empirical p
values of < 5 × ­10−8 at the genome-wide significant level
although 187 and 71 genetic variants showing a suggestive association (p < 1 × ­10−5) were detected in the logistic

regression analysis and Fisher’s exact test, respectively(Additional file 1: Tables S2 − 3).
Next, we conducted two-way ANOVA for the GWASs
of eNO-independent and neurogenic activities in the
additive and dominant genetic models to test for associations between these activities and the interaction
between CIVD response and SNPs. Manhattan plots of
GWASs for eNO-independent and neurogenic activities are shown in Fig. 5. We detected 10 SNPs associatedwith CIVD-related phenotypes (Table 2 and Additional
file 1: Tables S4 − 6). The GWASs with two-way ANOVA
revealed that interactions between CIVD response andfive SNPs [rs17113836 of LINC01470 (long intergenic
nonprotein coding RNA 1470) or four SNPs of COL4A2(collagen type IV alpha 2 chain)] were associated withneurogenic activity at phase 2 at the genome-wide significant level (p value range 1.89 × ­10−8 to 3.24 × ­10−8) in the
additive model (Table 2 and Additional file 1: Fig. S1a).
These SNPs were also associated with eNO-independent activity at phase 2 at a borderline significant level (pvalue range 7.02 × ­10−7 to 9.36 × ­10−6). In the dominantmodel, five SNPs [rs931740 of PRLR (prolactin receptor)
and four SNPs at 14q23.3] met the criteria mentioned in“Methods” (p value ranges 1.03 × ­10−8 to 4.78 × ­10−8 forneurogenic activity; 3.29 × ­10−6 to 9.39 × ­10−6 for eNOindependent activity) (Table 2 and Additional file 1: Fig.
S1b).
We examined the wavelet amplitudes of eNO-independent and neurogenic activities across all phases(repeated measurements) in groups of subjects differentiated by CIVD response and rs73619297 of COL4A2(Fig. 6a) or rs931740 of PRLR (Fig. 6b). Interestingly,
the homozygotes of the major alleles in subjects whoYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Page 8 of 16

Fig. 4 Principal component analysis of genotype data from East Asian populations (a) and the JPQ cohort (b). Data were plotted according to
the first (horizontal axis) and second (vertical axis) principal components. CDX: Chinese Dai in Xishuangbanna, China; CHB: Han Chinese in Beijing,
China; CHS: Southern Han Chinese; JPQ: Japanese in Kyushu, Japan; JPT: Japanese in Tokyo, Japan; KHV: Kinh in Ho Chi Minh City, Vietnam

had no apparent CIVD response (NonCIVD–G/G for
rs73619297 and NonCIVD–T/T for rs931740) had thelowest mean amplitudes of eNO-independent activity across the phases, whereas the mean amplitudes ofneurogenic activity in these homozygotes were consistently the highest amplitudes. Using a GEE model, we
also tested the significance of differences between therepeated measurement data of wavelet amplitudes insubjects with the NonCIVD–G/G or NonCIVD–T/T
and other CIVD-genotype categories for each SNP. To

minimize effects of possible population stratification, weused PC1 and PC2 estimated from PCA analysis based on
the imputed data (Additional file 1: Fig. S2) as covariates.
GEE analyses with adjustment for PC1 and PC2 indicatedthat, in the additive model, the eNO-independent activity of subjects with the NonCIVD–G/G for rs73619297
differed significantly from that of the remaining subjects(p = 0.0056), whereas such a difference was not detected
for neurogenic activity (p = 0.1955). The eNO-independent and neurogenic activities differed significantlyYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Neurogenic activity

Page 9 of 16

Endothelial NO-independent activity

4 SNPs of rs931740
14q23.3of PRLR

rs931740
4 SNPs ofof PRLR
14q23.3
4 SNPs of
rs17113836
COL4A2of LINC01470

4 SNPs of
COL4A2
rs17113836 of
LINC01470

Fig. 5 Manhattan plots of p values [plotted as − ­log10(p)] in genome-wide association studies of CIVD responses with the additive model

Table 2 Summary of SNPs showing the significant association of CIVD with neurogenic and endothelial NO-independent activities
according to two-way ANOVA
Genetic model

RefSNP ID

Positiona

Gene

Allele frequency

CADD score

Neurogenic
F

Additive

Dominant

eNO-independent

P

F

P

−8

28.6

7.02 × ­10−7

rs17113836

5: 152,344,712

LINC01470

A: 0.72
T: 0.280.3

38.3

1.89 × ­10

rs73619297

13: 110,967,247

COL4A2

G: 0.62
A: 0.38

0.7

21.1

3.24 × ­10−8

13.2

9.36 × ­10−6

rs73619299

13: 110,967,456

COL4A2

C: 0.62
T: 0.38

2.8

21.1

3.24 × ­10−8

13.2

9.36 × ­10−6

rs7321461

13: 110,967,990

COL4A2

A: 0.62
C: 0.38

0.1

21.1

3.24 × ­10−8

13.2

9.36 × ­10−6

rs3832899

13: 110,968,318

COL4A2

T: 0.62–: 0.38

1.5

21.1

3.24 × ­10−8

13.2

9.36 × ­10−6

rs931740

5: 35,222,040

PRLR

T: 0.59
C: 0.415.3

39.9

1.03 × ­10−8

24.6

3.29 × ­10−6

rs10144241

14: 65,638,439

14q23.3

C: 0.57
A: 0.43

11.3

35.6

4.78 × ­10−8

22.1

9.39 × ­10−6

rs10133408

14: 65,638,443

14q23.3

T: 0.57
G: 0.43

14.1

35.6

4.78 × ­10−8

22.1

9.39 × ­10−6

rs11850650

14: 65,639,582

14q23.3

T: 0.57
C: 0.43

7.4

35.6

4.78 × ­10−8

22.1

9.39 × ­10−6

rs11158583

14: 65,641,686

14q23.3

G: 0.57
C: 0.43

9.3

35.6

4.78 × ­10−8

22.1

9.39 × ­10−6

Abbreviations: CADD Combined Annotation Dependent Depletion, CIVD Cold-induced vasodilation, eNO endothelial nitric oxide, GWAS Genome-wide association
study, SNP Single-nucleotide polymorphisma

Chromosomal position in NCBI build GRCh37

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Page 10 of 16

Fig. 6 Wavelet amplitudes of neurogenic and endothelial NO-independent activities across four phases, including the baseline phase, in groupsof subjects (n = 94). These subjects were differentiated by CIVD response and the genotypes of rs73619297 of COL4A2 (a) or rs931740 of PRLR (b).
Phases 1, 2, and 3 represent “vasoconstriction,” “prior to CIVD,” and “CIVD,” respectively [9]. Data are means ± standard error of the mean. a.u., arbitrary
units

between subjects with NonCIVD–T/T for rs931740
and the remaining subjects in the dominant model(p = 7.3 × ­10−5 for eNO-independent activity; p = 0.0082
for neurogenic activity).
Network analysis of CIVD‑related categorical data

Classification according to CIVD response or wavelet
amplitudes of eNO-independent and neurogenic activities at phase 2 (see “Methods”) assigned study subjects
to 20 and 74 subjects with lower and higher eNO-independent activities, respectively, and 78 and 16 subjects
with lower and higher neurogenic activities, respectively.
We then constructed a network of categorical groups tovisualize the relationships among the three CIVD-related
factors (Fig. 7). Network analysis displayed that the study
subjects could be assigned to seven groups. The maingroup consisted of 72 individuals who exhibited a CIVD
response, had higher-than-average eNO-independent
activity, and had lower-than-average neurogenic activity.
The characteristics related to this group are referred tohereafter as the “typical CIVD phenotype.” Contrastingly,
four individuals showed the opposite pattern to that in
the main group, i.e., the absence of a CIVD response,
lower-than-average eNO-independent activity, andhigher-than-average neurogenic activity. The characteristics of this group are referred to as the “typical NonCIVD
phenotype” (Fig. 7).
We compared the relationships between the genotypesof the 10 candidate SNPs and the network of categorical
data (Additional file 1: Tables S7 and S8). It was difficult
to define the 72 subjects with the typical CIVD phenotype by the genotypes of candidate SNPs only. However, with the exception of one individual, all subjects
with the typical NonCIVD phenotype were homozygotes of major alleles at the candidate SNPs (Additional
file 1: Tables S7 and S8). As mentioned, the amplitudesof eNO-independent and neurogenic activities in subjects with NonCIVD–T/T for rs931740 and NonCIVD–
G/G for rs73619297 remained low and high, respectively,
across all phases relative to the amplitudes in the otherYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Page 11 of 16

Fig. 7 Network of CIVD-related factors (CIVD response, endothelial NO-independent activity, and neurogenic activity). “Presence” and “Absence”
indicate the presence and absence of CIVD response, respectively. “Low” and “High” indicate lower-than-average and higher-than-average
wavelet amplitudes, respectively, of endothelial NO-independent or neurogenic activities at phase 2. “Typical CIVD phenotype” group consistedof individuals with a CIVD response, lower-than-average neurogenic activity, and higher-than-average eNO-independent activity. “Typical
NonCIVD phenotype” group consisted of individuals without a CIVD response, higher-than-average neurogenic activity, and lower-than-average
eNO-independent activity

CIVD-genotype categorical groups (Fig. 6). Thus, theirmajor alleles may partly affect the blunted response ofeNO-independent and neurogenic activities.
In the 10 subjects without a CIVD response, there wasless fluctuation in the wavelet amplitudes of eNO-independent and neurogenic activities from the baseline tophase 3 in the 4 subjects with the typical NonCIVD phenotype relative to the corresponding amplitudes in theother subjects (Additional file 1: Fig. S3). This suggeststhat the response of skin blood vessels (i.e., vasoconstriction or vasodilation) in these four subjects is less susceptible to cold stimuli. We hypothesized that the blunted
response of blood vessels resulted in the maintenance ofskin temperature and tested this hypothesis by examining the minimum skin temperature of study subjects ina cold-water immersion experiment. The mean valuesof minimum skin temperature in subjects with the presence or absence of a CIVD response were 10.3 °C ± 0.4 °C
and 12.9 °C ± 2.0 °C, respectively. These values are theminimum finger temperatures in response to finger
cooling, and thereafter, the temperature increases as the
blood flow increases when CIVD is occurred (Fig. 1).
In the four individuals with the typical NonCIVD phenotype, the mean value of minimum skin temperature,

18.0 °C ± 3.3 °C, was substantially higher than that in theother individuals. The results suggest that the sensitivityof vasoconstriction and vasodilation to local cold stressare low in individuals with the typical NonCIVD phenotype because they can maintain the minimum skin temperature at relatively high levels, compared to individuals
with the other CIVD-related phenotypes, during finger
cold-water immersion.
Linkage disequilibrium analysis and functional predictionof candidate genes

We surveyed the LD between the candidate SNPs andtheir neighboring SNPs using LDproxy in the LDlink
database. According to this database, the candidate SNPs
were in LD with several SNPs in JPT (r2 > 0.8; Additional
file 1: Table S9): Two SNPs and one SNP were in LD with
rs931740 of PRLR and rs17113836 of LINC01470, respectively, whereas rs73619297 of COL4A2 and rs10144241
at 14q23.3 were in LD with six and twenty-three SNPs,
respectively. The candidate SNPs and the adjacent SNPs,
which were in LD with the query SNPs, were not locatedin the exonic region. In the JPQ cohort, each of the fourcandidate SNPs in COL4A2 and 14q23.3 were also instrong LD (D´ = 1.00; r2 = 1.00; Additional file 1: Fig. S4).

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

These results suggest that the associations of the focal
SNPs in COL4A2 or 14q23.3 with the CIVD-related traits
are not independent. According to the GWAS Catalog
and GRASP databases, no SNPs were significantly associated with vascular function or disease. Therefore, itremains difficult to determine the SNPs that are responsible for the CIVD response.
We also examined the CADD-scaled C-scores of candidate SNPs associated with CIVD to predict the effectsof the SNPs on protein function (Table 2). The candidate
SNPs had scores < 15 (CADD-scaled C-score range 0.1–
14.1), and the scores of PRLR and COL4A2 SNPs werelow (CADD-scaled C-score range 0.1–5.3), suggestingthat the focal SNPs do not have strong effects on proteinfunction.
Finally, we predicted the biological pathways of PRLR
and COL4A2 using the KEGG database. Accordingly,
PRLR and COL4A2 were involved in 5 and 10 KEGG
pathways, respectively (Table 3). KEGG analysis also
indicated that these genes were involved in the PI3K/Akt
signaling pathway [38, 39], which plays important roles invascular homeostasis and angiogenesis [40].

Discussion
In our previous study, we found that the eNO-independent activity in 18 subjects significantly increased at vasoconstriction (phase 1) and prior to CIVD (phase 2) [9].
To verify this finding, we performed the replication study
using a larger sample size (n = 94). In a finger cold-water
immersion experiment, we replicated the significant
increase in eNO-independent activity at phases 1 and 2.
We also replicated the significant decrease in neurogenic
Table 3 KEGG pathways of PRLR and COL4A2
Gene name

KEGG ID

Pathway

PRLR

hsa04060

Cytokine-cytokine receptor interaction

hsa04080

Neuroactive ligand-receptor interaction

hsa04151

PI3K-Akt signaling pathway

hsa04630

JAK-STAT signaling pathway

hsa04917

Prolactin signaling pathway

hsa04151

PI3K-Akt signaling pathway

hsa04510

Focal adhesion

hsa04512

ECM-receptor interaction

hsa04926

Relaxin signaling pathway

hsa04933

AGE-RAGE signaling pathway in diabetic complications

hsa04974

Protein digestion and absorption

hsa05146

Amoebiasis

hsa05165

Human papillomavirus infection

hsa05200

Pathways in cancer

hsa05222

Small cell lung cancer

COL4A2

Page 12 of 16

activity at these phases, as was observed in the previous
study [9]. These results indicate that eNO-independent and neurogenic activities contribute to the CIVD
response during finger cold exposure. Furthermore, we
found that 10 of 94 study subjects had no apparent CIVD
response.
In humans, NO reportedly plays an important role inthe dilation of skin blood vessels induced by local warming of the skin [41–43]. However, the cold-water immersion experiments in our previous [9] and present studies
indicated that eNO-dependent activity may not have asubstantial effect on the CIVD response. Given that NOdependent vasodilation may act primarily on the superior mesenteric arteries (~ 650 µm in diameter) in rats
whereas NO-independent vasodilation may be relevant
to small arteries, e.g., resistance arteries (~ 200 µm in
diameter) [10], we previously hypothesized that CIVD islargely unaffected by eNO-dependent activity [9] because
the inner diameter of human arteriovenous anastomoses (AVAs), which may play an important role in CIVD,
is relatively small (~ 10–150 µm) [44]. Although a diameter of capillaries is smaller than that of AVAs, capillary
generally consists of a thin wall of endothelial cells andthe contribution to blood flow regulation is lower. Therefore, the relaxation of AVAs by eNO-independent activity, rather than eNO-dependent activity, is more likely tocontribute to CIVD.
In the present study, we detected 10 candidate SNPs in
GWASs of CIVD using two-way ANOVA. These SNPs
were not shared with genetic variants moderately associated with three CIVD-related categorical data in logistic
regression analyses or Fisher’s exact test (p < 1 × ­10−5).
This inconsistency might be due to the difference of statistical method used. In addition, given that the effectof each SNP on CIVD is not large, larger sample size
which can increase statistical power is required to detecta more definitive association between SNPs and CIVD
response. One and four of the 10 candidate SNPs identified by the GWASs with two-way ANOVA are locatedin PRLR and COL4A2, respectively, whereas rs17113836
is in LINC01470. Because LINC01470 is not well characterized, the functional relevance of this long noncoding
RNA to CIVD could not be inferred. Regarding the SNPsof the functional genes, we found that the study subjects
with NonCIVD–T/T for PRLR rs931740 or NonCIVD–
G/G for COL4A2 rs73619297 showed blunted changes
in eNO-independent and neurogenic activities during
finger cold exposure. Notably, the main group also contained the homozygotes of the major alleles of rs931740
and rs73619297. Thus, we can hypothesize that the major
alleles of PRLR and COL4A2 SNPs have effects in concert with other genetic variants across the genome on the
blunted CIVD-related activities, although the individualYasukochi et al. Journal of Physiological Anthropology(2023) 42:2

effects of each SNP were not large. Alternatively, epigenetics might drive the phenotypic plasticity in response
to local cold exposure. Functional analysis or an epigenome-wide association study will help validate theseresults.
PRLR encodes a receptor belonging to a type I cytokine
receptor superfamily that binds prolactin, which is ananterior pituitary hormone. Prolactin and its receptor (PRLR) can regulate a wide variety of physiological functions in the nervous system [45]. Additionally,
prolactin can suppress the expression of inducible NO
synthase (NOS) induced by proinflammatory cytokines[46]. According to our KEGG analysis, PRLR is involvedin the JAK/STAT signaling pathway, which is known to
be induced by JAK2 via a prolactin–PRLR interaction[47, 48]. In addition, the JAK2/STAT3 signaling pathway
can regulate tumor angiogenesis through modulation ofseveral angiogenesis-related genes [49, 50]. Given thatthe interaction between prolactin and PRLR molecules
affects regulation of NO-independent vasodilation in rats[51], PRLR might contribute to the development and vasodilation of new or existing vessels in response to finger
cold exposure.
COL4A2 encodes the alpha-2 chain of type IV collagen(collagen IV), which is a main component of the basement membrane (BM) and essential for BM stability [52].
Col4a1/Col4a2 knockout mice reportedly die due to various deficiencies, including neuronal ectopias and misarrangement of the capillary network during angiogenesis[52, 53]. Given that (i) the scaffolding of components,
including collagen IV, in the BM plays a key role in vascular homeostasis and development [54] and (ii) familial mutations of COL4A2 can be risk factor for cerebral
abnormalities, such as porencephaly and small-vessel disease [55], the four candidate SNPs of COL4A2 identified
in the present study might affect vascular dynamics during finger cold-water immersion through changes in theaffinity for COL4A1 (collagen type IV alpha 1 chain) in
heterotrimer formation.
According to a KEGG analysis, PRLR and COL4A2
are related to the PI3K/Akt signaling pathway [38, 39],
which is known to promote production of NO, a vasodilator, through phosphorylation of endothelial NOS inendothelial cells [40, 56, 57]. However, NO is a potent
vasodilator of relatively large vessels [10, 11], and eNOdependent vasodilation may not be a key factor in the
CIVD response [9]. Endothelium-derived hydrogen
peroxide has been reported as an EDHF that mediates vasodilation of small resistance vessels in humans[58, 59], and EDHF-mediated relaxation in response toacetylcholine can occur in NOSs (neuronal, inducible,
and endothelial NOSs) in a system-dependent manner[60]. Although the relationship between the PI3K/Akt

Page 13 of 16

signaling pathway and EDHF has yet to be elucidated,
this pathway could potentially be involved in not only
NO-dependent vasodilation but also EDHF-mediatedvascular relaxation [61, 62]. Thus, we can hypothesizethat PRLR and COL4A2 contribute to CIVD throughthe relaxation of vascular smooth muscle cells mediatedby EDHF [11]. Although further studies may reveal the
relationship between EDHF-mediated vasodilation and
CIVD, we cannot currently rule out the possibility that
PRLR and COL4A2 participate in novel CIVD-related
biological pathways.
We acknowledge several limitations in the present
study. First, the sample size is not sufficiently large toshow a significant association between SNPs and CIVDrelated traits with high statistical power. This was due tothe difficulty in sampling from a large number of subjects
for GWASs; however, our study, including a finger coldwater immersion experiment, had the largest sample sizein the CIVD research field to date. To validate the candidate genetic variants related to CIVD identified in our
GWASs, replication GWASs in other Japanese cohorts orother ethnic groups are required. Second, the functional
relevance of the candidate SNPs or genes in relation to
CIVD remains unclear. Lastly, the candidate SNPs identified could be in LD with neighboring genetic variants;
thus, we cannot rule out the possibility that such neighboring genetic variants may have effects on CIVD. To
address these issues, functional analysis using genetic
engineering technologies, such as genome editing,should be used to examine the effects of SNPs or genes
on phenotypes.

Conclusions
Our finger cold-water immersion experiment revealed
that ~ 10% of the 94 Japanese study subjects did not show
an apparent CIVD response. Furthermore, the GWASsof CIVD identified 10 candidate genetic variants associated with CIVD-related traits in the subjects, and two
candidate functional genes, PRLR and COL4A2, werelikely involved in vascular function. Our findings suggestthat the focal SNPs may affect the blunted SkBF response
associated with eNO-independent and neurogenic activities, although SNP with a large effect on CIVD has yet to
be detected. It will be worthwhile investigating whether
PRLR and COL4A2 are related to EDHF-mediated vasodilation in peripheral tissues. Our results suggest thatgenetic polymorphisms are attributable to the interethnic variability of the CIVD response through differences in the allelic frequencies of SNPs associated with
CIVD. Using the largest-scale dataset of wavelet analysis
produced to date, we are the first to report GWASs for
CIVD; thus, this study improves our understanding of
the mechanisms underlying CIVD.

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Abbreviations
ANOVA	Analysis of variance
AVA	Arteriovenous anastomoseBM	Basement membrane
CADD	Combined Annotation Dependent Depletion
CDX	Chinese Dai in Xishuangbanna, China
CHB	Han Chinese in Beijing, China
CHS	Southern Han Chinese
CIVD	Cold-induced vasodilation
COL4A1	Collagen type IV alpha 1 chain
COL4A2	Collagen type IV alpha 2 chaineNO	Endothelial nitric oxide
EDHF	Endothelium-derived hyperpolarizing factor
GEE	Generalized estimating equation
GRASP	Genome-Wide Repository of Associations Between SNPs andPhenotypes
GWAS	Genome-wide association study
JPT	Japanese in Tokyo, Japan
KEGG	Kyoto Encyclopedia of Genes and Genomes
KHV	Kinh in Ho Chi Minh City, Vietnam
LD	Linkage disequilibrium
LINC01470	Long intergenic nonprotein coding RNA 1470
MAF	Minor allele frequency
NO	Nitric oxide
NOS	NO synthase
PCA	Principal component analysis
PRLR	Prolactin receptor
QC	Quality control
SkBF	Skin blood flow
SNP	Single-nucleotide polymorphism
VCF	Variant call format
1KGP	1000 Genomes Project

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40101-​023-​00319-2.
Additional file 1: Fig. S1. Wavelet amplitudes of neurogenic andendothelial NO-independent activities at phase 2 (prior to CIVD) in groupsof subjects differentiated by their CIVD response and genotypes of (a)rs73619297 of COL4A2 or (b) rs931740 of PRLR. Data are means ± standard
error of the mean. a.u., arbitrary units. Fig. S2. Principal component analysis of imputed genotype data from the JPQ cohort. Fig. S3. Wavelet amplitudes of endothelial NO-independent, endothelial NO-dependent, andneurogenic activities in four phases, including the baseline phase, in 10
study subjects with no CIVD response. Phases 1, 2, and 3 represent “vasoconstriction,” “prior to CIVD,” and “CIVD,” respectively [1]. a.u., arbitrary units.
Fig. S4. Physical position and linkage disequilibrium of candidate SNPsaround COL4A2 or at 14q23.3 in 94 Japanese subjects. The diagram wascreated using Haploview version 4.2. The numbers in diamonds represent
r2 values (×100). The haplotype block was defined using the method of
Gabriel et al. [2]. SNPs shown in bold are CIVD-associated SNPs. Table S1.
Characteristics of all study subjects including one female. Table S2. SNP
with p < 1×10-5 for CIVD-related categorical traits (CIVD response and
eNO-independent activity, neurogenic activity) in the logistic regression
analysis with the additive model. Table S3. SNP with p < 1×10-5 for CIVDrelated categorical traits (CIVD response and eNO-independent activity,
neurogenic activity) in the Fisher’s exact test with the additive model.
Table S4. Allele frequency of candidate SNPs detected in GWASs of CIVD
response with neurogenic and endothelial NO-independent activities.
Table S5. SNP with p < 1×10-5 for eNO-independent and neurogenic
activities in the additive model using two-way ANOVA. Table S6. SNP
with p < 1×10-5 for eNO-independent and neurogenic activities in thedominant model using two-way ANOVA. Table S7. Relationships among
the genotypes of five candidate SNPs defined in the dominant modeland network of three CIVD-related factors. Table S8. Relationships among
the genotypes of five candidate SNPs defined in the additive model andnetwork of three CIVD-related factors. Table S9. Correlation estimates offour candidate SNPs with their adjacent SNPs using LDproxy.

Page 14 of 16

Acknowledgements
The authors thank all the subjects for their participation in the study. Computations were partially performed on the NIG supercomputer at ROIS National
Institute of Genetics.
Authors’ contributions
T. S. and S. K. conceived the experiments. T. S., T. K., and S. K. designed theexperiments. T. S., T. K., Y. N., and M. U. collected and analyzed the data of theexperiments. T. S., T. K., and S. K. interpreted the data. Y. Y. contributed to thedesign of genetic analysis and interpretation of the data and drafted themanuscript. T. S. and S. K. revised the manuscript and provided critical feedback. All authors have read and approved the final manuscript.
Funding
This research was supported by JSPS KAKENHI Grant Number JP16H02529.
Availability of data and materials
Scientifically motivated request for data sharing will be considered. If justified,applicable parts of the data can be made available in anonymized format.

Declarations
Ethics approval and consent to participate
The study protocol complied with the Declaration of Helsinki and wasapproved by the local ethnics committee at Kyushu University (685–02).
Informed consent was obtained from all individual participants included in
the study.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Author details
1
Department of Genome Analysis, Institute of Biomedical Science, Kansai
Medical University, 2‑5‑1 Shin‑Machi, Hirakata, Osaka 573‑1010, Japan.
2
Department of Mechanical Engineering, Faculty of Engineering, Kyushu
University, 744 Motooka, Nishi‑Ku, Fukuoka 819‑0395, Japan. 3 Departmentof Mechanical Engineering, Graduate School of Engineering, Kyushu University, 744 Motooka, Nishi‑ku, Fukuoka 819‑0395, Japan. 4 Graduate Schoolof Systems Life Science, Kyushu University, 744 Motooka, Nishi‑ku, Fukuoka
819‑0395, Japan.
Received: 1 November 2022 Accepted: 10 February 2023

References
1. Lewis T. Observations upon the reactions of the vessels of the humanskin to cold. Heart. 1930;15:177–208.
2. Daanen HAM. Finger cold-induced vasodilation: a review. Eur J Appl
Physiol. 2003;89:411–26.
3. O’Brien C. Reproducibility of the cold-induced vasodilation response in
the human finger. J Appl Physiol. 2005;98:1334–40.
4. Cheung SS, Mekjavic IB. Cold-induced vasodilatation is not homogenous or generalizable across the hand and feet. Eur J Appl Physiol.
2007;99:701–5.
5. Kvandal P, Landsverk SA, Bernjak A, Stefanovska A, Kvernmo HD, Kirkebøen KA. Low-frequency oscillations of the laser Doppler perfusion signal
in human skin. Microvasc Res. 2006;72:120–7.
6. Bergersen TK, Hisdal J, Walløe L. Perfusion of the human finger during cold-induced vasodilatation. Am J Physiol Integr Comp Physiol.
1999;276:R731–7.
7. BintiMd Isa K, Kawasaki N, Ueyama K, Sumii T, Kudo S. Effects of cold
exposure and shear stress on endothelial nitric oxide synthase activation.
Biochem Biophys Res Commun. 2011;412:318–22.

Yasukochi et al. Journal of Physiological Anthropology

8.
9.
10.
11.
12.

13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.

28.
29.
30.
31.
32.
33.

(2023) 42:2

Hodges GJ, Mallette MM, Cheung SS. Cutaneous neural activity andendothelial involvement in cold-induced vasodilatation. Eur J Appl
Physiol. 2018;118:971–8.
Sera T, Kohno T, Nakashima Y, Uesugi M, Kudo S. Low-frequency oscillations of finger skin blood flow during the initial stage of cold-induced
vasodilation at different air temperatures. J Physiol Anthropol. 2020;39:37.
Hwa JJ, Ghibaudi L, Williams P, Chatterjee M. Comparison of acetylcholine-dependent relaxation in large and small arteries of rat mesentericvascular bed. Am J Physiol Circ Physiol. 1994;266:H952–8.
Shimokawa H. Reactive oxygen species in cardiovascular health anddisease: special references to nitric oxide, hydrogen peroxide, and Rhokinase. J Clin Biochem Nutr. 2020;66:83–91.
Kvandal P, Stefanovska A, Veber M, Kvermmo HD, Kirkebøen KA. Regulation of human cutaneous circulation evaluated by laser Doppler flowmetry, iontophoresis, and spectral analysis: importance of nitric oxide andprostaglandines. Microvasc Res. 2003;65:160–71.
Hodges GJ, Mallette MM, Tew GA, Saxton JM, Moss J, Ruddock AD, et al.
Effect of age on cutaneous vasomotor responses during local skin heating. Microvasc Res. 2017;112:47–52.
Kawai Y, Mimori T, Kojima K, Nariai N, Danjoh I, Saito R, et al. Japonica
array: Improved genotype imputation by designing a population-specific
SNP array with 1070 Japanese individuals. J Hum Genet. 2015;60:581–7.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al.
PLINK: a tool set for whole-genome association and population-based
linkage analyses. Am J Hum Genet. 2007;81:559–75.
O’Connell J, Gurdasani D, Delaneau O, Pirastu N, Ulivi S, Cocca M, et al.
A general approach for haplotype phasing across the full spectrum of
relatedness. PLoS Genet. 2014;10:e1004234.
Delaneau O, Marchini J, McVeanh GA, Donnelly P, Lunter G, Marchini JL,
et al. Integrating sequence and array data to create an improved 1000
Genomes Project haplotype reference panel. Nat Commun. 2014;5:3934.
Howie BN, Donnelly P, Marchini J, Rioux J, Xavier R, Taylor K, et al. A flexible
and accurate genotype imputation method for the next generation of
genome-wide association studies. PLoS Genet. 2009;5:e1000529.
The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature 2010;467:1061–73.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence
Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Bandelt HJ, Forster P, Rohl A. Median-joining networks for inferring
intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Clayton D. snpStats: SnpMatrix and XSnpMatrix classes and methods. R
package version 1.20.0. 2015.
R Core Team. R: a language and environment for statistical computing. https://​www.r-​proje​ct.​org/. 2016.
R Core Team. R: a language and environment for statistical computing. https://​www.r-​proje​ct.​org/. 2019.
RStudio Team. RStudio: integrated development environment for R.
http://​www.​rstud​io.​com/. RStudio, Inc. 2019.
Kato N, Takeuchi F, Tabara Y, Kelly TN, Go MJ, Sim X, et al. Meta-analysis of
genome-wide association studies identifies common variants associatedwith blood pressure variation in east Asians. Nat Genet. 2011;43:531–8.
Yasukochi Y, Sakuma J, Takeuchi I, Kato K, Oguri M, Fujimaki T, et al.
Longitudinal exome-wide association study to identify genetic susceptibility loci for hypertension in a Japanese population. Exp Mol Med.
2017;49:e409.
Matsukura M, Ozaki K, Takahashi A, Onouchi Y, Morizono T, Komai H,
et al. Genome-wide association study of peripheral arterial disease in a
Japanese population. Xu Z, editor. PLoS One. 2015;10:e0139262.
Klarin D, Lynch J, Aragam K, Chaffin M, Assimes TL, Huang J, et al.
Genome-wide association study of peripheral artery disease in the Million Veteran Program. Nat Med. 2019;25:1274–9.
Liang KY, Zeger SL. Longitudinal data analysis using generalized linear
models. Biometrika. 1986;73:13–22.
Hanley JA, Negassa A, Edwardes MD, Forrester JE. Statistical analysis ofcorrelated data using generalized estimating equations: an orientation.
Am J Epidemiol. 2003;157:364–75.
Halekoh U, Højsgaard S, Yan J. The R package geepack for generalized
estimating equations. J Stat Softw. 2006;15:1–11.
Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of
LD and haplotype maps. Bioinformatics. 2005;21:263–5.

Page 15 of 16

34. Machiela MJ, Chanock SJ. LDlink: a web-based application for exploring
population-specific haplotype structure and linking correlated alleles ofpossible functional variants. Bioinformatics. 2015;31:3555–7.
35. Leslie R, O’Donnell CJ, Johnson AD. GRASP: analysis of genotype-phenotype results from 1390 genome-wide association studies and corresponding open access database. Bioinformatics. 2014;30:185–94.
36. Buniello A, MacArthur JAL, Cerezo M, Harris LW, Hayhurst J, Malangone
C, et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic
Acids Res. 2019;47:D1005–12.
37. Rentzsch P, Witten D, Cooper GM, Shendure J, Kircher M. CADD:
predicting the deleteriousness of variants throughout the humangenome. Nucleic Acids Res. 2019;47:D886–94.
38. Chen CC, Stairs DB, Boxer RB, Belka GK, Horseman ND, Alvarez JV, et al.
Autocrine prolactin induced by the Pten-Akt pathway is required forlactation initiation and provides a direct link between the Akt andStat5 pathways. Genes Dev. 2012;26:2154–68.
39. Liu Y, Zhang J, Chen Y, Sohel H, Ke X, Chen J, et al. The correlation androle analysis of COL4A1 and COL4A2 in hepatocarcinogenesis. Aging.
2020;12:204–23.
40. Shiojima I, Walsh K. Role of Akt signaling in vascular homeostasis andangiogenesis. Circ Res. 2002;90:1243–50.
41. Kellogg DL, Liu Y, Kosiba IF, O’Donnell D. Role of nitric oxide in thevascular effects of local warming of the skin in humans. J Appl Physiol.
1999;86:1185–90.
42. Holowatz LA, Thompson-Torgerson CS, Kenney WL. Altered mechanisms of vasodilation in aged human skin. Exerc Sport Sci Rev.
2007;35:119–25.
43. Johnson JM, Kellogg DL. Local thermal control of the human cutaneous circulation. J Appl Physiol. 2010;109:1229–38.
44. Walløe L. Arterio-venous anastomoses in the human skin and their rolein temperature control. Temperature. 2016;3:92–103.
45. Patil MJ, Henry MA, Akopian AN. Prolactin receptor in regulation ofneuronal excitability and channels. Channels. 2014;8:193–202.
46. Corbacho AM, Macotela Y, Nava G, Eiserich JP, Cross CE, Martı́nez de laEscalera G, et al. Cytokine induction of prolactin receptors mediates
prolactin inhibition of nitric oxide synthesis in pulmonary fibroblasts.
FEBS Lett 2003;544:171–5.
47. Bachelot A, Binart N. Reproductive role of prolactin. Reproduction.
2007;133:361–9.
48. Radhakrishnan A, Raju R, Tuladhar N, Subbannayya T, Thomas JK, GoelR, et al. A pathway map of prolactin signaling. J Cell Commun Signal.
2012;6:169–73.
49. Zhao M, Gao F-H, Wang J-Y, Liu F, Yuan H-H, Zhang W-Y, et al. JAK2/
STAT3 signaling pathway activation mediates tumor angiogenesis byupregulation of VEGF and bFGF in non-small-cell lung cancer. Lung
Cancer. 2011;73:366–74.
50. Gao P, Niu N, Wei T, Tozawa H, Chen X, Zhang C, et al. The roles of signaltransducer and activator of transcription factor 3 in tumor angiogenesis. Oncotarget. 2017;8:69139–61.
51. Gonzalez C, Rosas-Hernandez H, Jurado-manzano B, Ramirez-Lee MA,
Salazar-Garcia S, Martinez-Cuevas PP, et al. The prolactin family hormones regulate vascular tone through NO and prostacyclin production
in isolated rat aortic rings. Acta Pharmacol Sin. 2015;36:572–86.
52. Pöschl E, Schlötzer-Schrehardt U, Brachvogel B, Saito K, Ninomiya Y,
Mayer U. Collagen IV is essential for basement membrane stability but
dispensable for initiation of its assembly during early development.
Development. 2004;131:1619–28.
53. Kuo DS, Labelle-Dumais C, Gould DB. COL4A1 and COL4A2 mutations
and disease: insights into pathogenic mechanisms and potential therapeutic targets. Hum Mol Genet. 2012;21:R97-110.
54. Marchand M, Monnot C, Muller L, Germain S. Extracellular matrix scaffolding in angiogenesis and capillary homeostasis. Semin Cell Dev Biol.
2019;89:147–56.
55. Verbeek E, Meuwissen MEC, Verheijen FW, Govaert PP, Licht DJ, Kuo
DS, et al. COL4A2 mutation associated with familial porencephaly andsmall-vessel disease. Eur J Hum Genet. 2012;20:844–51.
56. Dimmeler S, Fleming I, Fisslthaler B, Hermann C, Busse R, Zeiher AM.
Activation of nitric oxide synthase in endothelial cells by Akt-dependent phosphorylation. Nature. 1999;399:601–5.

Yasukochi et al. Journal of Physiological Anthropology(2023) 42:2

Page 16 of 16

57. Quillon A, Fromy B, Debret R. Endothelium microenvironment sensing
leading to nitric oxide mediated vasodilation: a review of nervous andbiomechanical signals. Nitric Oxide. 2015;45:20–6.
58. Shimokawa H, Morikawa K. Hydrogen peroxide is an endotheliumderived hyperpolarizing factor in animals and humans. J Mol Cell
Cardiol. 2005;39:725–32.
59. Takaki A, Morikawa K, Murayama Y, Yamagishi H, Hosoya M, Ohashi J, et al.
Roles of endothelial oxidases in endothelium-derived hyperpolarizing
factor responses in mice. J Cardiovasc Pharmacol. 2008;52:510–7.
60. Takaki A, Morikawa K, Tsutsui M, Murayama Y, Tekes E, Yamagishi H, et al.
Crucial role of nitric oxide synthases system in endothelium-dependent
hyperpolarization in mice. J Exp Med. 2008;205:2053–63.
61. Ndiaye M, Chataigneau T, Chataigneau M, Schini-Kerth VB. Red wine polyphenols induce EDHF-mediated relaxations in porcine coronary arteries
through the redox-sensitive activation of the PI3-kinase/Akt pathway. Br J
Pharmacol. 2004;142:1131–6.
62. Schini-Kerth VB, Auger C, Kim JH, Étienne-Selloum N, Chataigneau T.
Nutritional improvement of the endothelial control of vascular tone bypolyphenols: role of NO and EDHF. Pflügers Arch. 2010;459:853–62.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research ? Choose BMC and benefit from:• fast, convenient online submission• thorough peer review by experienced researchers in your field• rapid publication on acceptance• support for research data, including large and complex data types• gold Open Access which fosters wider collaboration and increased citations• maximum visibility for your research: over 100M website views per year
At BMC, research is always in progress.
Learn more biomedcentral.com/submissions