Genetics and Molecular Biology

Reduced biotinidase activity is associated with a spectrum of deficiency ranging from total deficiency to heterozygous levels, a finding that is not always explained by the pathogenic variants observed in the BTD gene. The investigation of miRNAs, regulatory elements and variants in the 3’UTR region may present relevance in understanding the genotype-phenotype association. The aims of the study were to characterize the regulatory elements of the 3’UTR of the BTD gene and identify variants and miRNAs which may explain the discrepancies observed between genotype and biochemical phenotype. We evaluated 92 individuals with reduced biotinidase activity (level of heterozygotes = 33, borderline = 35, partial DB = 20 or total DB= 4) with previously determined BTD genotype. The 3’UTR of the BTD gene was Sanger sequenced. In silico analysis was performed to identify miRNAs and regulatory elements. No variants were found in the 3’UTR. We found 97 possible miRNAs associated with the BTD gene, 49 predicted miRNAs involved in the alanine, biotin, citrate and pyruvate metabolic pathways and 5 genes involved in biotin metabolism. Six AU-rich elements were found. Our data suggest variants in the 3'UTR of BTD do not explain the genotype–phenotype discrepancies found in Brazilian individuals with reduced biotinidase.

The BTD gene is composed by four exons, and its 3′UTR has 331 bp (ENST00000383778.5). The corresponding mRNA has two potential start codons (AUG) (Stanley et al., 2004, Pindolia et al., 2010. There are 17 different 3'UTR lengths with sizes ranging from 77 to 8226 pb, variable according to the transcript. The BTD gene has a constitutive expression pattern and healthy individuals present expression between 0.5 and 1.5 log10 transcripts per million (Figure 1). The threedimensional structure of biotinidase as predicted by in silico modeling consists of two domains (Pindolia et al., 2007).
Biotinidase deficiency (BD) is a metabolic disease, inherited in an autosomal recessive pattern, disabling the body to assimilate biotin from the diet and inhibiting biotin recycling (Baumgartner and Suormala, 2012). If not treated early, BD may lead to neurological and dermatological disorders (Wastell et al., 1988). BD may be total (activity <10% of normal) or partial (10-30%). There is an association between certain genotypes and the observed biochemical phenotype (total or partial), but in some patients, genotype and phenotype are mismatched. According to previous studies by our group, the association between the expected biochemical phenotype (according to genotype) and the actual biochemical phenotype occurs in 68.5% of cases, and variants in the 5′UTR of BTD do not seem to explain the variations found (Borsatto et al., 2014(Borsatto et al., , 2017(Borsatto et al., , 2019. Low activity of carboxylases can be found in BD and in Multiple Carboxylase Deficiency, a different disease caused by biallelic pathogenic variations in the HLCS gene, which encode the holocarboxylase synthetase enzyme (EC 6.3.4.10).
The aim of this study was to characterize the 3′UTR of the BTD gene in individuals with reduced biotinidase activity previously described by our group (Borsatto et al., 2014(Borsatto et al., , 2017(Borsatto et al., , 2019, and to identify which regulatory elements could influence the expression of biotinidase.

Material and Methods
The study was approved by the Research Ethics Committee of Hospital de Clínicas de Porto Alegre (n° 16-0480 and 12-0186), Brazil, and the subjects consented to participate by signing the Informed Consent Form.  Ninety-two individuals with reduced biotinidase activity were included: 33 with heterozygous level; 19 with borderline partial/heterozygous; 16 with borderline heterozygous/normal; 20 with partial deficiency; and 4 with total deficiency. These patients had the exons, exon-intron junctions, and 5′UTR of BTD previously sequenced, and were described by Borsatto et al. (2014) and Borsatto et al. (2017). The genotype and biochemical profile of the cohort is shown in Table 1, and details regarding the classification of the biochemical phenotype and BTD sequencing can be found in Borsatto et al. (2014) and Borsatto et al. (2017). Eighteen individuals had an inconsistent genotype-biochemical phenotype association (1-6, 24-33, 86, 87 - Table 1).
For genomic DNA extraction, blood samples were collected in EDTA-containing tubes and processed using the Easy-DNA gDNA Purification kit (Thermo Fisher). The 3′UTR of the BTD gene was amplified by PCR with specific primers. The products were purified with 20% PEG 8000/2.5M NaCl Table 1 -Genetic and biochemical profile of patients with reduced biotinidase activity included in the characterization of the 3'UTR.
To evaluate conservation of the 3′UTR of BTD, sequence alignments between different species were performed in MEGA software (version 7.0.26, Kumar et al., 2016), using the ClustalW algorithm (version 2.1, Thompson et al., 1994). Evolutionarily conserved regions were mapped in ECR Browser (Ovcharenko et al., 2004). The chromosomal position provided in the Atlas of UTR Regulatory Activity (AURA) (Dassi et al., 2014) was used to locate the 3′UTR of BTD gene.
The miRanda (Enright et al., 2003), mirSVR (Betel et al., 2010), and microRNA.org (Betel et al., 2008) algorithms were used for analysis of miRNAs target sites associated with BTD. The cutoff points for this analysis were a binding free energy of -25 Kcal -proposed as more stable by Seffens and Digby (1999), and the search for evolutionarily conserved targets (mainly 8-mer), as suggested by Garcia et al. (2011).
For polyadenylation analysis, the constitutive site was characterized according to the reference sequence curated by NCBI. APADB (Müller et al., 2014), APASDB (You et al., 2015) databases and the PolyA_SVM (Structural Support Vector Machine) algorithm of the RegRNA package (v. 2.0, Chang et al., 2013) were used to quantify sites usage and polyadenylation signals.
To identify other regulatory elements in the 3′UTR, the software RegRNA v. 2.0  and ARE Site 2 (Fallmann et al., 2015) were used. Secondary structures formed by miRNA-3′UTR interactions were obtained through the RNAfold Web server (Gruber et al., 2008).

Results
No variant was identified in the analysis of the 3′UTR of the BTD gene.
In the search of variants in genomic public databases, 43 variants were found in the AbraOM, 32 of them predicted as 'variant of uncertain significance' (VUS), and 11 as 'benign'. In the gnomAD, nine variants were found, all predicted by the ACMG as 'VUS'. In the LOVD database, three variants were found -one predicted as 'VUS' and two as 'benign'. The allele frequencies and the respective rsSNP as shown in Table 2.
In silico analysis of miRNAs yielded highly variable results. The number of miRNAs predicted in BTD gene were: 51 in miRGate database, 35 in miRTarBase, 5 in miRWalk, 4 in TarBase and 2 in TargetScanHuman (Table 3).
Seven miRNA target sites (Table 4) and one RNA binding protein (Musashi Binding Element) were identified. The mapped elements were presented in Figure 2. Forty-nine miRNAs were associated with genes that interact with the BTD gene in biotin metabolism (Table 5). The only miRNA shared between BTD and HLCS was the hsa-miR-222.
The three best-predicted secondary structure models are presented in Figure 3. The most appropriate secondary structure according to RNAfold analysis was the model of interaction between the 3′UTR of the BTD gene and hsa-miR-3934, with a binding free energy of -25.35 Kcal.
The polyadenylation signal used by the BTD gene coincides with the canonical AAUAAA hexamer. The dinucleotide that identifies the cleavage site was AA. Results from the APASdb database and the PolyA_SVM algorithm showed that the BTD gene has two major mapped polyadenylation sites. The first signal begins at position 2044 and has 32 pb; the second signal begins at position 2329 and has 32 pb. According to the APADB database, both polyadenylation sites of the BTD gene are located in the 3′UTR at positions chr3:15687323 (86.1% of usage) and chr3:15683749 (11.4% of usage).

Discussion
In this study, we investigated the presence of variants in the 3′UTR of the BTD gene in individuals with reduced biotinidase activity and, using bioinformatics tools, we discussed a possible relationship by regulatory elements with the expression of the BTD gene.
As far as we know, the 3′UTR had never been characterized in patients. As observed in the present cohort about 20% of the patients have discrepancies between expected BD according to genotype and type of BD according to enzyme activity.
The hypothesis for this investigation came from other diseases that present phenotype modification due to variants in the 3'UTR of the affected gene, as Glycogen Storage Disease Ia (Karthi et al., 2017) and Haemophilia A (Rosset et al., 2016). Modified regulatory elements may affect the interaction of the UTRs with proteins and microRNAs causing modulation of mRNA transcription, secondary structure, stability, localization, translation, and access to regulators like microRNAs (miRNAs), RNA-binding proteins (RBPs) and justify the discrepancies between genotype and phenotype (Steri et al., 2018;Skarp et al., 2020).
The high conservation of the 3'UTR was observed among the 92 patients analyzed proved by the 100% homology -no variants were found. Variant databases reinforced the conservation of the region through low frequencies of variants.
Subsequent investigations of the 3'UTR found several miRNAs and elements present in the region. Variations present in patients could justify differences in gene expression through factors related to 3'UTR.
The main predicted miRNAs associated with the BTD gene were: hsa-miR-7-5p, previously implicated in suppression of cell proliferation, induction of apoptosis, and angiogenesis (Li et al., 2016;Luo et al., 2018); hsa-miR-34a-5p, which is involved in cell proliferation and an important regulator of the central nervous system (Agostini and Knight, 2014;Jauhari and Yadav, 2019); and hsa-miR-145, identified in neonates and expressed specifically in the liver, where biotinidase expression is also higher (Fu et al., 2005;Noh et al., 2013).
The miR-7 cluster is known to be associated with genes related to the nervous system. Dostie et al. (2003) demonstrated that this miRNA may be unregulated in neuronal cells in spinal muscular atrophy, and involved in the neurological dysfunctions associated with Waisman Syndrome and Fragile X Syndrome. Untreated BD may lead to neurological problems and developmental delay. Thus, it is important to note that this miRNA, along with several potentially related factors, may be a candidate for investigation.
Hearing loss is a common sensorineural impairment in general populations. Experiments done in the inner ear of mice and humans have found differential expression of five miRNAs, among them miR-30, associated with different stages of ear development (Rudnicki and Avraham, 2012).
In the present analysis, miR-30 was associated with the BTD gene. Among patients with total BD, 75% of affected children have hearing loss (Wolf et al., 2002), with variable but usually irreversible severity.
Forty-nine miRNAs associated with genes that interact with the BTD were identified in the biotin metabolic pathway. These miRNAs have already been implicated in cell signaling, glycosylation pathways, and in arginine, biotin, tyrosine, and thiamine metabolism (Ortega et al., 2010). The PC gene that encodes pyruvate carboxylase, a biotin-dependent carboxylase, was found not only in the biotin metabolic pathway but also in alanine, citrate cycle, and pyruvate metabolic pathways (Rottiers and Näär, 2012).
Gene ontology analysis showed that these genes are involved in several biological processes, and act as coenzymes and in the metabolism of small molecules (Gene ontology: Fisher's exact with FDR multiple test correction: 9.95e-20 / 1.55e-15) (Thomas et al., 2003;Mi et al., 2013).
Among the most prominent results is the HLCS target gene. HLCS encodes the holocarboxylase synthetase that activates biotin-dependent carboxylases and catalyzes the binding of biotin to biotinidase. Experiments have shown that miR-539 decreases holocarboxylase synthetase levels, with the abundance of miR-539 being significantly higher at physiological biotin concentrations than in biotin-deficient and biotin-supplemented media, in all cell lines tested (Segura et al., 2013). The results of this study suggest that miR-539 may be one of several factors that detect biotin and regulate holocarboxylase synthetase levels. In the present study, this miRNA was not directly associated to the BTD gene, but to the holocarboxylase synthetase gene HLCS.
The SPCS1 and SPCS3 genes -subunits of the peptidase signal complex that act as hydrolases and participate in degradation of lysine (Kailes and Hartmann, 1996) -also stood out. The lysine present in the biotin-lysine complex (biocytin) is believed to be degraded through the action of this complex. The miRNAs associated with these genes may have an impact on expression of SPCS1 and SPCS3 and, consequently, on lysine degradation, preventing biotin recycled into its free form. In addition, hsa-miR-204 and hsa-miR-211, both predicted to be associated with SPCS1, are implicated in mechanisms of cell proliferation and metastasis in several types of cancer, including breast, colon, and lung cancer (Mazar et al., 2010). Esau et al. (2006) found that miR-122 allows the liver to function properly in adult mice. This miRNA is an important mechanism for regulation of genes involved in hepatic lipid metabolism. This corroborates the findings of Saha and Ruderman (2003) that observed negative effects on mice lipogenesis whereby a reduction in ACC gene expression, particularly ACC2, led to a decrease in malonyl CoA and subsequent increase in fatty acid oxidation. As biotinidase acts as a cofactor for several carboxylases, miRNAs may be involved in feedback regulation of this system. This miRNA was not found to be associated with BTD, but appears to be involved with citrate and pyruvate metabolism genes.
Based on the assumption that a single miRNA can regulate several target genes, miR-31-3p and miR-34a-5p were associated with the BTD gene and with the PCCA and PCCB genes, which encode subunits of the enzyme propionyl-coA-carboxylase, one of the biotin-dependent carboxylases. Dysfunction in these genes can lead to propionic acidemia, a disease characterized mainly by neurological and cardiac damage. Rivera- Barahona et al. (2017) found that these miRNAs are deregulated in the liver of mice; more specifically, overexpression of the miR-34 family is observed in patients with cardiac involvement, and is associated with other neurodegenerative diseases.

Conclusions
The present study was pioneer in the analysis of the 3′UTR of BTD gene in individuals with reduced biotinidase activity. Although the sequencing of this region has not found variants, it described their evolutive conservation.
The study of the 3′UTR in individuals with reduced biotinidase activity allowed us to conclude that variants in this region do not explain the genotype-phenotype discrepancies found in Brazilian patients. However, several factors as miRNAs sites and regulatory elements have been identified, which may influence the expression patterns of the BTD gene. To date, there are no strongly validated interactions between miRNAs and the BTD gene. Thus, its experimental validation remains as a perspective for future research.