Skip to main content
  • Original article
  • Open access
  • Published:

Upregulation of cathepsin L gene under mild cold conditions in young Japanese male adults



Physiological thermoregulatory systems in humans have been a key factor for adaptation to local environments after their exodus from Africa, particularly, to cold environments outside Africa. Recent studies using high-throughput sequencing have identified various genes responsible for cold adaptation. However, the molecular mechanisms underlying initial thermoregulation in response to acute cold exposure remain unclear. Therefore, we investigated transcriptional profiles of six young Japanese male adults exposed to acute cold stress.


In a climatic chamber, the air temperature was maintained at 28°C for 65 min and was then gradually decreased to 19°C for 70 min. Saliva samples were obtained from the subjects at 28°C before and after 19°C cold exposure and were used for RNA sequencing.


In the cold exposure experiment, expression levels of 14 genes were significantly changed [false discovery rate (FDR) < 0.05] although the degree of transcriptional changes was not high due to experimental conditions or blunted transcriptional reaction in saliva to cold stress. As a result, differential gene expression analyses detected the cathepsin L (CTSL) gene to be significantly upregulated, with FDR < 0.05 and log2 fold change value > 1; thus, this gene was identified as a differentially expressed gene. Given that the cathepsin L protein is related to invasion of the novel coronavirus (SARS-CoV-2), mild cold stress might alter the susceptibility to coronavirus disease-19 in humans. The gene ontology enrichment analysis for 14 genes with FDR < 0.05 suggested that immune-related molecules could be activated by mild cold stress.


The results obtained from this study indicate that CTSL expression levels can be altered by acute mild cold stress.


Anatomically modern humans emerged in Africa ~300,000 to 200,000 years ago [1]; they then dispersed worldwide from Africa many tens of thousands of years ago, but the timing remains controversial [2,3,4,5,6]. This migration event is called “Out-of-Africa” and these migrating human populations have genetically adapted to various local environments [7, 8]. Physiological thermoregulatory ability is one of the crucial factors for adapting to local environments, especially to cold environments outside Africa. In response to cold conditions, skin blood vessels are constricted to minimize heat loss from the human body, while metabolic heat production can compensate the heat loss occurring after vasoconstriction [9]. Humans have two main types of adipose tissue, white (WAT) and brown (BAT) adipose tissues, along with a third adipose tissue called the beige or brite adipose tissue that is present in WAT (brown-like thermogenic adipocytes). BAT is involved in thermoregulatory nonshivering thermogenesis (NST) through adrenergic stimulation [10,11,12,13,14]. NST is defined as “heat production due to metabolic energy transformation by processes that do not involve contraction of the skeletal muscles,” whereas shivering thermogenesis is defined as “an increase in the rate of heat production during cold exposure due to increased contractile activity of the skeletal muscles not involving voluntary movements and external work” [15]. Thermoregulatory ability can be affected by lifestyles such as exercise and seasonal acclimatization [9, 16] whereas repeated cold exposure can induce insulative cold adaptation without altered metabolic response [17].

Mutations in NST-related genes may confer resistance to cold environments. Recent studies have revealed genetic factors associated with the NST response to cold stress. The uncoupling protein 1 (UCP1) gene is highly expressed in BAT related to adaptive adrenergic NST [18], resulting from the uncoupling of mitochondrial respiration from ATP synthesis. Single-nucleotide polymorphisms (SNPs) or haplotypes of UCP1 can affect adaptive responses to cold environments as well as the development of obesity through energy expenditure in human populations [19, 20] and even in Japanese individuals [21,22,23]. Further, Japanese individuals with the AA genotype of SNP rs1057001 in the tribbles pseudokinase 2 (TRIB2) gene showed relatively high expression levels of several thermogenic-related genes [24], suggesting that this genotype is adaptive to cold climates. A previous study has reported that positive natural selection may have operated on a regulatory genetic variant near the transient receptor potential cation channel subfamily M member 8 (TRPM8) gene in Eurasian populations for cold adaptation [25], whereas strong signals of positive selection were detected in genomic regions containing carnitine palmitoyltransferase 1A (CPT1A) and LDL receptor-related protein 5 (LRP5) genes on chromosome 11 in Northeastern Siberian populations [26].

