Unexplored regulatory sequences of divergently paired GLA and HNRNPH2 loci pertinent to Fabry disease in human kidney and skin cells: Presence of an active bidirectional promoter.

Fabry disease (FD) is a rare hereditary disorder characterized by a wide range of symptoms caused by a variety of mutations in the galactosidase α (GLA) gene. The heterogeneous nuclear ribonucleoprotein (HNRNPH2) gene is divergently paired with GLA on chromosome X and is thought to be implicated in FD. However, insufficient information is available on the regulatory mechanisms associated with the expression of HNRNPH2 and the GLA loci. Therefore, the current study performed bioinformatics analyses to assess the GLA and HNRNPH2 loci and investigate the regulatory mechanisms involved in the expression of each gene. The regulatory mechanisms underlying GLA and HNRNPH2 were revealed. The expression of each gene was associated with a bidirectional promoter (BDP) characterized by the absence of TATA box motifs and the presence of specific transcription factor binding sites (TFBSs) and a CpG Island (CGI). The nuclear run-on transcription assay confirmed the activity of BDP GLA and HNRNPH2 transcription in 293T. Methylation-specific PCR analysis demonstrated a statistically significant variation in the DNA methylation pattern of BDP in several cell lines, including human adult epidermal keratinocytes (AEKs), human renal glomerular endothelial cells, human renal epithelial cells and 293T cells. The highest observed significance was demonstrated in AEKs (P<0.05). The results of the chromatin-immunoprecipitation assay using 293T cells identified specific TFBS motifs for Yin Yang 1 and nuclear respiratory factor 1 transcription factors in BDPs. The National Center for Biotechnology Information-single nucleotide polymorphism database revealed pathogenic variants in the BDP sequence. Additionally, a previously reported variant associated with a severe heterozygous female case of GLA FD was mapped in BDP. The results of the present study suggested that the expression of the divergent paired loci, GLA and HNRNPH2, were controlled by BDP. Mutations in BDP may also serve a role in FD and may explain clinical disease diversity. Copyright: © Al-Obaide et al.

Introduction

Fabry disease (FD; OMIM 301500) is an X-linked inherited error lysosomal storage disorder, caused by mutations in the galactosidase α (GLA) gene that encodes the α-galactosidase A enzyme (α-GAL A) (1-6). Patients with α-GAL A deficiencies exhibit a harmful accumulation of globotriaosylceramide (GL-3/Gb3) in capillary endothelial, renal (podocytes, tubular, glomerular endothelial, mesangial and interstitial cells), cardiac (cardiomyocytes and fibroblasts) and nerve cells (7-10). It is hypothesized that the accumulation of GL-3 causes a wide spectrum of pathogenic symptoms that include progressive renal insufficiency, cardiac involvement, and neuropathology (11-14). End-stage renal disease and life-threatening cardiovascular or cerebrovascular complications limit the life-expectancy of untreated patients by 20 and 10 years, respectively, compared with the general population (1). Almost all FD-associated complications are non-specific and clinically indistinguishable from similar abnormalities that occur in the context of more common disorders, for example neuropathic pain, sweating, gastrointestinal symptoms, and pulmonary symptoms (15-21). FD is treated via enzyme replacement therapy (ERT), which substitutes the missing or altered, partially functional α-GAL A (18). Additionally, DGJ (1-deoxygalactonojirimycin), a pharmacological chaperone, is used to treat amenable α-GAL A missense mutations with adverse side effects (20). The aforementioned therapeutic techniques do not reverse all symptoms experienced by patients with FD (18,19,22). Therefore, it is undetermined whether FD symptoms are solely associated with the accumulation of GL-3 or whether other genetic mechanisms may be involved.

