The skeletal remains of 13 adult Yàmanas from Tierra del Fuego were selected from the collection of the Museum of Anthropology "G. Sergi" of the Sapienza University of Rome (Italy), which includes fifteen complete or almost complete skeletons (13 adults, 1 juvenile, 1 infant) in good state of preservation. Great part of the collection was recovered in 1883-1884 by the Italian explorer Giacomo Bove during one of his expeditions in South America, and then acquired by the University of Rome (1886) through a donation for the rising museum curated by Giuseppe Sergi . The selection criteria included the availability of femur and lumbar vertebras and estimated age of 25 years old or older as proof of achievement of adult bone structure.
We retrospectively selected a population of 217 patients expressing F-FDG BAT among a total of 6454 patients who underwent 8004 consecutive 18F-FDG PET/CT scans from January 2007 to June 2010 at the Istituto Nazionale dei Tumori Regina Elena (Rome, Italy, 41.81°N 12.45°W). The modality of BAT expression detection is described elsewhere. Among the 217 patients with BAT depots, individuals younger than 25 years old or who underwent PET/CT scans for malignant diseases and received the last treatment for malignant diseases within one year before the scan were excluded, as where subjects whose diagnosis and/or treatment timings were unknown. We also excluded those with diabetes, renal failure, and under steroids or osteoporosis treatment at the time of PET/CT. Moreover, a gender, age and BMI matched control population of BAT negative subjects was selected among the same population following the same exclusion criteria stated above. To control for BMI difference, given the reported lean body mass of the Fuegians, we only selected patients with a BMI ≤25 kg/m2.
Among the patients enrolled in the present study, a subpopulation was subsequently selected to be age, BMI and gender matched to those estimated for the fuegian skeletal remains adopting the following criteria: available Dual energy X-ray absorptiometry (DXA) scan performed within 1 year of the PET scan with a Hologic scanner.
Positron emission tomography/computed tomography (PET/CT) of living humans
PET/CT scans were performed using a Biograph 16 High Rez PET/CT scanner (Siemens AG, Munich, Germany). Patients came from the temperate metropolitan area of Rome and were instructed to fast overnight for at least 12 h before the scan and to abstain from carbohydrates and very fatty foods consumption, nicotine, caffeine, or alcohol intake for the preceding 24 hours. Room temperature water was allowed at all times. Since their arrival to the hospital, the patients were in an air-conditioned environment at about 22°C. After the injection of 18F- fluorodeoxyglucose (FDG) the patients rested at 24°C for 1 hour before undergoing the PET/CT scan. Data on gender, BMI, age, plasma glucose level and malignancy status (active: PET/CT scan positive for malignancy, not active= PET/CT scan negative for malignancy) were obtained for all patients. Parameters measured from PET/CT scans included, in addition to BAT characteristics: cross-sectional area (CSA, mm2), cortical bone area (CBA, mm2), endocortical bone area (EndA, mm2), muscle area (MA, mm2) taken at proximal-shaft diaphysis, 1 cm distal to the lesser trochanter. These were calculated with the use of ImageJ 1.52i, Wayne Rasband, National Institutes of Health, USA.
Computed tomography (CT) of Fuegians femurs
For each Fuegian skeleton, the right femur was analyzed using a commercial CT scanner (Revolution, GE Healthcare). Before the scan, femurs were oriented in standardized anatomical planes. Diaphyseal structural properties were calculated from DICOM files with the use of ImageJ 1.52i (Wayne Rasband, National Institutes of Health, USA). Parameters measured included: cross-sectional area (CSA, mm2), cortical bone area (CBA, mm2), and endocortical bone area (EndA, mm2) taken at proximal-shaft diaphysis, 1 cm distal to the lesser trochanter.
DXA of femur and lumbar spine of Fuegians
Lumbar spine (LS, L1-L4), femoral neck (FN) and total hip (TH) BMD of Fuegians femur and lumbar vertebrae remains were measured through DXA (QDR Discovery Acclaim, Hologic Inc., Waltham, MA, USA). Before DXA scanning, femurs of Fuegians were oriented in standardized anatomical planes and the lumbar vertebrae placed in a rack that allowed them to be aligned, using rice as a soft tissue proxy.
Alignment data from Raghavan et al., 2015  were downloaded from http://www.cbs.dtu.dk/suppl/NativeAmerican/. Data corresponding to 8 Fuegians and 14 distinct Native Americans individuals exposed to cold temperatures were identified and extracted. Read statistics per sample and alignment methods are shown in the Supplementary Materials for Raghavan et al. .
Since the individual Fuegians libraries have a variable endogenous content (from 0.6% to 23.8% of the total number of reads mapped to the human genome) and a low sequencing depth (average depth from 0.003 to 1.7) , samples with less than 5 million of mapper reads were removed prior to further analyses. The remaining 5 fuegian individuals included 3 Yaghan (~47M of average mapped reads and ~21% of endogenous content), 1 Selknam (81.5M of mapped reads and ~16% of endogenous content) and 1 Kaweskar (~15M of mapped reads and ~1.4% of endogenous content).Genotypes were called both per-sample and in a multi-sample approach by using GATK (v.184.108.40.206)  and modelled after the GATK Best Practices Workflows and the parameters used in Raghavan et al., 2015. Briefly, variants were called per-sample using HaplotypeCaller in GVCF mode and prepared for filtering (tools involved: CombineGVCFs, GenotypeGVCFs). Variants were extracted, filtered and ricalibrated for QualByDepth, RMSMappingQuality, MappingQualityRankSumTest, ReadPosRankSum, FisherStrand, ReadPosRankSumTest and HaplotypeScore using the HapMap 3.3, and dbSNP138 resources with priors 15 and 2, respectively.
The filtered variants were imported and annotated using R packages vcfR (v.1.12.0), GenomicFeatures (v.1.42.1) and VariantAnnotation (v.1.36.0) and the gencode annotation v.19 as reference. SNPs located in coding, promoter and 5’/3’ UTR regions of a panel of 28 genes directly involved in BAT metabolism from Sazzini et al., 2014 , were selected and further investigated (S1).
Candidate variants were further characterized based on possible functional effect as predicted by in silico analysis using MutationTaster and RegulationSpotter .These are online applications performing several in silico tests on both DNA and protein level ultimately estimating the impact of the identified variant, such as functional consequences of synonymous or intronic mutations up to amino acid substitutions, deletion and insertion of sequences, or variants within the intron-exon border.
Data relative to the study cohort of living humans was collected from a database previously used to report BAT prevalence in central Italy . The project had been approved by the Institutional Review Board of the Dipartimento di Prevenzione e Diagnostica Oncologica, UOC di Medicina Nucleare, Istituto Nazionale Tumori Regina Elena, Roma. The data were analyzed anonymously. Written informed consent was given by the patients or from the next of kin or caretakers for their information to be stored in the hospital database and used for research. All methods were carried out in accordance with relevant guidelines and regulations.
Statistical analysis was performed using SPSS 25.0 (SPSS, Inc., Chicago, IL). Data are expressed as mean ± standard deviation for normally distributed variables. Variables were tested for normality of distribution using Shapiro-Wilk test. Variables that were not normally distributed were log-transformed. Comparisons between groups were performed using the paired Student's t-test or Wilcoxon test as appropriate. Relationships between variables were measured by Pearson’s correlation coefficient. Stepwise regression modeling was used to determine predictors of CBA. A chi-square test was performed to compare variant frequencies between Fuegians, Siberians/North Americans, and other modern populations from gnomAD database (https://gnomad.broadinstitute.org/). An α error of 0.05 was considered the threshold for statistical significance