The present study was approved by the Institutional Review Board of the Faculty of Medicine, Chulalongkorn University (Bangkok, Thailand (approval no. 264/62) and conducted in accordance with the 1964 Helsinki Declaration and its later amendments, as well as all relevant ethical guidelines and regulations. The rights and interests of all participants were fully protected. Fifteen-year-old monozygotic twin girls with generalized hyperpigmentation presented in April 2019 at Khanu Woralaksaburi Hospital, Kamphaeng Phet, Thailand were referred to Naresuan University Hospital, Phitsanulok, Thailand. Both parents of the patients were aged 56 years. Skin biopsy was performed on the abdomen of both twins for histopathological study, diagnosis and fibroblast culture. Written informed consent was obtained from the parents of the patients. Tissue was fixed in 10% formaldehyde at room temperature for 24 h, followed by tissue processing which included dehydration through a graded ascending ethanol series, clearing in xylene, and embedding in paraffin. Paraffin-embedded tissues were then sectioned at 4-5 µm thickness for histological analysis. Hematoxylin and eosin staining was performed at room temperature following deparaffinization in xylene, rehydration through graded descending ethanol, staining with Harris's hematoxylin for 7.5 min, bluing with lithium carbonate, counterstaining with eosin for 1 min, followed by dehydration and clearing. Sections were imaged using an Olympus BX50 light microscope (Olympus Corporation), with image capture via Olympus cellSens Standard software v1.3 (evidentscientific.com/en/products/software/cellsens).

As part of the clinical work-up, both twins underwent a physical examination, electrocardiogram (ECG), and routine laboratory testing, including a complete blood count (CBC) and morning serum cortisol assessment, performed at Naresuan University Hospital, Phitsanulok, Thailand. A standard-dose (250 µg) adrenocorticotropic hormone (ACTH) stimulation test was performed in the morning between 08:00 and 09:00 to account for diurnal variation in cortisol secretion. A baseline blood sample was collected for measurement of serum cortisol and 17-hydroxyprogesterone prior to ACTH administration. Subsequently, 250 µg of cosyntropin was administered intravenously, followed by blood sampling at 30 and 60 min post-injection to assess stimulated cortisol and 17-OHP levels. Serum cortisol levels were analyzed using the electrochemiluminescence immunoassay method. The reference range was 6.2-19.4 µg/dl in the morning and 2.3-11.9 µg/dl in the afternoon. Serum 17-hydroxyprogesterone levels were analyzed using the radioimmunoassay (RIA) method, with a reference range of 0.10-1.50 ng/ml.

Genomic DNA was obtained from both twins and their parents, extracted from leucocytes using a Gentra Puregene Blood kit (Qiagen GmbH) and sheared using a Megaruptor 2 long hydropore to obtain DNA fragments of 15 kb. DNA fragments were analyzed using a Qubit fluorometer and pulsed-field gel electrophoresis. DNA fragments (2-5 µg) were prepared using the SMRTbell Express Template Prep kit 2.0 (cat. no. P/N 100-938-900) and SMRTbell Enzyme Clean-up kit 2.0 (cat. no. P/N 101-932-600; both Pacific Biosciences). Fragments <10 kb were eliminated using BluePippin (Sage Science). SMRTbell libraries were sequenced on the Sequel II system using Sequel II sequencing kit 2.0 (cat. no. P/N 101-820-200; Pacific Biosciences) and raw subreads were processed through the circular consensus sequencing (CCS) workflow using PacBio SMRTLink version 10.0 (pacb.com/products-and-services/analytical-software/smrt-analysis) to generate high fidelity (HiFi) long reads, in which both DNA strands were sequenced multiple times due to the circular nature of the template, achieving a minimum estimated quality value of 20 (phred-scaled, corresponding to an accuracy of 99%). All samples were aligned to the GRCh38 human reference genome using pbmm2 v1.9.0 (github.com/PacificBiosciences/pbmm2). The aligned CCS binary alignment map data were used for small variant, SV and copy number variant (CNV) detection.

FASTQ files, generated from long-read whole-genome sequencing, were aligned to GRCh38 using LAST v1418 tool (18), yielding a MAF file. The MAF file was then used by tandem-genotypes v1.1.0(19) to predict STR unit change of each read in relative to repeat unit in the reference genome. The repeat unit change was estimated in the twins, parents and 236 unaffected controls from the in-house database at Excellence Center for Genomics and Precision Medicine, King Chulalongkorn Memorial Hospital, Bangkok, Thailand. This database contains long-read sequencing data from patients with various genetic conditions, as well as from parents or siblings of certain probands. Importantly, the 236 samples designated as controls did not include individuals with pigmentation disorder. Next, STRs were identified using an in-house algorithm. Repeat patterns with 2-6 bases were selected to calculate distributions of repeat unit using kernel density estimation. Gaussian kernel with a bandwidth of five was selected to smooth the estimated density curve. Minimum points were used at the curve to cluster the data. For autosomal dominance mode, STRs exhibiting expansion only in the twins but not in their parents and controls, were searched. For autosomal recessive mode, expansion in both twins and parents but not in controls were searched. STRs were further reviewed for their relevance to the disease by examining genomic context, literature on disease association and predicted effects on gene function or expression. Analysis of single nucleotide variants (SNVs), insertions and deletions (Indels), SV, CNV and methylation profile is described in the Data S1.

Kinship coefficient was estimated using kinship-based inference from the KING package v2.3.2, based on genotype data derived from whole-genome sequencing (20). Relatedness was inferred based on estimated kinship coefficients using the following classification: Kinship coefficient >0.354 indicates duplicates or monozygotic twins, where (0.177, 0.354), (0.0884, 0177) and (0.0442, 0.0884) indicate first-second- and third-degree relatives, respectively.

Fibroblast cells were acquired from skin biopsies of two affected twins carrying the TYMS variant and four unaffected individuals who were not in the in-house database. The controls included three males and one female, all pediatric donors aged 1-5 years, with no history of pigmentation disorder or related skin conditions. The cells were cultured in DMEM supplemented with 10% fetal bovine serum (both Cytiva; HyClone) and 1X antibiotic-antimycotic (Gibco; Thermo Fisher Scientific, Inc.) and then incubated at 37˚C with 5% CO₂ (21). The use of primary human cells in cell culture experiments was approved by the Institutional Review Board of the Faculty of Medicine, Chulalongkorn University (approval no. 813/63).

RNA was extracted from fibroblasts and peripheral blood mononuclear cells (PBMCs) using the QIAamp RNA Blood Mini kit (Qiagen GmbH). PBMC controls included both parents of the twins and one unrelated 45-year-old adult female, all without any history of pigmentation disorders or related skin conditions. These PBMC controls were independent from the aforementioned fibroblast controls. cDNA was prepared from 1 µg RNA using the ImProm-II™ RT System (Promega Corporation) according to the manufacturer's protocol. The input RNA from fibroblasts was primed with an equal mix of anchored dT oligonucleotides for reverse-transcribing RNA to cDNA. TaqMan probes specific for TYMS (cat. no. Hs00426586_m1) and β-actin (Assay ID: Hs01060665_g1) were used, and all reactions were set up with Luna® Universal Probe qPCR Master Mix (New England Biolabs). A total of three replicates/sample were run for both TYMS and β-actin on the Applied Biosystems StepOnePlus Real-Time PCR system (Thermo Fisher Scientific, Inc.) and the results were analyzed using StepOne Software v2.3 (Thermo Fisher Scientific, Inc.). qPCR thermocycling conditions were as follows: Initial denaturation at 95˚C for 60 sec, followed by 40 cycles of denaturation at 95˚C for 15 sec and extension at 60˚C for 30 sec. The expression of the gene of interest (TYMS) was normalized against the control gene (β-actin), and the relative change between affected and controls was analyzed using the 2^-ΔΔCq method (22).

To assess TYMS expression, western blot analysis was performed. Briefly, protein lysate was prepared from fibroblast cells using radioimmunoprecipitation assay buffer supplemented with protease inhibitors (Cell Signaling Technology, Inc.). The protein concentration was determined using a Bicinchoninic Acid Protein Assay kit (Thermo Fisher Scientific, Inc.). Equal amounts of protein (50 µg/lane) were separated by 12% SDS-PAGE. The proteins were transferred onto a PVDF membrane using iBlot™ 2 Transfer Stacks (Invitrogen; Thermo Fisher Scientific, Inc.). The membrane was blocked with 5% non-fat dry milk and 1% BSA in TBS for 1 h at room temperature to prevent non-specific binding. Subsequently, the membrane was incubated overnight at 4˚C with the primary antibody against TYMS (1:500; cat. no. HL1236; Invitrogen; Thermo Fisher Scientific, Inc.). After washing with 0.1% TBST, the membrane was incubated with the corresponding anti-rabbit IgG, horseradish peroxidase-conjugated secondary antibody (1:2,000; cat. no. 7074; Cell Signaling Technology, Inc.) for 3 h at room temperature. TYMS protein bands were visualized at 34 kDa using an enhanced chemiluminescence detection SuperSignal™ West Femto Maximum Sensitivity Substrate (Thermo Fisher Scientific, Inc.), and chemiluminescent signals were visualized and quantified using the ImageQuant™ LAS 4000 system (version 1.3; GE Healthcare Life Sciences). The relative protein expression levels were normalized to the housekeeping gene β-actin (13E5) Rabbit mAb (1:1,000, Cat. No. 7087; Cell Signaling Technology, Inc.).

The twins had a kinship coefficient of 0.4934, confirming their monozygosity and ensuring that the sequencing data accurately corresponded to them. Both twins exhibited progressive skin darkening since birth, with strikingly similar pigmentation patterns. Both twins were healthy without any other symptoms. They were born to non-consanguineous parents with normal skin color and no family history of hyperpigmentation disorder. No history of drug intake or urine discoloration or photosensitivity was reported. The date of admission was April 2019 for adrenal insufficiency workup, which yielded a negative result. The patients were referred to a dermatologist. The first skin biopsy was performed in September 2019 for histopathological examination, followed by a second biopsy ~6 months later to obtain fibroblasts for cell culture. There was no subsequent follow-up.

On examination, both twins had generalized, diffuse, dark-black hyperpigmentation of the skin on the face, neck, body, hands and feet (Fig. 1A-D) and some reticulate hypopigmentation and hyperpigmentation were found on the abdomen. Hyperpigmentation was prominent on the distal finger and feet. Tongue and buccal mucosa also exhibited some hyperpigmentation. The hair and nails were normal.

Skin biopsy from the right trunk revealed increased melanin pigmentation in basal and suprabasal layers and presence of melanophages in upper dermis. These melanophages were identified as large, round or oval cells with abundant cytoplasmic melanin granules, giving them a dark brown/black appearance under the microscope. Additionally, no melanocyte proliferation was observed, based on the absence of increased melanocyte density or atypical melanocytic nests in the basal epidermis (Fig. 1E). Physical examination of both twins revealed no abnormality aside from hyperpigmentation. Electrocardiogram results were normal. Complete blood count showed mild anemia with hematocrit 31.7% and no macrocytosis. Cortisol levels in the morning and ACTH stimulation test results were normal, thus, excluding Addison's disease.

HiFi reads were mapped on the GRCh38 human reference genome. The mean breadth of coverage of the four samples was 97.19% (range, 95.08-98.94%). The average sequencing depth of the four samples was 28X (range, 26-31X).

Variants were filtered based on inheritance pattern. The number of concordant SNVs (twins have the same genotype) and Indels between twins from DeepVariant was 7,030,746. There were 29,137 autosomal dominant de novo variants, 236,591 autosomal recessive variants and 28,108,264 pairs of compound heterozygous variants. In the first approach, two de novo and two pairs of compound heterozygous variants were classified as likely pathogenic and frameshift mutations, according to Franklin by Genoox. However, those variants were regarded as false positives after visually inspecting BAM files with Integrative Genomics Viewer v2.15.4(23).

In the second approach, one autosomal recessive was identified in phosphatidylinositol glycan anchor biosynthesis class B (PIGB) (c.981T>G), which was a missense variant. This variant was absent in the in-house database but present in controls in gnomAD database. The gene constraint score of the c.981T>G variant indicated a tolerance to variation. Based on the guidelines of the American College of Medical Genetics and Genomics (ACMG), PIGB (c.981T>G) variant was classified as benign. Moreover, eight compound heterozygous variants which were missense variants in myosin heavy chain 13 (MYH13) (c.53G>A and c.791A>G), serpin family H member 1 (SERPINH1) (c.293G>A and c.1066G>A), myosin VIIA (MYO7A) (c.781G>A and c.4349A>G), UDP glucuronosyltransferase family 1 member A1 (UGT1A1) (c.1459C>T and c.964A>G), keratin 36 (KRT36) (c.122G>A and c.677G>T), tectonin β-propeller repeat containing 1 (TECPR1) (c.2270C>A and c.577G>A), tripartite motif containing 10 (TRIM10) (c.842G>A and c.977C>T), and missense and intronic variants in cholinergic receptor nicotinic epsilon subunit (CHRNE) (c.1375T>C and c.917+562T>C) were identified (Table SI). Despite their absence in the in-house database and controls in gnomAD database, the gene constraint score of all variants indicated tolerance to variation. Based on ACMG guidelines, compound heterozygous variants were classified as variants of uncertain significance.

SV analysis. The number of concordant SVs between twins from pbsv was 73,948 variants, of which 69,106 passed quality control and were filtered based on inheritance pattern, yielding 157 de novo and 2,267 autosomal recessive variants. No compound heterozygous variants were retained. In the first approach, no pathogenic or likely pathogenic variants were found. In the second approach, two autosomal recessive variants in the solute carrier family 22 member 31 (SLC22A31) and a region from sialic acid binding Ig like lectin 14 (SIGLEC14) to SIGLEC5 were obtained. After inspecting BAM files, these SVs were true positives; however, they were not rare (frequency >1% in the in-house genome database) in Thai individuals. Therefore, no SVs were considered to be associated with the phenotype.

CNV analysis. No CNVs were phenotype-associated. Copy numbers that were different from both parents were identified in 13 regions. After annotating the regions with ClassifyCNV, eight regions were predicted as uncertain significance (four of which were located within or overlapped protein-coding genes) and five were predicted as benign. After checking the copy number in these regions in the in-house database, all regions were ruled out since allele frequency was >1%.

The present study identified 211 DMRs, of which 161 resided in 143 genes and 50 DMRs were not in genes. The annotation showed that 51.66, 18.01, 12.32, 10.90, 5.69 and 1.42% were in the promoter (≤2kb), distal intergenic, intron, exon, promoter (2-3 kb) and 3' untranslated region. respectively. After performing biological function prediction, 35 significant terms were identified (Table SII). However, the significant enrichments were not related to skin pigmentation or melanogenesis.

After starting with 718,245 repeat patterns, filtering for STR motifs 2-6 base pairs in length, resulting in 240,677 patterns. Results from the in-house algorithm returned eight candidate patterns. Only one pattern (GATGGT) was consistent with autosomal recessive inheritance, whereas other patterns were not consistent between phenotype and mode of inheritance. There was a biallelic GATGGT repeat expansion with a range of 210-259 repeats in intron 3 of TYMS in both twins (Fig. 2A), whereas the parents exhibited a heterozygous expansion. The father had 106 repeat units in one allele and 230-245 repeat units expanded in the other allele. The mother had 93 repeat units in one allele and 217-224 repeat units expanded in the other allele compared with the reference genome (data not shown). The 236 controls had ≤172 repeat units (Fig. 2B). The majority of controls (47%; 112/236) had 106 repeat units. All controls with 172 repeat units carried the expansion on a single allele only. No homozygous controls were identified, unlike in the twins.

RT-qPCR was used to investigate the expression levels of TYMS in identical twin patients. RNA was extracted from cultured skin fibroblasts and PBMCs. The mean TYMS expression levels in the patients (PBMCs, 0.74 and 0.76; skin fibroblasts, 0.91 and 1.03) were within the range observed in control samples (PBMCs, 0.30-1.00; skin fibroblasts, 0.47-1.41; Figs. 3A and S1). No apparent differences were observed in either cell type. Protein was extracted from cultured skin fibroblasts. The western blot analysis showed that TYMS protein levels varied between the two probands, with twin A showing lower expression than twin B. Overall, the levels observed in both twins were within the range seen in control samples (Fig. 3B, S2 and S3).