The epigenetic modification also contributes to beige/brite adipocyte thermogenesis through the TET1 (tet methylcytosine dioxygenase 1) demethylase suppression by cold stimuli [27]. RNA sequencing (RNA-seq) and lipidomic analyses of BAT in mice exposed to cold conditions exhibited upregulation of several brown adipocyte genes related to thermogenesis; further, expression levels of genes involved in glycerolipid synthesis and fatty acid elongation were also remarkably enhanced by cold exposure [28]. RNA-seq analyses of humans showed that high UCP1 expression in the epicardial adipose tissue possibly downregulates gene expression related to immune responses such as T cell-related processes [29]. Further, transcriptome analysis using a microarray-based assay to investigate shared molecular adaptive responses between cold acclimation (14–15°C) and exercise training in the skeletal muscles of patients with type 2 diabetes showed overlapping effects on the transcriptional profiles involved in tissue remodeling [30]. A single-nucleus RNA-seq analysis of mice and human adipose tissues revealed that adipocyte subpopulations, which were identified by “P4 cells” in the study, likely regulate the thermogenic function of other adipocytes [31]. Although many studies on cold adaptation have been reported, little is known about the effects of acute cold stress on the transcriptional profiles of humans in vivo.

As mentioned above, previous genetic studies on cold adaptation have mainly focused on genetic factors responsible for the NST in BAT. However, vasoconstriction of the skin vessels plays a key role in the initial physiological response against acute cold exposure, and this response results in suppression of heat loss from the body to cold environments. To the best of our knowledge, the molecular mechanisms involved in initial thermoregulation remain unknown. Therefore, in this study, we performed the transcriptomic analysis of young Japanese male adults who were exposed to acute mild cold stress, aiming to elucidate the molecular mechanisms underlying physiological responses to cold environments.


Study subjects

In total, six young Japanese males aged > 20 years (mean age ± standard deviation, 21.8 ± 1.2 years) were recruited from the subject group that participated in our previous study on thermoregulatory and circulatory responses in hypobaric hypoxia and cold [32]. The anthropometric data of the study subjects are shown in Table 1. The subjects had abstained from drinking alcohol, doing exercise, and smoking for 24 h before the experiment, and had been prohibited from drinking caffeine and eating 2 h prior to the experiment, which was conducted in a climatic chamber of the Research Center for Environmental Adaptation (Fukuoka, Japan).

Table 1 Physical characteristics of six Japanese male subjects

Experimental design

The cold exposure experiment was carried out from October to December 2018. The study subjects who wore undershorts and short-sleeve T-shirts (~0.13 clo) participated in the experiment. The total duration of the experiment was 135 min, and the subjects rested in a supine position during the experiment. The air temperature (Tair) was maintained at 28°C for 65 min, gradually decreased to 19°C for 70 min, and then returned to 28°C for 30 min. Humidity was maintained at 50% RH (relative humidity) during the experiment. For transcriptome analyses, we collected the subjects’ saliva with a volume ≥ 2 mL using the Oragene® RNA kit (DNA Genotek, Ottawa, Canada) at 28°C before and after 19°C cold exposure.

RNA sequencing and differential gene expression analyses

We extracted total RNA from saliva samples, which contain buccal cells and blood leukocytes, using either the RNeasy micro kit (Qiagen, Hilden, Germany) or the NucleoSpin® RNA kit (Takara Bio, Shiga, Japan) according to the respective manufacturers’ protocols. The RNA quality for each sample is shown in Additional file 1: Table S1. The sequence library was prepared using the TruSeq® Stranded Total RNA Library Prep kit (Illumina, San Diego, CA, USA) after RNA integrity was determined using an Agilent 2100 Bioanalyzer or TapeStation (Agilent Technologies, Santa Clara, CA, USA). RNA-seq was conducted on an Illumina NovaSeq 6000 platform (Illumina) with 100-bp paired-end (PE) reads. The quality check, library preparation, and RNA-seq were performed by Macrogen (Seoul, South Korea). We removed the Illumina PE adapters and low-quality reads using fastp ver. 0.19.7 [33]. PE reads were mapped to the human reference genome (GRCh38) using STAR ver. 2.5.3 [34]. We estimated the gene expression levels from the mapped BAM (binary alignment map) files using RSEM ver. 1.3.0 [35].