The current study hypothesized that the causes of FD are not limited to α-GAL A enzyme malfunction as a result of pathogenic GLA genetic mutations. The GLA locus is paired in a divergent configuration with heterogeneous nuclear ribonucleoprotein (HNRNPH2) on chromosome X, sharing regulatory sequences at the 5'-ends. The National Center for Biotechnology Information (NCBI) database has previously suggested that HNRNPH2 (Gene ID: 3188) is a likely cause of FD. Previous studies have demonstrated that divergent genes are governed by bidirectional promoters (BDP), which are characterized by the absence of a TATA box, the presence of several transcription factor motifs including nuclear respiratory factor 1 (NRF1), Yin Yang 1 (YY1), GA binding protein transcription factor subunit α (GABPA) and zinc finger protein 143 (ZNF143), an abundance of GC and CpG islands (CGIs) (23-29). Divergent mutations in the BDPs of paired genes are associated with cancer (30-33), cerebral cavernous malformations (34) and nicotine initiation or addiction (35). Although heterozygous females with FD usually present an attenuated form of the disease, a recent study reported a severe case of FD observed in the GLA of a heterozygous female, likely occurring from a mutation that altered the CGI methylation status (36). The location of the aforementioned mutation is in the discovered BDP reported in this study. The current study therefore performed bioinformatics and additional experiments to investigate the potential presence of BDP and its association with GLA and HNRNPH2 expression.

Materials and methods

Genomics databases

