Mitochondrial DNA Mutation Analysis in Breast Cancer: Shifting From Germline Heteroplasmy Toward Homoplasmy in Tumors

Studies have suggested a potential role of somatic mitochondrial mutations in cancer development. To analyze the landscape of somatic mitochondrial mutation in breast cancer and to determine whether mitochondrial DNA (mtDNA) mutational burden is correlated with overall survival (OS), we sequenced whole mtDNA from 92 matched-paired primary breast tumors and peripheral blood. A total of 324 germline variants and 173 somatic mutations were found in the tumors. The most common germline allele was 663G (12S), showing lower heteroplasmy levels in peripheral blood lymphocytes than in their matched tumors, even reaching homoplasmic status in several cases. The heteroplasmy load was higher in tumors than in their paired normal tissues. Somatic mtDNA mutations were found in 73.9% of breast tumors; 59% of these mutations were located in the coding region (66.7% non-synonymous and 33.3% synonymous). Although the CO1 gene presented the highest number of mutations, tRNA genes (T,C, and W), rRNA 12S, and CO1 and ATP6 exhibited the highest mutation rates. No specific mtDNA mutational profile was associated with molecular subtypes of breast cancer, and we found no correlation between mtDNA mutational burden and OS. Future investigations will provide insight into the molecular mechanisms through which mtDNA mutations and heteroplasmy shifting contribute to breast cancer development.