Based on the read-count data at the gene level generated using the RSEM, genes with < 13 read counts across all samples were filtered out for further analysis. Normalization of count data and estimation of log2 fold change (LFC) between the experimental conditions were performed using the Bioconductor R package DESeq2 ver. 1.22.2 [36] via R software ver. 3.5.3 [37] in RStudio ver. 1.2.5019 [38]. The p value for the significance of gene expression differences was adjusted to the false discovery rate (FDR) using the Benjamini and Hochberg method [39] for multiple testing. The differentially expressed gene (DEG) with an FDR < 0.05 and |LFC| ≥ 1 was identified. We conducted principal component analysis (PCA) to compare transcriptomic similarity among samples. The detailed methods mentioned above have been described previously [40]. The molecular function of DEG was predicted using the information on the AmiGO 2 (, [41]) and Kyoto Encyclopedia of Genes and Genomes (KEGG; databases. Gene Ontology (GO) enrichment analysis was carried out using 11 upregulated and 3 downregulated genes with FDR < 0.05 in the Bioconductor R package topGO ver. 2.34.0 [42, 43].

Real-time reverse transcription quantitative PCR analysis

To confirm the expression levels of DEG, we performed real-time reverse transcription quantitative PCR (real-time RT-qPCR) analysis. The total RNA was isolated from the saliva samples mentioned above, according to combined protocols of the Oragene® RNA kit and mirVana™ miRNA isolation kit (Ambion, Austin, TX, USA) without adding the Lysis/Binding solution and miRNA Homogenate Additive [44, 45]. The quality of total RNA was evaluated using a Nanodrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) (Additional file 1: Table S1). Two saliva samples from an individual were removed for further analyses due to low sample volume.

Complementary DNA (cDNA) was synthesized from total RNA by reverse transcription with oligo(dT)20 primers, which are useful for excluding bacterial RNA in saliva samples [45], using the Superscript III First-Strand Synthesis System (Thermo Fisher Scientific). The experimental procedure was conducted according to the manufacturer’s protocol. Real-time RT-qPCR with two technical replicates for each of five biological replicates was performed using Applied Biosystems™ TaqMan® Fast Advanced Master Mix (Thermo Fisher Scientific) and gene-specific primer probes for the cathepsin L (CTSL, Hs00964650_m1) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH, Hs03929097_g1) genes, according to the manufacturer’s protocols. For real-time RT-qPCR, Thermal Cycler Dice® Real-Time System Lite TP760 (Takara Bio) was used with the following cycling conditions; 2 min at 50 °C, 20 s at 95 °C, and 50 cycles of 3 s at 95 °C and 30 s at 60 °C. Relative gene expression levels were determined using the comparative Ct (2-ΔΔCt) method [46] with the second derivative maximum method. The expression levels of the target gene were normalized to those of GAPDH.

Results and discussion

RNA-seq data of 12 saliva samples

We obtained twelve saliva samples from six study subjects to examine the transcriptomic changes after mild cold exposure. RNA-seq analysis for saliva samples generated a total of ~562 million reads (46.8 ± 5.4 million reads). After quality control using the fastp, 20.9 ± 2.8 million reads (16.3–27.3 million reads) were uniquely mapped to the human reference genome assembly. PCA (Fig. 1a) and Euclidean distance analysis (Fig. 1b) indicated that the expression vectors of samples from the same individuals tended to be clustered, irrespective of experimental conditions. This result is similar to that of our previous study on transcriptomic changes in response to hypobaric hypoxic exposure [40]. The transcriptomic profiles of 50 genes with the highest variance across samples also showed the clustering of vectors from each individual (Fig. 1c). These results indicate that transcriptional fluctuations in response to mild cold exposure at Tair of 19°C were limited within each individual.

Fig. 1
figure 1

Transcriptomic similarities among 12 samples from six Japanese participants in the cold exposure experiments. a Principal component analysis of RNA-seq data plotted according to the first (horizontal axis) and second (vertical axis) principal components. b Heatmap of Euclidean distance between expression vectors. Euclidean distance was calculated using the R dist function. The color intensity reflects the transcriptomic similarity between samples. c Hierarchical clustering of expression for 50 genes with the highest variance across samples. Colors represent normalized gene expression levels

The limited expression levels may be attributed to the following: (1) the intensity or duration of cold stress may have been insufficient to induce drastic transcriptional changes, (2) saliva samples containing buccal cells and leukocytes were less sensitive to cold exposure, (3) the timing of sample collection possibly decreases the statistical power to detect DEGs because we investigated transcriptional profiles after cold exposure but not during the exposure, and (4) seasonal acclimatization related to the low outside temperature in winter (October to December) may have affected the sensitivity of transcriptional responses to cold exposure in this experiment. Transcriptional profiles differ between WAT and BAT in humans [47] and mice [48] during cold exposure; thus, it is possible that the transcriptional profiles of saliva samples are largely different from those of adipose tissues under cold conditions and that the transcriptional dynamics are less sensitive to cold stress. To understand the molecular mechanisms underlying physiological responses to cold exposure, further analyses for transcriptomic profiles at several time points during different intensities of cold stress across various cell lines or tissues are required.