Genomic features of the divergently paired GLA and HNRNPH2 loci were searched in the following genomic databases: Ensembl (https://uswest.ensembl.org/index.html), NCBI-Gene (https://www.ncbi.nlm.nih.gov/gene/?term=), GeneCards (https://www.genecards.org/), EMBL-EBI (https://www.ebi.ac.uk/), UCSC genomic browser (https://www.genome.ucsc.edu/), PrESSTo (http://pressto.binf.ku.dk/about.php) and the Eukaryotic promoter database, EPD (https://epd.epfl.ch//index.php). The precise genomic map locations of the identified sequences were verified and uploaded to the hg38 human genome sequence using the ‘get DNA’ and ‘BLAT’ tools of the UCSC Genome Browser database. EMBOSS Needle (https://www.ebi.ac.uk/Tools/psa/emboss_needle/) was used for the pairwise sequence alignment of nucleotide sequences.

Regulatory sequences at the 5'-ends of GLA and HNRNPH2

The 5'-end genomic regions of the divergently paired GLA (ENSG00000102393) and HNRNPH2 (ENSG00000126945) loci were searched for regulatory sequences. The motif tool of the EPD database was used to search for TATA box motifs. Additionally, the JASPAR database (http://jaspar.genereg.net/) of transcription factor binding profiles (37,38) was used to determine the binding site motifs of YY1, NRF1, E2F transcription factor 1 (E2F1) and GABPA transcription factors in the BDP (23-29). CGI identification and plotting in the predicted BDP sequences were achieved using the EMBOSS cpgplot tool (https://www.ebi.ac.uk/Tools/seqstats/emboss_cpgplot/). The parameters used to search for CGIs were as follows: Obs/Exp CpG >0.6, G%+C% >50% and 100 nucleotides in length. The details of CGI prediction are included in the EMBOSS cpgplot manual.

Human cell lines, RNA and DNA

293T (ATCC® CRL-1573™) cells were purchased from The American Type Culture Collection and cultured in Eagle's minimum Essential medium (cat. no. 30-2003) supplemented with 10% heat-inactivated FBS (ATCC® 30-2020™), 100 U/ml penicillin and 100 µg/ml streptomycin. Cells were cultured in a humidified incubator in the presence of 5% CO2 at 37˚C. A purification kit (cat. no. 48700; Norgen Biotek Corp.) was used for the extraction of RNA and DNA from 293T cells. Human genomic DNA for methylation analysis isolated from adult epidermal keratinocytes (cat. no. 2119), renal glomerular endothelial cells (cat. no. 4009) and renal epithelial cells (cat. no. 4129) were purchased from ScienCell Research Laboratories, Inc.

Primer set design and quantitative (q)PCR analysis

Primer sets for qPCR-intercalating dyes were purchased from Integrated DNA Technologies, Inc. The PrimerQuest Tool was used for the custom design of nuclear run-on (NRO) and chromatin-immunoprecipitation (ChIP) assay primers. The methylation primer sets were designed with the MethPrimer program (39) and used for methylation specific-PCR (MSP) analysis. The real-time qPCR reactions were performed using the iTaq universal SYBR-Green reaction mix on a Bio-Rad CFX96 Real-Time system. The thermocycling conditions were as follows: Polymerase activation and DNA denaturation at 95˚C for 3 min followed by 35 (MSP assay) and 45 cycles (NRO or ChIP assays) of denaturation at 95˚C for 5 sec, and annealing and extension for 30 sec at 60˚C (NRO or ChIP assays) and 58˚C (MSP assay), respectively. Replicate PCRs were run on the same sample where target and reference are amplified in separate wells (40). At least three independent experiments were performed for each sample (41).

ChIP analysis

ChIP was carried out using the EpiQuik™ Chromatin Immunoprecipitation kit (Epigentek Group Inc.) in accordance with the manufacturer's protocol. 293T cells were fixed using 1% formaldehyde solution for 10 min at room temperature and sonicated samples were subsequently immunoprecipitated at room temperature for 90 min using NRF1 (ab175932) and YY1 (ab38422) antibodies purchased from Abcam. Normal mouse IgG (1 mg/ml) (negative control) and 1 mg/ml anti-RNA Polymerase II (positive control) were provided as part of the kit; 1 µl of normal mouse IgG, 1 µl of anti-RNA polymerase II, and 5 µg of anti-NRF1 and anti-YY1 were added to each well. The wells were covered with Parafilm M and incubated at room temperature for 120 min. Protein-DNA complexes were de-crosslinked and the DNA was purified using the F-Spin column. The primer sets were designed to produce amplicons that host the NRF1 motifs (forward, 5'-AGCTGAGGAACCCAGAACTA-3' and reverse, 5'CAATCCATTGTCCAGTGCTCTA-3') and YY1 motifs (forward, 5'GTCATGAGCGTCCACCATTT-3' and reverse, 5'-CCTCTTTCGTTCTCTGCTTTCC-3'). The Fold enrichment method was used to analyze ChIP-qPCR data. The 2-∆∆CT value was calculated from (Ct IP)-(Ct mock). Normalization was completed using the IgG Ct value. Duplicate PCRs were run on the same sample and Ct average data were used (40). Four independent experiments were performed for each sample (41).

Bisulfite DNA treatment and MSP analysis

The methylamp DNA Modification kit (Epigentek Group Inc.) was used to investigate the methylation status of predicted BDP sequences in accordance with the manufacturer's protocol. Genomic DNA from adult epidermal keratinocytes (cat. no. 2119), renal glomerular endothelial cells (cat. no. 4009), renal epithelial cells (cat. no. 4129) and 293T cells was used for MSP analysis. Purified genomic DNA (100 ng) was treated using the Methylamp DNA Modification kit, after which the converted DNA was cleaned, captured and eluted using R6 (Modified DNA Elution) solution and an F-Spin column. Eluted DNA was analyzed using the iTaq universal SYBR-Green reaction mix (Bio-Rad Laboratories, Inc.). The methylated primer sets were designed using the MethPrimer program (39) and used for CGI methylation analysis in the BDP sequence. The MSP primers for methylated and unmethylated regions of the 323 bp CGI-2 were as follows: M pair (forward, 5'-TTTTTTTAAACGGTTATAGCGAGAC-3' and reverse, 5'-CTTAATTTACCAAATAACCCGTA-3'), U pair (forward, 5'-TTTTTTAAATGGTTATAGTGAGATGG-3' and reverse, 5'-AATACAACACCTTAATAATCCCAAA-3'). The qPCR was performed as aforementioned. The percentage of sample methylation was calculated using the following equation: Percentage methylation =100/[1+2∆Ct (meth-unmeth)]. ∆Ct (meth-unmeth) was calculated by subtracting the Ct values of methylated CGI signals from the Ct values of the unmethylated CGI signal (42,43). Each sample was run in duplicate for qPCR analysis (40). Three independent experiments were performed for each sample (41).

NRO analysis

NRO is the method of choice when measuring the transcriptional activity of nuclear nascent mRNA transcripts (44). NRO analysis was used to investigate the expression and the quantification of GLA and HNRNPH2 nascent RNA transcription. 293T cells were used for the preparation of nuclei. The preparation of cell cultures, nuclei collection NRO transcription, nuclear RNA extraction, immunoprecipitation of bromouridylated NRO-RNAs, nascent nuclear RNA extraction, cDNA preparation and NRO cDNA quantification was performed using the bromouridine immunocapture nuclear run-on RT-qPCR method (44). NRO transcription was performed in fully resuspended nuclei by gentle pipetting at 30˚C for 30 min using the transcription buffer reaction cocktail master mix, which was composed of BrUTP, ATP, GTP, CTP and UTP. A total of 4,000,000 nuclei were used per sample for NRO transcription. Nuclear RNA was extracted using the MEGAclear transcription clean-up kit (Thermo Fisher Scientific, Inc.; cat no. AM1908) according to the manufacturer's protocol. Immunoprecipitation of bromouridylated NRO-RNAs was performed using 2 µg/tube of mouse monoclonal IgG1 anti-BrdU antibody IIB5 (sc-32323) (Santa Cruz Biotechnology, Inc.). Nascent nuclear RNA was extracted using the RNAzol method (GeneCopoeia, Inc., cat. no. E01010A) and precipitated by adding an equal volume of isopropanol and 20 µg RNase-free glycogen (1 µl of 20 mg/ml stock) to each sample. The High-Capacity cDNA Reverse Transcription kit (Thermo Fisher Scientific, Inc.; cat. no. 4368814) was used for the quantitative conversion of the extracted NRO RNA to single-stranded cDNA in a single 20 µl reaction by the thermal cycler with the following thermocycling conditions: 10 min at 25˚C, 30 min at 37˚C, 5 min at 85˚C. Amplified DNA was prepared from NRO-cDNA using the designed primer sets for HNRNPH2-E1 (forward, 5'-AGTAGTTCTGGTCGTCGT CTA-3' and reverse, 5'-ACACACCAACCTCTAACGATAC-3') and GLA-E1 (forward, 5'-AGGTTACCCGCGGAAATTTAT-3' and reverse, 5'-GAAACGAGGGCCAGGAAG-3'). Normalization was performed using hypoxanthine phosphoribosyltransferase 1 (HPRT1) primer sets (forward, 5'-TGAGGATTT GGAAAGGGTGT-3' and reverse, 5'-GAGCACACAGAGGG CTACAA-3'). The normalization and relative expression of target sequences was determined using HPRT1 as a reference gene. Duplicate PCRs were run on the same sample and Ct average data were used to calculate the 2-ΔΔCT value (40). Three independent experiments were performed for each sample (41).

Statistical analysis

Microsoft Excel 2010 (Microsoft Corporation) and GraphPad Prism 7.01 (GraphPad Software, Inc.) were used to analyze the data of various parameters. Two-sample comparison was determined by performing a Student's t-test. One-way ANOVA with Tukey's multiple comparisons test/post hoc test was used to evaluate four independent groups simultaneously and to test statistical differences between every possible pair of all groups. P<0.05 was considered to indicate a statistically significant difference.

Results

GLA and HNRNPH2 loci

The genomic setting of the paired GLA (ENSG00000102393; NCBI-Gene ID: 2717) and HNRNPH2 (ENSG00000126945; NCBI-Gene ID: 3188) genes is divergent, with a head-to-head formation on chromosome X. The GLA locus is situated on the reverse strand, whereas the HNRNHP2 locus appears on the forward strand (Fig. 1A). The Ensembl database revealed that GLA and HNRNPH2 shared a 991 nucleotide sequence between the two loci at the 5'-end (Fig. 1B). Pairwise sequence alignment analysis performed using the EMBOSS Needle tool also determined that the shared sequence between GLA and HNRNPH2 exhibited 100% similarity.

Figure 1

Mapped locations of the divergently paired head to head genes, GLA and HNRNPH2. (A) Genomic locations of the two paired genes at Xq22.1 are presented. GLA and HNRNPH2 are composed of eight and two exons respectively. The dotted elliptical contour demonstrates the shared sequence of GLA (ENSG00000102393) and HNRNPH2 (ENSG00000126945). (B) The 991 nucleotides shared between GLA complement and HNRNPH2 reverse complement. GLA, galactosidase α; HNRNPH2, heterogeneous nuclear ribonucleoprotein.

Shared regulatory sequences between GLA and HNRNHP2 loci

The 991 nucleotide sequence shared between GLA and HNRNHP2 loci lacked a TATA box motif and was rich in binding sites for several transcription factors that are found in BDP, including YY1, NRF1, GABPA and E2F1 (Table I). Additionally, bioinformatics analysis demonstrated that the shared sequence was CG rich, containing CGIs. EMBOSS Cpgplot software revealed one CGI when the following options were searched: Window size, 100; minimum sequence length, 200 bases; minimum Obs/Exp CpG, >0.6; %C+%G, >50.00%. However, three CGIs were revealed in the shared sequence when the search tool was set to a minimum sequence length of 100 bases (Fig. 2).

Table I

Transcription factor binding sites predicted using the JASPAR tool in the GLA and HNRNPH2 bidirectional promoter.

TFBSs	Number	Similarity Score
YY1	10	0.81-0.93
NRF1	8	0.80-0.89
GABPA	10	0.80-0.88
E2F1	12	0.80-0.94

YY1, Yin Yang 1; NRF1, nuclear respiratory factor 1; FABPA, GA binding protein transcription factor subunit α; E2F1, E2F transcription factor 1.

Figure 2

Three CGIs identified in the BDP 991 nucleotide sequence. The analysis was performed using the EMBOSS Cpgplot search tool. Search options for the CGIs were as follows: Window size, 100; minimum sequence length, 100 bases; minimum Obs/Exp CpG, >0.6; percent C + percent G, >50.00. BDP, bidirectional promoter; CGIs, CpG islands.

Two promoter prediction software tools determined two transcriptional start sites (TSSs) along the 991 nucleotide sequence (Table II), located at the 5'-ends of GLA and HNRNPH2. Furthermore, the EPD eukaryote promoter database revealed two main TSSs in the GLA promoter and three main TSSs in the HNRNPH2 promoter within the predicted BDP. Bioinformatics analysis suggested that the predicted 991 nucleotide sequence had the molecular features of a BDP. The predicted BDP was therefore presumed to direct the transcription of the GLA locus in a divergent manner along with its counterpart locus, HNRNHP2.

Table II

Two TSSs predicted by two promoter prediction tools in the bidirectional promoter.

Promoter prediction tool	TSS position	Score[a]
Promoter 2.0 prediction server	200	0.632
Neural network promoter prediction	688	0.97

aNeural network output activity (score) for each input subsequence. TSS, transcription start sites.

Methylation status of CGIs in kidney cell lines

The CGIs associated with bidirectional promoters are prone to DNA methylation (25). The methylation status of the previously elucidated CGI-2 at position: 241..563 (Fig. 2) was determined in four cell types, including human adult epidermal keratinocyte cells (AEK), human renal glomerular endothelial cells (HRG), human renal epithelial cells (REC) and 293T cells. The results revealed variable levels and significant differences in the extent of DNA methylation in AEK, HRG, REC and 293T cells, the highest result of which was observed in AEKs (Fig. 3A; P<0.05). RECs exhibited the lowest level of methylation compared with HRG and 293T cells, which themselves demonstrated comparable methylation percentages. Multiple comparison results presented statistical differences between HRG and AEK, HRG and REC, AEK and REC, and REC and 293T (P<0.05), but not between HRG and 293T and between AEK and 293T (P>0.05).

Figure 3

Experimental validation of the GLA and HNRNPH2 bidirectional promoter. (A) DNA methylation assayed using the methylation specific-PCR method revealed a significant variation in CGIs DNA methylation in the predicted bidirectional promoter sequence of the tested cell lines. (B) Identification of YY1 and NRF1 motifs in the GLA and HNRNPH2 bidirectional promoter sequence, as determined by performing a ChIP assay in 293T cells. (C) The bidirectional promoter sequence mapped in the overlapping 5’ sides was composed of exon 1 GLA and HNRNPH2 sequences. (D) Bidirectional transcription of the GLA and HNRNPH2, as determined by performing a nuclear run-on assay in 293T cells, P>0.05. AEK, human adult epidermal keratinocyte; HRG, human renal glomerular endothelial cells; REC, human renal epithelial cells; YY1, Yin Yang 1; NRF1, nuclear respiratory factor 1; ChIP, chromatin immunoprecipitation.

ChIP assay predicts NRF1 and YY1 motifs in the BDP

The JASPAR 2018 tool (37,38) identified the occurrence of YY1 and NRF1 binding sites in the shared 991 nucleotide sequence (Table I). The results of the ChIP assay, which was conducted using 293T cells, identified YY1 and NRF1 binding sites in the BDP sequence (Fig. 3B). The designed primers produced amplicons (~100 bp) containing YY1 and NRF1 motifs. One-way analysis showed significant statistical differences in four tested groups, P<0.05. But multiple comparison results displayed statistical differences between every possible pair of all groups (P<0.05) but not between YY1 and NRF1, P>0.05.

Quantification of nascent GLA and HNRNPH2 mRNA transcripts by NRO

Nascent NRO bromouridine-labeled mRNA derived from the nuclei of 293T cells was measured via RT-qPCR. The primer sets used during PCR were designed to amplify the adjacent sequences specific to GLA and HNRNPH2 at the 5'-ends (Fig. 3C). The data obtained using the relative 2-∆∆CT method (40) revealed an insignificant difference in divergent transcription activity between GLA and HNRNPH2 when using primer sets designed for exon 1 of GLA and HNRNPH2, P>0.05 (Fig. 3D).

Pathogenic variants of the BDP

The NCBI-dbSNP database revealed 2,657 and 1,517 variants in the GLA and HNRNPH2 loci, respectively. Additionally, pathogenic variants were identified in the BDP sequence at chromosome X:101407622-101408612. Several of these variants are presented in Table III. For example, the position of rs104894829 (C>T) was at chromosome X:101407773, whereas rs104894847 (C>G) was at chromosome X:101407846. Hossain et al (36) reported a 37-year-old female carrier of a heterozygous mutation with severe FD symptoms. The position of the C>A variant was at chrX:101407903 in the BDP sequence.

Table III

Pathogenic variants identified in the bidirectional promotor.

Variant ID	Chromosome position	Type
NA	chrX:101407903	C>A
rs104894829	chrX:101407773	C>T
rs104894831	chrX:101407786	G>A
rs104894835	chrX:101407803	T>C
rs104894836/rs28935192
rs104894847	chrX:101407738	A>C
rs104894848	chrX:101407846	C>G
rs1569306022	chrX:101407710	C>G
rs1569306168	chrX:101407716	C>T
	chrX:101407809	A>G

Data from NCBI-dbSNP and partially adapted from Ref. (36). NA, not available; chr, chromosome.

Discussion

FD is a clinically heterogeneous, slow and progressive disease that shares common symptoms in the general population (15-21). Although patients with FD exhibit a wide range of clinical symptoms, their detailed genetic presentation has not yet been fully elucidated. Pathogenic mutations in GLA have been reported as the main genetic cause of FD (11-13). A common clinical presentation of patients with FD is neuropathic pain, with both sexes exhibiting symptoms into adulthood (45,46). However, neuropathic pain usually presents in the early years of childhood (47,48). The majority of patients with FD may experience chronic or episodic pain, which are known as FD crises or acroparaesthesiae (49-52).

There is an indirect association between pain and the alternate splicing of mRNA. Donaldson and Beazley-Long (53) reviewed the known alternative splice variants that are important to pain. They suggested that these sites may serve as appropriate therapeutic targets. Furthermore, their study demonstrated that pain may result from indirect genetic conditions that arise from defects in the alternative RNA splicing process of several genes. Proteins encoded by the HNRNP gene family, including HNRNPH2, are RNA binding proteins that are associated with the mRNA splicing process. Previous studies have demonstrated the association between HNRNP genes and various clinical symptoms, such as pain. In this respect, HNRNPH1 and HNRNPF can serve as post-transcriptional regulators of opioid receptor expression (54,55). HNRNPH2 and HNRNPF genes encode a similar protein structure and bind to similar sequences in order to influence gene expression (56). Since HNRNPH2 is a member of the HNRNP gene family, this may indicate its involvement in the pain experienced by patients with FD. The NCBI-Gene database notes that HNRNPH2 may be implicated in FD. It is hypothesized that aberrations in HNRNPH2 expression can cause defects in mRNA splicing. Nearly all human multiexon genes are subject to alternative splicing (57), which has been reported in GLA (58-60). A recent study has suggested that DNA methylation not only affects transcription, but also regulates alternative splicing (61). Accordingly, the current study hypothesized that defects in HNRNPH2 expression caused by BDP mutations can influence symptoms associated with pre-mRNA splicing malfunctions, including pain. Thus, defects in HNRNPH2 expression can alter the normal mRNA splicing process of multiexon genes including GLA and other genes associated with diseases and pain (53).

At present, there is insufficient information available on the regulatory mechanisms associated with the expression of paired, divergent HNRNPH2 and GLA loci on chromosome X. The present study provided bioinformatic and experimental data that demonstrated the occurrence of a BDP regulatory sequence at the 5'-end of GLA and HNRNPH2. Bioinformatics analysis identified a shared 991 nucleotide sequence in the BDP, which was associated with bidirectional transcription. The main features of this sequence, which are characteristic of BDPs (23-29), was the absence of TATA-box binding motifs, alongside the presence of a 323 bp CGI-2 and the occurrence of specific transcription factor motifs, including, YY1 and NFR1. The results of the present study also indicated the presence of YY1 and NFR1 transcription factors' motifs in 323 bp CGI-2. The methylation status of four human cell lines was also evaluated. The NRO assay performed in the current study confirmed the divergent expression of GLA and HNRNPH2 in the nascent nuclear transcripts of 293T cells.

Several mammalian promoters demonstrate divergent transcription. In the human genome, it is estimated that >10% of genes are divergently transcribed, where a single promoter is shared with their transcriptional start sites (23). BDPs are known to direct divergent gene expression, which is the case for GLA and HNRNPH2. Additionally, bidirectional promoters are associated with various diseases, for example cancer (30-33), cerebral cavernous malformations (34) and nicotine initiation or addiction (35). The present data also provided a reason for the unexplained clinical manifestations observed in patients with FD. Although asymptomatic females with FD are likely an exception, it has been reported that heterozygous females may suffer from significant multisystemic clinical symptoms (62).

Finally, although the results of the present study indicated the potential involvement of the BDP, the HNRNPH2 gene and the GLA gene in FD, the precise mechanism that regulates the bidirectional transcription of GLA and HNRNPH2 is yet to be understood. Additional experimental studies using novel methods, including high-throughput assays (63) are required for the better understanding of the architecture and cis-regulatory elements of GLA and HNRNPH2. Furthermore, investigation into the role of chromatin (64) on the transcriptional status of the two genes may be required. Another limitation to the present study is the lack of tissues or cells from patients with FD.

In conclusion, the molecular characteristics of the GLA and HNRNHP2 bidirectional promoter may define the role of upstream pathogenic variants of the divergently paired GLA and HNRNPH2 genes in FD. Accordingly, defects in the BDP can simultaneously disturb the expression of GLA and HNRNPH2 and cause diverse clinical manifestations associated with FD.