INTRODUCTION
Breast cancer is the most common cancer type in women worldwide (1); in 2018, 2.09 million cases and 627,000 deaths were estimated around the world (2). In Mexico, this malignancy is currently the leading cause of cancer-related deaths in women (3)(4)(5). Breast cancer is a heterogeneous disease in terms of the molecular features and clinical outcome; knowledge of the genetic causes is still incomplete. To comprehensively address this issue, the analysis of nuclear genes has been frequently done (4,(6)(7)(8)(9)(10). However, mutations in mitochondrial genes have also been reported to play important roles in human cancer development (11). Experimental evidence in mouse models demonstrated that G13997A and 13885insC mutations in the reduced form of the nicotinamide adenine dinucleotide dehydrogenase subunit 6 (ND6) gene are associated with the transition of poor metastatic tumor cells into a highly metastatic phenotype (12).
Mitochondria are indispensable for numerous cellular processes, such as respiratory energy metabolism, cell proliferation, apoptosis, senescence, and immune response (13). Mitochondrial DNA (mtDNA) is a circular DNA molecule of 16,569 bp containing 2 rRNA, 22 tRNA, and 13 genes encoding protein subunits for the mitochondrial complexes of the oxidative phosphorylation system (OXPHOS). The mitochondrial genome is characterized by many DNA molecules per cell and a high mutation rate that creates a state known as heteroplasmy (a mixture of mutant and normal mtDNA in cells) (14,15). The relative abundance of certain de novo mutations might be pathogenic, although the heteroplasmic threshold effect and the molecular mechanisms underlying the mutation selection during disease development have not yet been established (15). However, the association among several mtDNA germline variants and a broad spectrum of human malignancies has been widely documented (16)(17)(18). For instance, the T16189C (D-Loop) and G10398A (A114T at ND3 gene) germline variants have been commonly associated with an increased risk of cancer (17,(19)(20)(21). The important role of germline mtDNA variants in tumors has been noted in groups of variants inherited without any recombination, called haplogroups, conferring risk for cancer and as modifiers for potential metastasis and for the response to treatment in affected patients (22).
Genomic evidence suggests that mtDNA germline variants and tumor mutations are also involved in breast cancer development. mtDNA mutations have been reported in almost 60% of breast tumors, varying in nature (mainly single base substitutions and indels) and heteroplasmy levels (23)(24)(25)(26); however, whether these mutations are drivers or passengers remains largely uncharacterized (27,28). mtDNA sequencing from matched normal-tumor breast tissues has shown that a high proportion of the somatic mutations are singletons arising in a single patient and that many of those were under clonal expansion in the primary or metastatic tumors (17,22,(28)(29)(30)(31)(32)(33). To explore the landscape of mtDNA mutations in breast cancer from Mexican women and to know if mtDNA mutations and haplogroups are associated with the molecular subtypes of breast tumors or with overall survival (OS), we performed whole mtDNA sequencing in matched-paired primary breast tumors and peripheral blood.

Patient Data
We included matched tumor samples and peripheral blood from 92 women with breast cancer. Patients were recruited from the Instituto de Enfermedades de la Mama FUCAM (Mexico City, Mexico), and all patients provided written informed consent. Patients with chemotherapy, radiotherapy, and neoadjuvant treatment before the surgery were excluded. Based on the current international criteria and using immunohistochemical (IHC) markers, breast tumors were classified as luminal A (LA), luminal B (LB), HER2-positive (H2), and triple negative (TN) (34). Patients were followed up for at least 30 months (30-124) after the clinical diagnosis.

Sample Processing and DNA Extraction
Peripheral blood was obtained before surgical treatment using EDTA vacutainer tubes, and the buffy coat (leukocytes) was separated out by centrifugation at 1850 × g. Fresh tumor tissues were obtained during surgical treatment and stored at −80 • C until further analysis. Breast tumor samples with >80% of tumor cells were included. Total DNA from leukocytes and tumor samples were extracted using a Maxi Kit (Qiagen, Hilden, Germany) and AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany), respectively, according to the manufacturer's instructions. The nucleic acids were quantified by spectrophotometry (NanoDrop, Thermo Fisher Scientific, Waltham, MA, United States) and stored at −20 • C until library preparation.

Whole mtDNA Sequencing
To avoid contamination with nuclear mitochondrial segments, the whole mtDNA was enriched using long PCR with two overlapping primer sets (Supplementary Table 1

Bioinformatic Analysis
We included mitochondrial genome data with depth coverage with at least 300X per base. A quality control (QC) evaluation was performed on the raw sequence data using the FastQC algorithm (35), and reads with Phred quality >25 were accepted for further analysis. The Burrows-Wheeler alignment (BWA-MEM) algorithm was used to align all sequences against the revised Cambridge Reference Sequence (GRCh38, GCF_000001405.38) (36). Duplicated reads on mtDNA were removed. Then, reads were realigned and quality scores were recalibrated using Genome Analysis Toolkit (GATK) (37). Haplogroups were assigned on the mtDNA variants detected in blood samples and using Haplogrep software (v2. 1.19) and Phylotree databases (38,39). The Strelka2 algorithm with default parameters (40) was used to identify somatic mutations by comparing germline and matched mtDNA variants. Only variants labeled PASS were considered high-quality variants (40). To reduce false positives, a somatic mutation was considered when the mutant allele fraction (MAF) was <0.001 in normal tissue and >0.01 in tumor tissue in all samples where it occurred. Germline variants were defined if MAF > 0.01 in the blood sample regardless of the proportion in tumors and if this occurred in at least one patient. Variants not fulfilling these criteria were considered as probable mutations. These thresholds were estimated based on the depth coverage. All variants displaying mutant allele abundance from 1% to <95% were deemed heteroplasmic. Homoplasmic was defined as a mutant allele >95% (32). Mutations were reviewed in the MITOMAP database to determine if they had been reported in breast cancer or in other malignancies (41). To predict their role as a driver or passenger, we used the Variant Effect Predictor (VEP) software with the default parameters (42). To assess whether an mtDNA profile was associated with the tumor subtype, we performed an unsupervised hierarchical clustering analysis accounting for all mtDNA mutations using R (version 3.6.0). The distribution of mtDNA somatic mutations among breast cancer subtypes was visualized using Venn diagrams (43).

Statistical Analysis
Mutations were counted according to gene distribution, and frequencies were obtained for all analyzed patients. To assess the association between mtDNA mutations and haplogroups and clinical and molecular characteristics among groups, chisquare and Fisher's exact tests were calculated when appropriate. Mann-Whitney U tests were used to evaluate differences in heteroplasmy levels between the blood and tumor tissues. An OS analysis was performed, including patients who were followed up to 124 months after diagnosis and initial treatment, defining OS as the time between initial diagnosis and death. mtDNA mutational burden was stratified into two groups: patients without any mutation and patients with mutations, as well as according to the mutational burden mean (patients with a low or high mutational burden). To evaluate whether somatic mtDNA mutational burden is associated with patient outcome, Kaplan-Meier survival curves were obtained and the hazard risk was estimated using a Cox proportional hazard model adjusting for variables that have been reported to influence the OS (age, tumor stage, clinical stage, and hormone receptor status). P-values < 0.05 were considered statistically significant. Statistical analysis was conducted using SPSS (Version 25, IBM Corp., Armonk, NY, United States). Plots were generated using R (version 3.6.0).

Study Population
We included 92 unrelated patients with a mean age of 53.83 ± 11.5 years (range 33 to 92 years; Supplementary Table 2).

Distribution of Haplogroups
From 564 germline variants identified in the normal tissues, 296 were used by Haplogrep for haplogroup typing. After haplotyping the mtDNA of all patients, haplogroup A was the most frequent (44.6%), followed by haplogroups B (22.8%), C (11.9%), D (12%), and L (5.4%). Haplogroups H and J were found in 2.2 and 1.1% of patients, respectively. Stratification analysis by IHC subtypes revealed that haplogroups C and D were present only in LA and LB subtypes, and TN tumors carried haplogroups A and B. No statistically significant differences in the haplogroup distribution was found among the tumor subtypes (P > 0.05).  Figure 2). None of the 709 variants were detected in tRNA-A, tRNA-D, or tRNA-S1. Transitions and transversions accounted for 67% (459) and 33% (226) of the SNVs, respectively ( Figure 1A).

Mitochondrial DNA Variant Distribution
Of the total variants detected in tumor tissues, 324 were defined as germline variants and 173 as somatic mutations. Two hundred eleven variants were considered as probable mutations.

Distribution of Somatic Mutations
We detected 384 variants displaying MAF < 0.01 in blood samples, but >0.01 in tumor samples. The 211 variants showing MAFs between 0.001 and 0.01 in normal tissues were named as probable mutations, and the remaining 173 were identified as somatic mutations (<0.001 MAF in normal tissue). Although 14 (6.6%) of the probable mutations presenting very low MAF in normal tissue (<0.008) were homoplasmic in tumors and 47 (22.3%) of all probable mutations were detected in more than two breast tumors (being T8430A: 9.8% and A8439C: 6.5% the most frequent), toward reducing biases derived from the technology and bioinformatic tools used, our further analysis was focused only on those variants classified as somatic mutations.
Regarding the heteroplasmy status of the somatic mutations, 37 (21.4%) were homoplasmic, whereas the remaining 136 (78.6%) showed a high range of heteroplasmy levels. The T15115C, C16270A, and C16527T mutations, which were detected in two patients, exhibited a homoplasmic status in one patient but heteroplasmic in the other one.
From all somatic mutations, 56 (32.4%) are not recorded in the MITOMAP database (41), including two of the most frequent mutations detected in our cohort (G3219T and T3631A).

Functional Effect Prediction Analysis of Somatic Mutation
To determine whether the somatic mutations are either driver or passenger mutations, we annotated them using VEP software (42). Of 102 coding mutations, 34 (33.3%) were synonymous  Table 1). Given the nature of the mtDNA molecules, the T5726C and G5777A mutations (located in tRNA-N and tRNA-C genes, respectively) overlap in L-strand replication origin 1 (MT-OLR1) sequences and have a double function (Supplementary Table 3). Functional prediction analysis suggested that non-coding mutations (regulatory regions) might modify sequences associated with mtDNA replication and transcription processes ( Table 2).

Somatic Mutations Associated With Breast Cancer and Other Cancer Types
To determine if the mtDNA somatic mutations identified in this study have previously been associated with cancer, all mutations were examined in the MITOMAP database. We found that 25 (14.6%) of these mutations have been identified in several types of human cancer. Only G13333A and A15836C mutations, which were detected in one patient (1.1%), have been recorded in breast tumors ( Table 2).

Distribution of Somatic Mutations by Breast Tumor Subtype
We assessed the association between mitochondrial mutations and clinical characteristics of breast cancer. No significant differences were found in terms of age, status, or tumor subtype among patients carrying any number of mutations and with no mutations (P > 0.05).
No statistical significance was found in terms of the distribution of somatic mutations and breast cancer subtypes (P = 0.86; Supplementary Table 4).

Tumor Mutational Burden Association and Overall Survival
To determine the use of mtDNA tumor mutational burden as a potential predictor of survival in breast cancer, a Kaplan-Meier analysis was performed (Supplementary Figure 6). Regarding OS, our analysis revealed no statistically significant differences among positive versus negative mtDNA somatic mutation groups (HR = 1.637; 95% CI = 0.401-6.683), or when comparing a low (≤2 mutations) or high (>2 mutations) mutational burden (HR = 1.385; 95% CI = 0.329-5.833). The threshold was established according to the mutational burden mean (n = 2) of somatic mutations.

DISCUSSION
Breast cancer, the leading cancer type in women around the world, is characterized by its complexity and its highly heterogeneous genetic background, resulting from the interaction between known environmental factors and mutations in the nuclear genome and mtDNA. The mitochondrial genome function is relevant for the transformation processes; tumor cells are able to escape from different regulatory mechanisms, but they cannot evade energy flow mechanisms. Thus, cancer cells could gain a functional mtDNA mutation as a strategy to adjust energy metabolism during adaptation to oncogenic conditions (44)(45)(46). To evaluate the landscape of somatic mutation in breast tumors and to determine whether the mtDNA mutational burden is correlated with the OS of breast cancer, we sequenced 92 whole mtDNA paired peripheral blood-tumor samples from Mexican women.

Native American Haplogroups in Mexican Women With Breast Tumors
Seven different haplogroups were identified in our samples. The Native American haplogroups (A, B, C, and D) were more common than the major European (H and J) and African (L) haplogroups. Although the present study is limited by the selected population (women with breast cancer), we found similar frequencies of the A, B, and C haplogroups with those reported in the Mexican-Mestizo (mixed) populations and Mexican-Americans living in the United States (47,48). Considering our findings and those previously published, the Native American haplogroups' frequencies agree with the demographic history of the Mexican-Mestizo population that is composed of Native-American (53%), European (42%), and African (5%) ancestry (48,49). We found a higher frequency of the D haplogroup (12%) than in previously published works about the Mexican-Mestizo population (6.2%) (50,51). The D haplogroup has been associated with predisposition to breast cancer in the Chinese population and experimental data revealed that the D5 branch haplogroup promotes tumorigenesis through AKT activation, mediated by a high concentration of reactive oxygen species (ROS) (52,53). However, we need to include more patients and to perform a case-control analysis to determine whether the haplogroup D has clinical significance in Mexican women with breast cancer.

High Heteroplasmy Levels of Germline Variants in Tumors and Heteroplasmy Shifting
In this study, we were able to distinguish acquired mutations from germline variants and to detect heteroplasmy shifting by analyzing normal-tumor matched tissues. On the one hand, <33% of the germline variants exhibiting low heteroplasmy levels (1-2%) were detected in more than two subjects, and 29.4% of singletons germline variants have been annotated in the public databases, suggesting that these could be rare variants. On the other hand, for several germline variants (as A1736C, T4248C, and C16290T), the mutant allele was enriched in tumors in comparison with normal tissues (Figure 2). A663G, a marker of haplogroup A, was the most common germline variant. No evidence of an association among haplogroup A, or this variant with breast cancer risk, has been documented. However, the 663G allele has been reported as associated with various mitochondrial diseases (54)(55)(56)(57). Additionally, A663G has been found as a rare variant in papillary thyroid carcinoma and nasopharyngeal carcinoma (frequency <1%) (58,59). In silico predictions suggest that the change of A > G could destabilize the structure of the rRNA by interrupting the hydrogen bond between U and A, then altering its function (54,60). Nevertheless, its role in cancer development needs to be addressed.
Differences of heteroplasmy levels among tissues (blood <20%, tumors 1-95%) observed in this study could be biased by contamination by a different proportion of germline mtDNA during sequencing (69) or that the unregistered singletons germline variants showing <2% of heteroplasmy load arose during the genomic sequencing process (70). Our current strategy does not allow us to identify whether the heteroplasmy comes from the enrichment of mutated mitochondria or from mtDNA copy numbers (mtDNA-CN) variations. Diverse analyses have highlighted the great heterogeneity of the mtDNA-CN in different tissues and cancer types (71)(72)(73)(74), but studies comparing mtDNA-CN among blood and normal tissue adjacent to the tumor are scarce and have found inconclusive results (28,(74)(75)(76)(77).
mtDNA heteroplasmy is a normal condition of randomized processes and regulated homeostatic mechanisms, such as the mitochondrial segregation during the cell division, and mitochondrial fusion and fission process (78,79). The changes in allelic frequencies have been suggested to be the response to selective pressures generated under physiological and pathological environmental conditions (17). However, it could be also as a result of changes in mtDNA-CN that occurred during the mtDNA replication (which is cell cycle independent) (80)(81)(82), or under oxidative stress induced by endogenous and exogenous factors (hormones, age, dietary, etc.) (83). The correlation among heteroplasmy and mtDNA-CN has scarcely been explored, and higher levels of heteroplasmy were correlated with lower mtDNA-CN in the central nervous system, but not in other tissues (73). In breast cancer, controversial results are found (83,84), and it has been suggested that both processes are not interacting with each other and play different roles in the development of breast tumors (83). Although we have to consider that the shift toward homoplasmy observed in the present work could be influenced by mtDNA-CN rather than an enrichment of mutated mitochondria, heteroplasmy shifting has been reported in kidney and thyroid carcinomas (85), as well as breast cancer (28). Studies indicate that tumor cells show a neutral evolution of their mtDNA and even have the ability to tolerate pathogenic mutations; however, some mutant alleles reaching a critical threshold can contribute to cancer development and progression (85). The molecular mechanisms by which these alleles (whose heteroplasmy levels shift in the tumor) contribute to carcinogenesis is still unclear.

High Frequency of Heteroplasmic mtDNA Mutations in Breast Tumors
In our study, 24 (26.1%) breast tumors were negative for somatic mutations and 68 (73.9%) were positive; over 78% of these mutations were heteroplasmic. Studies of breast cancer and diverse types of tumors have reported different proportions of tumors carrying mtDNA mutations (73.7 and 45.7%, respectively) (26,28,69,86). The discrepancies among these findings could be partially explained by limitations derived from the sequencing methods, or due to the different thresholds used to define a somatic mtDNA mutation (28,86). We found a high frequency of heteroplasmic mtDNA mutations and coding mutations (mainly non-synonym mutations), as well as a high mutation rate in the D-loop region and tRNA genes, as previously reported in several types of tumors (86,87).
Coding mtDNA mutations could alter the protein function and OXPHOS system, and non-coding mutations might affect fundamental biological processes such as mtDNA replication, transcription, and structural mtDNA organization, thus conferring advantages for tumor cell proliferation (88,89).
It has been suggested that deleterious pathogenic mutations, arising during the carcinogenic process are eliminated when protein function and viability cells are compromised, but selected if they enable cell proliferation (44,(90)(91)(92)(93). In fact, it has been proposed that due to the ND5 proton channel function, mutations in this gene might contribute to a deficient yield of the OXPHO system, leading to a high rate of oncogenic mutations in nuclear DNA and mtDNA (94,95). Studies describing mutations in metastatic lesions, which were not found in primary tumors, support the hypothesis of the relevance of mtDNA mutations in cancer progression (90,96). Li et al., by analyzing 26 types of cancer, proposed that regions and functional units are under negative and positive selection processes, suggesting the important role of mtDNA mutations in tumors (86).
Various authors stated that missense mutations of mtDNA are non-tumorigenic since they do not determine the impairment of the respiratory chain (28,97). Likewise, since tRNA mtDNA pathogenic mutations are rarely fixed, the deleterious effect of these mutations might be functionally compensated for by gaining new mutations or altering mtDNA-CN (98). Other authors have proposed random processes defining the landscape of mtDNA mutations in cancer, suggesting that the heteroplasmic status arises in tumor progenitor cells by chance, without conferring any physiological or tumorigenic advantages (87). Primary tumors may have multiple clonal subpopulations that emerged during the numerous selective events required to complete the tumorigenesis process (99); thus, we have to consider that in addition to mtDNA-CN changes, clonal selection occurring during breast cancer progression could influence the levels of heteroplasmy observed in the present study and in other types of tumor (80,100). More experiments are needed to understand the mechanisms underlying the mtDNA germline variants selection toward a homoplasmic state in breast cancer development.
As several of the recurrent mutations identified in this analysis (carried by >3% of our patients) have been related to other human diseases (58,62,101,102), these mutations could also play a role in breast tumor development. However, it is necessary to validate our findings in an external cohort, including a population with a different haplogroup background. As well, we have to take into consideration that driver and passenger mutations may acquire oncogenic potential when, next to them, mtDNA gains more variants. In animal models, it was observed that the mtDNA containing two mutations in the ND6 gene (G13997A and 13588insC) enhances the metastatic potential of tumor cells (12). Functional experiments could improve our understanding of the oncogenic biological processes mediated by these mutations (103,104). The use of transmitochondrial hybrids might be helpful to address these issues (12).
The differentiation between a somatic mutation and a germline variant is a very complex task during the study of the mitochondrial genome. Technical and bioinformatics limitations might derive in misclassification of the variants. An interesting observation was the detection of the two coding variants T8430A (L22H) and A8439C (Q25P) as part of the probable mutation group. 8430A allele has been proposed to be associated with elevated ROS concentrations and reduced ATP synthesis (67, 105, 106). The 8439C allele has been described as a mutation in breast cancer and in silico analysis suggested that this mutation is highly deleterious and affects the ATP8 function (107). Greater depth coverage of the mtDNA analysis could precisely determine the nature of the probable mutations (germline variants or somatic mutation) found in this work. As well, additional studies on T8430A and A8439C are needed to decipher their role in breast cancer development and progression.

mtDNA Mutations and Clinical Relevance
In addition to the description of the landscape of the mtDNA mutations in breast tumors, we also evaluated their clinical relevance in tumor classification and prognosis. Over 85% of the somatic mutations were found with low frequency (<1%), probably due to individual variations (28). The LA and LB tumor subtypes displayed the highest number of mutations, sharing 4.1% of the mutations. However, we found neither an mtDNA mutation profile nor haplogroups associated with tumor subtypes. Few studies have explored the association between mtDNA mutational burden and haplogroups with breast tumor molecular subtypes, and their results are still inconclusive (51). Based on nuclear DNA mutation data showing that hormonereceptor-negative tumors have a higher mutational burden than hormone-receptor-positive tumors (108), we expected to observe the same behavior regarding mtDNA mutations. We suggest that the elevated mutational burden of hormonereceptor-positive tumors (LA and LB) could be explained by superoxide radical formation during estrogen metabolism, which promotes mtDNA alterations (103,109). A recent study revealed a differential expression profile of mitochondria-related genes among molecular subtypes of breast cancer, perhaps derived from the partial contribution of the mtDNA mutations in breast tumor biology (110). Increasing sample size and balancing tumor subtypes could support our observations.
Overall survival analysis showed no differences among patients carrying a high versus low mtDNA mutational burden in tumors, which is in contrast with the data reported in breast cancer (28). A high mtDNA mutational burden has been associated with a worse outcome in pulmonary adenocarcinoma, acute myeloid leukemia, and pancreatic ductal adenocarcinoma (18,111,112); the opposite was observed in oral squamous cell carcinoma, and even acute myeloid leukemia, and breast cancer (28,113,114). We cannot exclude potential bias due to the nature of the studied mutations, the clinical characteristics of the patients included, or the ethnic/genetic background of the populations studied. Our small sample size and potential population stratification, since all patients were recruited in Mexico City, might also bias our results. As well, we cannot ignore the fact that other alterations, such as the mtDNA-CN (suggested as a potential prognostic biomarker for breast cancer), mito epigenetic processes, and nuclear genes that are involved in the mitochondrial biogenesis, may be clinically relevant (18,71,115,116). Thus, studies are necessary for larger cohorts of patients to determine the significance of the mtDNA mutational burden as a biomarker in breast cancer.
This study describes the distribution of mtDNA germline variants and mutations in breast tumors in a population of Mexican women. We found a high mutation rate in the D-loop region and tRNA genes. Heteroplasmy analysis suggested that negative and positive selection processes are shaping the landscape of mtDNA mutations in breast cancer, but functional experiments are needed to further understand the oncogenic biological processes.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: European Nucleotide Archive (https://www.ebi.ac.uk/ena) (Accession: PRJEB40354).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Instituto Nacional de Medicina Genómica. The patients/participants provided their written informed consent to participate in this study.

FUNDING
This study was funded by grants from the Consejo Nacional de Ciencia y Tecnología (CONACyT) through the program Fondo Sectorial de Investigación en Salud y Seguridad Social (S0008-2016-1-272618) and the Instituto Nacional de Medicina Genómica (C2_19/2009 and 24/2016/I). CP-A was supported by CONACyT (CVU 821714). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
We thank the Unit of Sequencing (USeq) at INMEGEN for technical support and the Department of Computational Genomics at INMEGEN for the bioinformatic support. Finally, we are extremely grateful to all the patients who participated in the study.