A differentially expressed gene in response to cold exposure

To explore the significant changes in gene expression levels before and after 19°C cold stress, we conducted differential expression analysis. This analysis detected 11 and 3 significantly (FDR < 0.05) upregulated and downregulated genes, respectively (Table 2). However, the |LFC| value of 13 genes was less than unity (equivalently, their expression levels after cold exposure were less than twice compared to those before the exposure), whereas the value of the cathepsin L (CTSL) gene alone was higher than unity (Fig. 2; LFC = 1.03, FDR = 0.049). Therefore, we identified CTSL as the DEG in the present study. Although the remaining genes did not reach the |LFC| threshold, some of the genes with smaller changes had sufficient read counts compared with CTSL, and these changes seemed to be sufficiently reliable. Therefore, we identified the 13 genes as moderate DEGs. Differential expression analysis for the hypobaric hypoxic experiment in our previous study identified 30 upregulated DEGs [40], whereas the present study on cold exposure detected a single upregulated DEG (i.e., CTSL). The sample size of the previous study (three individuals) was smaller than that of the present study (six individuals); additionally, we confirmed no significant difference in the gene expression levels under no-stimulus conditions [40]. Consequently, the small number of the DEGs observed in the present study is unlikely due to the samples size and might be attributed to the possibilities described above (i.e., mild cold stress, short-term cold exposure, blunted transcriptional response in saliva to cold stress, or seasonal acclimatization).

Table 2 Fourteen protein coding genes with FDR < 0.05 in the differential gene expression analysis
Fig. 2
figure 2

Volcano plot for differential expression analysis between two experimental conditions, 28°C before and after 19°C cold exposure. The vertical and horizontal axes represent significant difference as –log10(FDR) and LFC, respectively. Broken lines represent the significance or LFC threshold

Next, we compared the expression levels of CTSL in each individual before and after cold exposure (Fig. 3). The comparative analysis indicated that the magnitude of gene expression changes after cold exposure varied with individuals. It would be intriguing to survey the genetic variants that impact the expression levels of CTSL in Japanese individuals. Promoter or enhancer regions are more likely to be the factor affecting transcriptional variations; however, we did not collect DNA samples. It will thus be desirable to examine genetic variants in the CTSL genomic region, including neighboring regions, to detect the genetic factors related to the differences in the transcriptional changes. The expression levels of the other 13 genes with FDR < 0.05 are shown in Additional file 2: Fig. S1. As with the CTSL, the interindividual variability of transcriptome profiles was observed across the genes.

Fig. 3
figure 3

Gene expression levels of CTSL with LFC = 1.03, identified using the DESeq2 program, in each of the study subjects. Normalized counts at the vertical axis represent the estimated expression abundance at the gene level

To confirm the difference in the expression levels of CTSL between the two experimental conditions (before and after 19°C cold exposure), a relative quantification analysis was conducted using the 2-ΔΔCt method after real-time RT-qPCR using saliva samples from five individuals. The mean values of relative quantification in the study subjects were higher after cold exposure than before the exposure with a mean fold change of 2.87 (Additional file 3: Fig. S2), suggesting that the expression levels of CTSL were upregulated after acute cold exposure. However, the difference in the mean values of relative quantification was not statistically significant (p = 0.44, by Wilcoxon signed-rank test). It is possible that the salivary RNA used for the RT-qPCR was degraded because this experiment was conducted at least 1 year later after the RNA-seq analysis. In fact, the gene expression level of subject 4 before cold exposure appeared to be curiously higher than those of other subjects (Additional file 3: Fig S2). This may result from the low quality of RNA samples, even though the transcriptional data were normalized. The results of RT-qPCR analysis showed similar trends to those of RNA-seq analysis; however, further analysis of real-time RT-qPCR with a larger sample size is required in the future.

Gene ontology enrichment analysis

To estimate which biological process is activated by cold stimuli, we performed topGO enrichment analysis using the DEG and moderate DEGs. The hierarchical tree of GO terms is shown in Additional file 3: Fig. S3. The enrichment analysis identified a total of eight significantly overrepresented GO terms in the biological process category (FDR < 0.05, by Fisher’s exact test). All the GO terms were related to migrations of white blood cell subtypes such as leukocyte migration (Table 3); thus, molecules involved in immune functions may be activated by the cold exposure experiments. In fact, it has been reported that acute mild cold stress may affect the proliferation of several white blood cell subtypes although effects of the proliferation on host defense remain unknown [49]. Incidentally, the expression change in the CD177 molecule (CD177) gene was the most statistically significant (FDR = 6.46×10-5) in the differential expression analysis. CD177 can be a neutrophil activation marker, and it has been reported that the expression levels were positively correlated with the severity of coronavirus disease-19 (COVID-19) [50]. Since the LFC value of CD177 (LFC = 0.999) was quite close to the threshold (LFC = 1.000), it is possible that CD177 is upregulated by cold stress. We will examine the transcriptional alteration of this gene under improved experimental conditions in the future.

Table 3 GO processes estimated by topGO enrichment analysis using 11 upregulated and 3 downregulated genes with FDR < 0.05

Functional prediction of the cathepsin L gene

The cathepsin L protein is a lysosomal cysteine proteinase involved in the catabolism of intracellular proteins such as collagen and elastin [51, 52]. The GO terms of the human CTSL included fibronectin binding, Toll-like receptor signaling pathway, adaptive immune response, cysteine-type endopeptidase activity, and protein binding. According to the KEGG database, CTSL plays major roles in immune responses including antigen processing and presentation (Table 4). The cathepsin L protein has various functional activities including degradation of cellular proteins and autophagic degradation [53].

Table 4 Kyoto encyclopedia of genes and genomes (KEGG) pathway of CTSL

Notably, the cathepsin L protein is reported to be involved in cleaving the spike protein of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [54,55,56]. The inhibition of cathepsin L can decrease viral infectivity, such as herpes simplex virus 2 [57]. In addition, the cathepsin L protein is required for SARS-CoV-2 invasion into the host cell [54, 56] and cold conditions may be a potential risk factor in the spread of COVID-19 [58]. One may hypothesize that cold stress increases susceptibility to COVID-19 infection in humans via the transcriptional changes in CTSL. Further analyses for transcriptional changes in response to different intensities of cold stress with different sampling points are required to examine whether consistent CTSL expression patterns are observed.

The present study had several limitations. First, the transcriptional profiles observed in the present study were derived from saliva samples only; thus, transcriptome analyses of other tissues such as skin and blood are required because the gene expression patterns among tissues are different. Second, the effects of middle- or long-term acclimatization to cold stress remain unclear. Third, as the significant gene expression change in CTSL was not replicated, a replication study in other Japanese cohorts is required to validate this expression change. Fourth, the functional relevance of CTSL upregulation to physiological changes in response to cold stress remains unclear and needs to be elucidated.


In conclusion, relatively mild cold stress did not greatly alter the transcriptional profiles of saliva samples from the study subjects, but the expression level of CTSL was significantly upregulated. This upregulation might be related to the immune responses induced by short-term cold stress. |LFC| values of 13 moderate DEGs were less than unity; however, transcriptional changes in some of the moderate DEGs might be involved in the physiological response to cold stress because some of these moderate DEGs showed a sufficient significant difference in transcriptional profiles compared with CTSL.

Availability of data and materials

The datasets presented in this article are not readily available because of ethical and privacy issues. However, scientifically motivated requests for data sharing will be considered by the reviewing ethical committee of our institute.



Binary alignment map


Brown adipose tissue


Complementary DNA


Coronavirus disease-19


Carnitine palmitoyltransferase 1A


Cathepsin L


Differentially expressed gene


False discovery rate


Glyceraldehyde-3-phosphate dehydrogenase


Gene ontology


Kyoto Encyclopedia of Genes and Genomes


Log2 fold change


LDL receptor-related protein 5


Nonshivering thermogenesis


Principal component analysis




Relative humidity


RNA sequencing


Reverse transcription quantitative PCR


Severe acute respiratory syndrome coronavirus 2


Single-nucleotide polymorphism

Tair :

Air temperature


Tet methylcytosine dioxygenase 1


Tribbles pseudokinase 2


Transient receptor potential cation channel subfamily M member 8


Uncoupling protein 1


White adipose tissue


  1. Hublin J-J, Ben-Ncer A, Bailey SE, Freidline SE, Neubauer S, Skinner MM, et al. New fossils from Jebel Irhoud, Morocco and the pan-African origin of Homo sapiens. Nature. 2017;546:289–92.

    Article  CAS  PubMed  Google Scholar 

  2. Wall JD. Inferring human demographic histories of non-African populations from patterns of allele sharing. Am J Hum Genet. 2017;100:766–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Reyes-Centeno H, Hubbe M, Hanihara T, Stringer C, Harvati K. Testing modern human out-of-Africa dispersal models and implications for modern human origins. J Hum Evol. 2015;87:95–106.

    Article  PubMed  Google Scholar 

  4. Groucutt HS, Petraglia MD, Bailey G, Scerri EML, Parton A, Clark-Balzan L, et al. Rethinking the dispersal of Homo sapiens out of Africa. Evol Anthropol. 2015;24:149–64.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Valverde G, Zhou H, Lippold S, De Filippo C, Tang K, Herráez DL, et al. A novel candidate region for genetic adaptation to high altitude in Andean populations. PLoS One. 2015;10:1–22.

    Article  Google Scholar 

  6. Hershkovitz I, Weber GW, Quam R, Duval M, Grün R, Kinsley L, et al. The earliest modern humans outside Africa. Science. 2018;359:456–9.

    Article  CAS  PubMed  Google Scholar 

  7. Fan S, Hansen MEB, Lo Y, Tishkoff SA. Going global by adapting local: a review of recent human adaptation. Science. 2016;354:54–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Ilardo M, Nielsen R. Human adaptation to extreme environmental conditions. Curr Opin Genet Dev. 2018;53:77–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Maeda T. Relationship between maximum oxygen uptake and peripheral vasoconstriction in a cold environment. J Physiol Anthropol. 2017;36:42.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Enerbäck S, Jacobsson A, Simpson EM, Guerra C, Yamashita H, Harper ME, et al. Mice lacking mitochondrial uncoupling protein are cold-sensitive but not obese. Nature. 1997;387:90–4.

    Article  PubMed  Google Scholar 

  11. Ravussin E, Galgani JE. The implication of brown adipose tissue for humans. Annu Rev Nutr. 2011;31:33–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Anunciado-Koza R, Ukropec J, Koza RA, Kozak LP. Inactivation of UCP1 and the glycerol phosphate cycle synergistically increases energy expenditure to resist diet-induced obesity. J Biol Chem. 2008;283:27688–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Cannon B, Nedergaard J. Nonshivering thermogenesis and its adequate measurement in metabolic studies. J Exp Biol. 2011;214:242–53.

    Article  PubMed  Google Scholar 

  14. Cannon B, Nedergaard J. Brown adipose tissue: function and physiological significance. Physiol Rev. 2004;84:277–359.

    Article  CAS  PubMed  Google Scholar 

  15. IUPS Thermal Commission. Glossary of terms for thermal physiology. Pflügers Arch. 1987;410:567–87.

    Article  Google Scholar 

  16. Nishimura T, Motoi M, Egashira Y, Choi D, Aoyagi K, Watanuki S. Seasonal variation of non-shivering thermogenesis (NST) during mild cold exposure. J Physiol Anthropol. 2015;34:1–6.

    Article  Google Scholar 

  17. Wakabayashi H, Wijayanto T, Kuroki H, Lee JY, Tochihara Y. The effect of repeated mild cold water immersions on the adaptation of the vasomotor responses. Int J Biometeorol. 2012;56:631–7.

    Article  PubMed  Google Scholar 

  18. Golozoubova V, Cannon B, Nedergaard J. UCP1 is essential for adaptive adrenergic nonshivering thermogenesis. Am J Physiol Metab. 2006;291:E350–7.

    CAS  Google Scholar 

  19. Lappalainen T, Salmela E, Andersen PM, Dahlman-Wright K, Sistonen P, Savontaus ML, et al. Genomic landscape of positive natural selection in Northern European populations. Eur J Hum Genet. 2010;18:471–8.

    Article  PubMed  Google Scholar 

  20. Hancock AM, Clark VJ, Qian Y, Di Rienzo A. Population genetic analysis of the uncoupling proteins supports a role for UCP3 in human cold resistance. Mol Biol Evol. 2011;28:601–14.

    Article  CAS  PubMed  Google Scholar 

  21. Yoneshiro T, Ogawa T, Okamoto N, Matsushita M, Aita S, Kameya T, et al. Impact of UCP1 and β3AR gene polymorphisms on age-related changes in brown adipose tissue and adiposity in humans. Int J Obes. 2013;37:993–8.

    Article  CAS  Google Scholar 

  22. Nishimura T, Katsumura T, Motoi M, Oota H, Watanuki S. Experimental evidence reveals the UCP1 genotype changes the oxygen consumption attributed to non-shivering thermogenesis in humans. Sci Rep. 2017;7:5570.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Matsushita H, Kurabayashi T, Tomita M, Kato N, Tanaka K. Effects of uncoupling protein 1 and ß3-adrenergic receptor gene polymorphisms on body size and serum lipid concentrations in Japanese women. Maturitas. 2003;45:39–45.

    Article  CAS  PubMed  Google Scholar 

  24. Nakayama K, Iwamoto S. An adaptive variant of TRIB2, rs1057001, is associated with higher expression levels of thermogenic genes in human subcutaneous and visceral adipose tissues. J Physiol Anthropol. 2017;36:16.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Key FM, Abdul-Aziz MA, Mundry R, Peter BM, Sekar A, D’Amato M, et al. Human local adaptation of the TRPM8 cold receptor along a latitudinal cline. PLoS Genet. 2018;14:1–22.

    Article  CAS  Google Scholar 

  26. Cardona A, Pagani L, Antao T, Lawson DJ, Eichstaedt CA, Yngvadottir B, et al. Genome-wide analysis of cold adaptation in indigenous siberian populations. Kayser M, editor. PLoS One. 2014;9:e98076.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Damal Villivalam S, You D, Kim J, Lim HW, Xiao H, Zushin PJH, et al. TET1 is a beige adipocyte-selective epigenetic suppressor of thermogenesis. Nat Commun. 2020;11:1–14.

    Article  CAS  Google Scholar 

  28. Marcher AB, Loft A, Nielsen R, Vihervaara T, Madsen JGS, Sysi-Aho M, et al. RNA-seq and mass-spectrometry-based lipidomics reveal extensive changes of glycerolipid pathways in brown adipose tissue in response to cold. Cell Rep. 2015;13:2000–13.

    Article  CAS  PubMed  Google Scholar 

  29. Chechi K, Vijay J, Voisine P, Mathieu P, Bossé Y, Tchernof A, et al. UCP1 expression–associated gene signatures of human epicardial adipose tissue. JCI Insight. 2019;4:e123618.

    Article  PubMed Central  Google Scholar 

  30. Nascimento EBM, Hangelbroek RWJ, Hooiveld GJEJ, Hoeks J, Van Marken Lichtenbelt WD, Hesselink MHC, et al. Comparative transcriptome analysis of human skeletal muscle in response to cold acclimation and exercise training in human volunteers. BMC Med Genomics. 2020;13:124.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Sun W, Dong H, Balaz M, Slyper M, Drokhlyansky E, Colleluori G, et al. snRNA-seq reveals a subpopulation of adipocytes that regulates thermogenesis. Nature. 2020;587:98–102.

    Article  CAS  PubMed  Google Scholar 

  32. Shin S, Yasukochi Y, Wakabayashi H, Maeda T. Effects of acute hypobaric hypoxia on thermoregulatory and circulatory responses during cold air exposure. J Physiol Anthropol. 2020;39:28.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  PubMed  Google Scholar 

  35. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. R Core Team. R: a language and environment for statistical computing. 2019.

    Google Scholar 

  38. RStudio Team. RStudio: integrated development environment for R. RStudio, Inc; 2019.

    Google Scholar 

  39. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.

    Google Scholar 

  40. Yasukochi Y, Shin S, Wakabayashi H, Maeda T. Transcriptomic changes in young Japanese males after exposure to acute hypobaric hypoxia. Front Genet. 2020;11:559074.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22:1600–7.

    Article  CAS  PubMed  Google Scholar 

  43. Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for Gene Ontology. R package version 2.34.0; 2018.

    Google Scholar 

  44. Patel RS, Jakymiw A, Yao B, Pauley BA, Carcamo WC, Katz J, et al. High resolution of microRNA signatures in human whole saliva. Arch Oral Biol. 2011;56:1506–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Ostheim P, Tichý A, Sirak I, Davidkova M, Stastna MM, Kultova G, et al. Overcoming challenges in human saliva gene expression measurements. Sci Rep. 2020;10:11147.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  47. Chondronikola M, Volpi E, Børsheim E, Porter C, Saraf MK, Annamalai P, et al. Brown adipose tissue activation is linked to distinct systemic effects on lipid metabolism in humans. Cell Metab. 2016;23:1200–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Hao Q, Yadav R, Basse AL, Petersen S, Sonne SB, Rasmussen S, et al. Transcriptome profiling of brown adipose tissue during cold exposure reveals extensive regulation of glucose metabolism. Am J Physiol Metab. 2015;308:E380–92.

    CAS  Google Scholar 

  49. Castellani JW, Brenner IK, Rhind SG. Cold exposure: human immune responses and intracellular cytokine expression. Med Sci Sport Exerc. 2002;34:2013–20.

    Article  CAS  Google Scholar 

  50. Lévy Y, Wiedemann A, Hejblum BP, Durand M, Lefebvre C, Surénaud M, et al. CD177, a specific marker of neutrophil activation, is associated with coronavirus disease 2019 severity and death. iScience. 2021;24:102711.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Du X, Chen NLH, Wong A, Craik CS, Brömme D. Elastin degradation by cathepsin v requires two exosites. J Biol Chem. 2013;288:34871–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Uchiyama Y, Waguri S, Sato N, Watanabe T, Ishido K, Kominami E. Cell and tissue distribution of lysosomal cysteine proteinases, cathepsins B, H, and L, and their biological roles. Acta Histochem Cytochem. 1994;27:287–308.

    Article  CAS  Google Scholar 

  53. Dana D, Pathak SK. A review of small molecule inhibitors and functional probes of human cathepsin L. Molecules. 2020;25:698.

    Article  CAS  PubMed Central  Google Scholar 

  54. Yousefi B, Valizadeh S, Ghaffari H, Vahedi A, Karbalaei M, Eslami M. A global treatments for coronaviruses including COVID-19. J Cell Physiol. 2020;235:9133–42.

    Article  CAS  PubMed  Google Scholar 

  55. Zhao M-M, Yang W-L, Yang F-Y, Zhang L, Huang W-J, Hou W, et al. Cathepsin L plays a key role in SARS-CoV-2 infection in humans and humanized mice and is a promising target for new drug development. Signal Transduct Target Ther. 2021;6:134.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Gomes CP, Fernandes DE, Casimiro F, da Mata GF, Passos MT, Varela P, et al. Cathepsin L in COVID-19: from pharmacological evidences to genetics. Front Cell Infect Microbiol. 2020;10:777.

    Article  Google Scholar 

  57. Hopkins J, Yadavalli T, Agelidis AM, Shukla D. Host enzymes heparanase and cathepsin L promote herpes simplex virus 2 release from cells. Longnecker RM, editor. J Virol. 2018;92:e01179–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Mecenas P, Bastos RTDRM, Vallinoto ACR, Normando D. Effects of temperature and humidity on the spread of COVID-19: a systematic review. Samy AM, editor. PLoS One. 2020;15:e0238339.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


The authors deeply thank all the subjects for their participation in the study. We owe special thanks to Dr. Midori Motoi for sample management. Computations were partially performed on the NIG supercomputer at ROIS National Institute of Genetics.


This work was supported by the JSPS KAKENHI Grant-in-Aid for Scientific Research (B) [Grant Number 18H02518, to TM, HW, and YY] and Grant-in-Aid for Scientific Research (C) [Grant Number 18K06441, to YY].

Author information

Authors and Affiliations



TM contributed to the conception and design of the study. YY contributed to the research design; performed the molecular lab work, analysis, and interpretation of the data; and drafted the manuscript. TM and SS contributed to the cold exposure experiments, saliva sampling, and acquisition of the data. TM, HW, and SS contributed to the interpretation of the data and the revision of the manuscript. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Yoshiki Yasukochi.

Ethics declarations

Ethics approval and consent to participate

The study protocol complied with the Declaration of Helsinki and was approved by the Clinical Research Ethics Review Committee of Mie University Hospital (Approval number U2019-016) and the Ethics Committee of the Faculty of Design, Kyushu University (Approval number 269). Written informed consent was obtained from all subjects prior to their enrollment.

Consent for publication

All participants gave written informed consent for publication after a complete explanation of this study.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The quality of total RNA used in RNA-seq and real-time RT-qPCR analyses.

Additional file 2: Fig. S1.

Gene expression levels of 13 differentially expressed genes with |LFC| < 1, identified using the DESeq2 program, in each of the study subjects. Normalized counts at the vertical axis represent the estimated expression abundance at the gene level.

Additional file 3: Fig. S2.

Relative quantification levels of CTSL measured using real-time RT-qPCR and the comparative Ct method. The bold black bars represent median values. The vertical axis represents the mean relative quantification of CTSL transcripts.

Additional file 4: Fig. S3.

Hierarchical tree of GO terms in the biological process category. Significantly overrepresented GO terms (FDR < 0.05) are shown in red characters.

Rights and permissions

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 the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other 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 this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yasukochi, Y., Shin, S., Wakabayashi, H. et al. Upregulation of cathepsin L gene under mild cold conditions in young Japanese male adults. J Physiol Anthropol 40, 16 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: