Rare Variants in Autophagy and Non-Autophagy Genes in Late-Onset Pompe Disease: Suggestions of Their Disease-Modifying Role in Two Italian Families

Pompe disease is an autosomal recessive disorder caused by a deficiency in the enzyme acid alpha-glucosidase. The late-onset form of Pompe disease (LOPD) is characterized by a slowly progressing proximal muscle weakness, often involving respiratory muscles. In LOPD, the levels of GAA enzyme activity and the severity of the clinical pictures may be highly variable among individuals, even in those who harbour the same combination of GAA mutations. The result is an unpredictable genotype–phenotype correlation. The purpose of this study was to identify the genetic factors responsible for the progression, severity and drug response in LOPD. We report here on a detailed clinical, morphological and genetic study, including a whole exome sequencing (WES) analysis of 11 adult LOPD siblings belonging to two Italian families carrying compound heterozygous GAA mutations. We disclosed a heterogeneous pattern of myopathic impairment, associated, among others, with cardiac defects, intracranial vessels abnormality, osteoporosis, vitamin D deficiency, obesity and adverse response to enzyme replacement therapy (ERT). We identified deleterious variants in the genes involved in autophagy, immunity and bone metabolism, which contributed to the severity of the clinical symptoms observed in the LOPD patients. This study emphasizes the multisystem nature of LOPD and highlights the polygenic nature of the complex phenotype disclosed in these patients.


Introduction
Glycogen storage disease type II (GSDII; Online Mendelian Inheritance in Man (OMIM) 232300; Pompe disease or acid maltase deficiency) is an autosomal recessive disorder caused by mutations in the acid alpha-glucosidase (GAA) gene encoding a lysosome enzyme that hydrolyses glycogen to glucose [1].
Mutations in GAA lead to a wide spectrum of clinical phenotypes, ranging from a severe infantile form, presenting with cardiomyopathy and muscular hypotonia, to a relatively mild, LOPD form characterized by a slowly progressing proximal muscle weakness, often involving respiratory muscles [1,2]. More recently, many additional symptoms associated with LOPD came to light: dysarthria and dysphagia, osteoporosis, scoliosis, sleep apnoea, small fibre neuropathy, hearing loss, impaired gastric function, urinary tract and anal sphincter involvement, pain and fatigue, as well as a risk of cardiac arrhythmia and cerebral and intracranial aneurysms [3]. These clinical findings emphasize the multisystem nature of LOPD.
In LOPD patients, the levels of GAA enzyme activity and the severity of clinical pictures may be highly variable among individuals, even in those who harbour the same combination of GAA mutations [4,5]. This indicates that there is no definite correlation between enzyme activity levels and severity of phenotype in LOPD patients [5][6][7][8][9][10][11][12][13][14][15]. This evidence leads to the hypothesis that the expression of GAA mutations could be modified by genetic background and non-genetic factors (life-style, nutrition and environment) [16,17]. To identify the modifier factors, which might affect the severity of symptoms of LOPD patients, studies on large informative families are advisable. Indeed, GAA-mutated family members have the closest possible genetic background and might better define the genotype-phenotype correlation, including a wide spectrum of clinical parameters.
In this manuscript, we report on a detailed clinical, morphological and genetic study of 11 adult LOPD siblings belonging to two large Italian families carrying compound heterozygous GAA mutations [18]. To better manage the wide spectrum of clinical and biochemical data, we calculated a severity index that was useful to simplify the genotypephenotype relationships in both LOPD families.

Clinical Evaluation of the LOPD Patients
In this study, we focused our analysis on 11 LOPD siblings from two families (family 1: II-4, II-5, II-7, II-8, II-9, II-12, II-13; 4 females and 3 males; family 2: II-I, II-3, II-4, II-5; 2 females and 2 males) from Southern Italy (ages 50-67 years for family 1, and 41-51 years for family 2) ( Figure 1A,B). Family 1 was already published [18]. To reduce the genetic heterogeneity, we included in this study only the patients of the second generation who carried the same combination of GAA mutations. Subjects II-1 and II-6 had a different GAA genetic background [18], a less severe phenotype than other siblings did, and they refused both ERT and WES analysis. A complete description of all the clinical data collected in this study is reported in Table 1. All data were collected at baseline before ERT treatment. The immune response to ERT was assessed after the first two infusions and at three, six and twelve months of treatments, according to the ERT protocol, and the data shown in the Table 1 were those recorded at 12 months.
Mean age at first symptoms was 41.7 years ± 6.7 standard deviation (SD) (range 29-53 years), with a higher average age at onset for females (43.5 ± 8.5 years) versus males (39.6 ± 3.2 years). The mean time lag between early symptoms and diagnosis was 12.5 ± 4.5 years. Analysis of motor function revealed mild reduced muscle strength in all patients, primarily affecting paraspinal and pelvic muscles, with manual muscle testing (MMT) mean score of 83.54% ± 8.61 and Gait, Stairs, Gower, Chair (GSCG) score of 11.09 ± 6.6. The 6-min walking test (6MWT) revealed a reduced physical functional capacity, achieving a distance of 241.66 ± 119 m for females and 279.2 ± 122.2 m for males.
The severity of muscle weakness varied between the patients regardless of gender A complete description of all the clinical data collected in this study is reported in Table 1. All data were collected at baseline before ERT treatment. The immune response to ERT was assessed after the first two infusions and at three, six and twelve months of treatments, according to the ERT protocol, and the data shown in the Table 1 were those recorded at 12 months. Mean age at first symptoms was 41.7 years ± 6.7 standard deviation (SD) (range 29-53 years), with a higher average age at onset for females (43.5 ± 8.5 years) versus males (39.6 ± 3.2 years). The mean time lag between early symptoms and diagnosis was 12.5 ± 4.5 years. Analysis of motor function revealed mild reduced muscle strength in all patients, primarily affecting paraspinal and pelvic muscles, with manual muscle testing (MMT) mean score of 83.54% ± 8.61 and Gait, Stairs, Gower, Chair (GSCG) score of 11.09 ± 6.6. The 6-min walking test (6MWT) revealed a reduced physical functional capacity, achieving a distance of 241.66 ± 119 m for females and 279.2 ± 122.2 m for males.
The severity of muscle weakness varied between the patients regardless of gender and duration of the disease. A significant association between age at onset and muscle weakness was previously demonstrated in affected siblings of family 1 [18]. In accordance to this finding, looking at family 2 (Table 1), patients II-4 and II-5 showed worse motor capacity with earlier onset of symptoms (respectively at 38 and 29 years old) compared to other family members; in particular, subject II-5 was able to walk only with double support for a few meters (10 m at 6MWT). According to the Fatigue Severity Scale (FSS, Italian version), all our patients obtained a score higher than 5 (severe fatigue) [19]. The mean forced vital capacity (FVC) value in the upright position of all patients was 76.81% ± 21.76; the FVC was ≥80% of the predicted normal in 5 out of 7 siblings of family 1 8,9,12,13) and only in II-1 of family 2 ( Table 1). The other subjects analysed had reduced FVC upright values, which dropped in the supine position. However, patients II-4 and II-5 of family 2 (Table 1) were not able to perform spirometry in the supine position due the severe drop in respiratory performance, for which they needed non-invasive ventilation (NIV). These findings were suggestive of a restrictive respiratory pattern ranging from a mild to a moderately severe degree [20]. A meaningful correlation between respiratory impairment and age at onset was already reported in family 1 [18]; this is also true in family 2, in which the youngest members (II-4 and II-5) showed the most severe respiratory conditions compared to those of the other siblings (Table 1). Body mass index (BMI, kg/m 2 ) was assessed since it represents a significant comorbidity in determining reduced motor capacity and worsening respiratory function. Although the patients with metabolic myopathies usually have normal or reduced BMI due to progressive muscle wasting, our patients showed an increased BMI mean value (27.74 ± 5.57) (Table 1). Furthermore, numerous studies demonstrated a high prevalence of osteoporosis in patients with LOPD [21]. To explore this aspect, we measured in all patients the circulating level of vitamin D and performed the analysis of bone mass density (BMD). The vitamin D levels were reduced in all subjects. In particular, patients II-12 and II-13 of family 1 and II-2 of family 2 (Table 1) showed a mean value of serum vitamin 25(OH)D of 19.25 ± 13.11, which denoted a status of insufficiency/deficiency. BMD was suggestive of osteopenia in most patients with reduced levels of vitamin D, whereas subjects F1-II-8 and F2-II-3 had normal T-score values despite vitamin D deficiency. In the patient II-4 of family 1, BMD indicated osteoporosis in accordance with her personal history of spontaneous vertebral collapse of two dorsal vertebrae. Furthermore, bone fractures occurred with low trauma and at a young age, with delayed healing, in patients F1-II-8 at 30 years old (proximal radius fracture) and in F2-II-5 at 42 years old (neck femur fracture). Additionally, LOPD is a predisposing condition to dilative arteriopathy and cerebral aneurysm formation [22], while structural cardiac alterations do not seem to be present when Pompe disease begins in adulthood [23]. We observed from brain magnetic resonance imaging (MRI) dilated intracranial internal carotid and basilar artery dolichoectasia in all siblings of family 1. Vascular anomalies in patient II-8, who carried a paramagnetic prosthesis, were confirmed through Angio-TC. Echocardiography revealed mild mitral prolapse in all subjects of family 1, except for II-13. Members of family 2 did not show any cerebral-vascular anomaly or cardiac valvular diseases. At the time of first evaluation of the patients, the GAA activity on Dried Blood Spot (DBS) was assessed, revealing a significantly reduced activity (mean values 0.56 ± 0.78 mol/h/L). Surprisingly, in subject II-4 of family 1, the dosage was in the normal range. A high number of lymphocytes containing Periodic Acid-Schiff (PAS)-positive granules were present on blood smears in the majority of the patients, with mean percentage values of 50.36% ± 21.36% (Table 1). Histological analysis of the muscle samplings detected a percentage of vacuolated fibres of 3.09% ± 5.02%, higher in the patients II-9 of family 1 (11.18%) and II-5 of family 2 (15.31%). Although the vacuolated fibres are representative of autophagy impairment in skeletal muscles, they do not seem to correlate with duration of the disease, age at biopsy nor with the residual GAA enzyme activity [18,24]. Finally, response to ERT and adverse reaction to therapy were considered. Specific serum IgG against rh-GAA ranged 850 to 6400 after one-year therapy. Patient II-5 of family 1 developed a severe adverse reaction after 18 months of ERT, disclosing diffuse cutaneous rash and shortness of breath during infusion. There were increased levels of specific serum IgG (51,000 mg/dL) and C3 and C4 fractions of the complement, while the C1 inhibitor was decreased (0.1 g/L) and the IgE antibodies against rh-GAA were absent.

Generation of a Phenotype Severity Index (SI) for LOPD
To simplify the large amount of clinical results and define the relationship between clinical variability and genetic factors, we elaborated a phenotype SI. The minimum score for SI was 1, which indicates a diagnosis of LOPD without additional significant clinical signs, while the maximum SI score was 47. Results are reported in Table 2. The grade of severity for a single item was assigned: grade 0 indicated normal findings; age at onset (years old) was graded from 1 to 6 ( . BMI (Kg/m 2 ) score was determined according to the definition of a normal weight, overweight and obesity as previously described (grade 0: 19.5-24.9; grade 1: 25-29.9, grade 2: 30-35.9, grade 3: 36-40, grade 4: >40). Reduced vitamin D 25-OH serum levels (ng/mL) were classified in accordance to the common definitions of insufficiency and deficiency (grade 0: 30-100; grade 1: 20-30; grade 2: <20), and for the BMD were considered values of the T-score (g/cm 2 ) of three different districts if suggestive of osteopenia (grade 1) or osteoporosis (grade 2). We assigned a grade of 1 to bone fractures (grade 0: no history of bone fractures; grade 1: presence of bone fractures), basilar artery dolichoectasia (grade 0: normal findings at brain MRI; grade 1: presence of BAD) and to mitral valve prolapse (grade 0: normal mitral valve; grade 1: presence of mitral valve prolapsed). Equally was attributed a score of 1 to the relief of serum IgG anti-rhGAA (grade 0: absence of IgG; grade 1: presence of IgG).
For the FSS, all our patients obtained a score higher than 5 (severe fatigue) [19] and therefore the inclusion of this measure would not have changed the severity index calculation.
As reported in Table 3, the clinical variables significantly associated with a worse severity score were global 6MWT (p = 0.03), FVC (p = 0.01), ∆FVC (p = 0.02), VDD (p = 0.01), BMD femoral neck (p = 0.03) and GSGC (p = 0.01). Considering the small sample size, 6MWT was analysed as global because the sex distribution was well distributed between the two subgroups. The association between 6MWT, GSGC, FVC, ∆FVC and the severity index score was also confirmed by Pearson correlation.

Comprehensive Genetic Analysis of the LOPD Patients
At the molecular level, the two families carried compound heterozygous mutations in GAA. Specifically, the 7 siblings of family 1 carried the two mutations c.118C > T (p.R40X) (reported as very severe in the Pompe Mutation Database) and c.2647-7G > A (p.N882fs) (potentially mild) [18], while the 4 siblings of family 2 carried the two mutations c.-32-13T > G (potentially mild) and c.1124G > T p.(Arg375Leu) (potentially less severe). Although the GAA c.2647-7G > A (p.N882fs) splicing variant of family 1 was associated to variable expression of the mutated allele [18], this mutation alone cannot explain the complex phenotype observed in these patients.
Considering that affected siblings belonging to the same family share the same GAA genomic region, we used the WES approach to search for other gene variants/mutations that might influence the severity of the LOPD phenotype.
We adopted a prioritization scheme to identify novel variants in modifier genes, similar to the approach taken in our recent studies ( Figure 2) [27,28]. We annotated 101,756 high-quality variants from which we selected splicing variants (including splice acceptors and donor sites) and exonic variants, including non-synonymous variants (NSV), stop gain/loss and short coding insertions or deletions (indels; I) (16,230 variants); synonymous variants were excluded. Subsequently, we excluded the very common variants (minor allele frequency (MAF) ≥ 0.2) and selected the variants reported as highly deleterious (Combined Annotation Dependent Depletion (CADD) score ≥ 20) (2,806 variants), corresponding to those variants predicted to be amongst the 1% of the most deleterious variants in the genome. To dissect the different aspects of the disease manifestation, we performed dichotomization of the clinical signs as reported in the SI score (Table 2). Each patient was classified for each specific symptom as positive/negative or mild/severe. Subsequently, we performed two different analyses to identify the variants/genes shared between the two families (Analysis I, Figure 2) and the variants shared among relatives in each family (Analysis II, Figure 2). We selected 589 genes, including 174 shared between the two families, 160 specific to family 1 and 255 specific to family 2. To further reduce the number of variants, we performed a functional prioritization of the candidate genes using a STRING database analysis, as described in Material and Methods.
We identified 16 rare and deleterious variants (14 non-synonymous 1 non-frameshift insertion and 1 splicing) lying in 12 genes, involved across different interconnected pathways potentially deregulated in LOPD, including autophagy/lysosomal (RILP, FNIP2, TRAPPC11, PLIN2, IRS1, KL, LRP4, RUNX1), immunity (RILP, FNIP2, TRAPPC11, IRS1, KL, LRP4, RUNX1, FAM26F), bone metabolism (TRPV5, TRPV6, IRS1, KL, LRP4, RUNX1) and skeletal muscle development and/or disease (TRAPPC11, PLIN2, IRS1, LRP4, RUNX1, DCHS1) (Tables 4 and 5 and Figure 3). Three out of the 12 genes, Rab Interacting Lysosomal Protein (RILP), Folliculin Interacting Protein 2 (FNIP2) and Perilipin 2 (PLIN2), were shared between the two families. RILP and FNIP2 were associated to a SI score > 20 and PLIN2 was associated with obesity in the two families. Polymorphic variants in the DCHS1 gene, which causes dominant Mitral Valve Prolapse type 2 (MVP2) [26], were identified in all patients with MVP as single or compound heterozygous (patients II-5, II-7 and II-12), suggesting their role as MVP susceptibility factors. (Analysis II, Figure 2). We selected 589 genes, including 174 shared between the two families, 160 specific to family 1 and 255 specific to family 2. To further reduce the number of variants, we performed a functional prioritization of the candidate genes using a STRING database analysis, as described in Material and Methods. We identified 16 rare and deleterious variants (14 non-synonymous 1 non-frameshift insertion and 1 splicing) lying in 12 genes, involved across different interconnected pathways potentially deregulated in LOPD, including autophagy/lysosomal (RILP, FNIP2,   Severity score > 20  Figure 3). Three out of the 12 genes, Rab Interacting Lysosomal Protein (RILP), Folliculin Interacting Protein 2 (FNIP2) and Perilipin 2 (PLIN2), were shared between the two families. RILP and FNIP2 were associated to a SI score > 20 and PLIN2 was associated with obesity in the two families. Polymorphic variants in the DCHS1 gene, which causes dominant Mitral Valve Prolapse type 2 (MVP2) [26], were identified in all patients with MVP as single or compound heterozygous (patients II-5, II-7 and II-12), suggesting their role as MVP susceptibility factors.   In the patient II-5 of family 1, who developed adverse response to ERT, we identified a splicing variant c.525 + 2T > G into the intron 2 of FAM26F gene, which codifies for a recently identified tetraspanin-like membrane glycoprotein that is involved in various immune modulating responses. A 100% probability was calculated by bioinformatics analysis that this variant completely disrupted the canonical donor splicing site into intron 2 of In the patient II-5 of family 1, who developed adverse response to ERT, we identified a splicing variant c.525 + 2T > G into the intron 2 of FAM26F gene, which codifies for a recently identified tetraspanin-like membrane glycoprotein that is involved in various immune modulating responses. A 100% probability was calculated by bioinformatics analysis that this variant completely disrupted the canonical donor splicing site into intron 2 of the FAM26F gene, suggesting that it was a deleterious change generating a frame shift (p.S174fs) and a predicted truncated protein of 224 amino acids. A reduced expression of the wild-type transcript was disclosed in patient II-5 combined with the presence of the aberrant spliced FAM26F transcript ( Figure 5A,B). the FAM26F gene, suggesting that it was a deleterious change generating a frame shift (p.S174fs) and a predicted truncated protein of 224 amino acids. A reduced expression of the wild-type transcript was disclosed in patient II-5 combined with the presence of the aberrant spliced FAM26F transcript ( Figure 5A,B). Considering that the insufficiency/deficiency of vitamin D is the most important biochemical/clinical biomarker identified in our families, we investigated the allelic contribution in our patients of variants associated with low levels of vitamin D in genes involved in its metabolism. As reported in Table 6, we found in the two families an excess of CC alleles for the two variants rs2228570 (Vitamin D receptor (VDR) gene, c.T2C (p.M1T), also known as the FokI variant) and rs4588 (Vitamin D-binding protein (GC/DBP), c.C1307A (p.T436K)), which were reported as associated with low vitamin D levels in several populations [29].  Considering that the insufficiency/deficiency of vitamin D is the most important biochemical/clinical biomarker identified in our families, we investigated the allelic contribution in our patients of variants associated with low levels of vitamin D in genes involved in its metabolism. As reported in Table 6, we found in the two families an excess of CC alleles for the two variants rs2228570 (Vitamin D receptor (VDR) gene, c.T2C (p.M1T), also known as the FokI variant) and rs4588 (Vitamin D-binding protein (GC/DBP), c.C1307A (p.T436K)), which were reported as associated with low vitamin D levels in several populations [29].

Expression Studies
The immunofluorescence study on muscle biopsies showed that the GAA protein defined the profile of the muscle fibres and, in some of them, marked small sarcoplasmic spots, some of which were also p62 positive ( Figure 6). p62, also known as Sequestosome 1 (SQSTM1), is involved in autophagy-dependent elimination of ubiquitinated proteins, and autophagy inhibition leads to the accumulation of p62-positive aggregates. RILP was expressed in the subsarcolemmal spots and at the level of the sarcoplasmic micro vacuoles in both RILP mutated and non-mutated biopsies of LOPD patients; however, the reaction appeared less intense in mutated tissues. In healthy subjects (CNT1), a faint subsarcolemmal continuous signal was observed in the majority of fibres, while CNT2 showed almost no fluorescent signal (Figure 7). The PLIN2 signal was intense and continuous at sarcolemma in mutated and wild type muscle fibres, while it appeared faint and discontinuous in external controls (CNT1, CNT2). PLIN2-positive sarcoplasmic macro and micro vacuoles were also evident in the mutated muscle, while in both the wild type and external control tissues smaller and scattered sarcoplasmic PLIN2-positive spots were seen ( Figure 8).
The FNIP2 reaction intensely labelled the vacuoles and micro vacuoles that also expressed the lysosome marker LAMP2 in the mutated biopsies. In the wild type as well as in external control tissues, the FNIP2 reaction marked weakly and diffusely the sarcolemma and the sarcoplasm of some fibres (Figure 9). The c.G850A (p.G284S) variant in RILP was associated with the most severe SI score in the two families. IF analysis in HeLa cells showed that the RILP284G protein accumulated near the nucleus (a strong signal intensity was observed) and co-localized with the LAMP2 antibody, suggesting lysosome and autophagosome localization. In contrast, a low signal intensity was observed for the RILP284S protein variant, associated with the SI score > 20, which was spread all over the cytoplasm and did not show a complete co-localization with LAMP2 and p62 antibodies ( Figure 10). These data suggest that this variant might perturb the binding of the protein to lysosomes.        The c.G850A (p.G284S) variant in RILP was associated with the most severe SI score in the two families. IF analysis in HeLa cells showed that the RILP284G protein accumulated near the nucleus (a strong signal intensity was observed) and co-localized with the LAMP2 antibody, suggesting lysosome and autophagosome localization. In contrast, a low signal intensity was observed for the RILP284S protein variant, associated with the SI score > 20, which was spread all over the cytoplasm and did not show a complete colocalization with LAMP2 and p62 antibodies ( Figure 10). These data suggest that this variant might perturb the binding of the protein to lysosomes. The localization of the RILP wt is indicated by intense green spots localized near the nucleus. The LAMP2/CY3 reaction was quite superimposable to that of RILP. The RILP mutated protein is visualized as less intense green spots spread all over the cytoplasm. An incomplete co-localization with LAMP2 and p62 is highlighted.
In HeLa cells, the PLIN2 wt and mutated proteins (PLIN2206K and PLIN2309C) localized in the peripheral spots. The PLIN2 wt and PLIN2206K proteins also localized in the nucleus with different intensity (Figure 11). The localization of the RILP wt is indicated by intense green spots localized near the nucleus. The LAMP2/CY3 reaction was quite superimposable to that of RILP. The RILP Figure 10. Immuno-localization of the RILP wt (RILP-284G-GFP) and mutant (RILP-284S-GFP) proteins in the HeLa cells wt and mutant. RILP proteins were fused to the green fluorescent protein (GFP), and LAMP2 and p62 proteins were stained using anti-LAMP2/CY3 and anti-p62/CY3. F-actin fibres were stained with rhodamin phalloidin and nuclei were stained with Hoechst. Slides were visualized with a Nikon Confocal Microscope A1R at 60× magnification. mutated protein is visualized as less intense green spots spread all over the cytoplasm. An incomplete co-localization with LAMP2 and p62 is highlighted. In HeLa cells, the PLIN2 wt and mutated proteins (PLIN2206K and PLIN2309C) localized in the peripheral spots. The PLIN2 wt and PLIN2206K proteins also localized in the nucleus with different intensity (Figure 11). Figure 11. Immuno-localization of the PLIN2 wt (PLIN2-GFP) and mutant (PLIN2-206K-GFP) and (PLIN2-309C-GFP) proteins in the HeLa cells wt and mutant. PLIN2 proteins were fused to the green fluorescent protein (GFP), and LAMP2 and LC3 proteins were stained using anti-LAMP2/CY3 and anti-LC3/CY3. Nuclei were stained with Hoechst. Slides were visualized with a Nikon Confocal Microscope A1R at 60× magnification (LAMP2) and with Nikon Eclipse Ni-E 40x magnification (LC3-II). The localization of the PLIN2 wt and mutant proteins is indicated by green spots spread all over the cytoplasm. PLIN2 wt protein is localized in peripheral spots as well as weakly in the nucleus. The PLIN2-206K mutant protein is expressed intensely in the nucleus in several cells, while the PLIN2-309C mutant protein was observed only as peripheral spots.
In summary, we discovered a complex genetic background characterized by a set of rare variants in different genes worthy of future studies to confirm their role in the severity of the LOPD etiopathogenesis. These data highlights that the complexity of the LOPD phenotype is associated to a polygenic inheritance.

Discussion
Precision medicine is an emerging concept that stratifies the patients and their diseases by taking into account a wide array of individual data, including clinical and instrumental information, lifestyle, genetic background and individual drug responsiveness.
In the present study, we described two LOPD families from South Italy manifesting with a complex spectrum of disease symptoms, determining their clinical heterogeneity also within the same family.
The most frequent complaints were related to reduced muscle strength and fatigability, but some patients had to cope with respiratory insufficiency. These symptoms are usually the major contributors in developing severe complications during the disease course [30]. In this study, we generated a severity index score based on phenotypic manifestation of the disease and confirmed that impairment of motor capacity, assessed by GSCG scale as well as FVC and ∆-FVC values at baseline, were associated with a worse clinical outcome. Previous studies demonstrated that secondary factors might substantially influence the clinical course of the disease in LOPD patients, even if they carried the same GAA mutation [10]. In our cohort, alteration of bone metabolism, including vitamin D levels and BMD T-score of the femoral neck, showed a significant relationship with a higher value SI score. These results were supported by the evidence that low BMD is a frequent finding in patients with LOPD and may be causally related to decrease proximal muscle strength [21]. The essential role of vitamin D in muscle development and repair is well known, and the latest investigations demonstrated a genetic contribution of vitamin D to muscle functioning [31]. To our best knowledge, this is the first study demonstrating a significant association between low serum levels of Vitamin D and worse phenotype of LOPD patients.
We performed WES analysis of the affected siblings and identified 16 deleterious variants in 12 genes RILP, FNIP2, TRAPPC11, PLIN2, IRS1, KL, LRP4, RUNX1, FAM26F, TRPV5, TRPV6, and DCHS1. These genes are involved in interconnected pathways, including autophagy/lysosome, immunity, bone metabolism and skeletal muscle disease-related genes, suggesting that rare deleterious variants in these genes or more generally alterations of these pathways might act to worsen the LOPD phenotype.
Autophagy is a dynamic cellular process associated with the pathogenesis of a wide range of human disorders, including LOPD [32,33]. This implies that the deleterious variants in autophagy-related genes may affect disease onset and progression [33]. Recently, it has been suggested that any mutation affecting the autophagosome-lysosome machinery reduces GAA activity as a consequence of the down regulation of GAA synthesis and maturation [34,35]. This implies that the concurrence in each patient of variants in two or more genes acting in various steps of the autophagy-lysosomal pathway might explain differences in the clinical expressivity of the GAA mutations. Remarkably, in this study, we identified five LOPD patients (F1-II4, F1-II5, F1-II9, F2-II4 and F2-II5) showing the highest severity index score, deleterious variants affecting important domains of RILP and FNIP2 proteins involved in lysosome cellular positioning and autophagy [36].
Rab interacting lysosomal protein (RILP) is a 45 kDa lysosomal protein that acts as a downstream effector of Rab7 in phago-lysosome formation and in the transport of the endocytic cargo to the lysosome [37]. Active Rab7 on the phagosomal membrane recruits RILP, which in turn promotes a dynein-dynactin association with the phagosomes, mediating their retrograde moving towards the microtubule-organizing centre (MTOC), where late endosomes and/or lysosomes were clustered [37,38]. Remarkably, the c.G850A (p.G284S) variant, identified in both families investigated in this study, lays into the Rab binding domain and alters a highly conserved residue of the protein, suggesting that this change might disturb the formation of the RILP-Rab complex. These bioinformatics predictions were supported by a different pattern of expression and localization observed in skeletal muscle of the patients carrying the mutated allele and in HeLa cells. The RILP-284S protein was not efficiently recruited at autophagosome-lysosome and remains spread all over the cytoplasm, suggesting a perturbation of the autophagic flux.
Moreover, during nutrient insufficiency, the FLCN-FNIP complex modulates the formation of a FLCN-RILP-Rab complex for promoting the clustering of lysosomes around the nucleus [36].
Folliculin interacting proteins 2 (FNIP2) plays a crucial role in the regulation of the autophagic flux by modulating the formation of the FLCN-GABARAP complex [39]. Moreover, FLCN and FNIP proteins each contain both a longin and a differentially expressed in normal versus neoplastic cells (DENN) domain [40], which are implicated in the regulation of small GTPases and membrane trafficking. In fact, the FLCN-FNIP complex regulates both the Rag and Rab GTPase families [41,42], which in turn modulate the key mTORC1 signalling pathway and lysosomal distribution, respectively, in a manner dependent on amino acid availability [36]. Noteworthy, the two rare and deleterious variants identified in the two LOPD families alter the longin and DENN domains of the protein, suggesting that they might disturb the formation of this complex.
Taken together these evidences lead to hypothesize that alteration of the genes involved in lysosome positioning and autophagy might influence lysosome metabolism, worsening the effect of the GAA mutations in the two investigated families.
Separate consideration needs to be given to patient II-5 of family 2, who manifested the most severe phenotype observed in our cohort of patients, characterized by severe muscle and pulmonary impairments (6MWT = 10, GSCG = 26, FVC = 36% and vacuolated muscle fibres = 15%). Remarkably, this patient inherited, together with the two mutations in GAA gene and the two rare and deleterious variants in the RILP and FNIP2 genes, two compound heterozygous variants c.A1501C (p.N501H) and c.T5165A (p.L1722H), in the low-density lipoprotein receptor-related protein-4 (LRP4) gene. LPR4 is a major mediator of postsynaptic stability at the neuromuscular junction, a region critical for control of muscle contraction [43,44]. Recessive mutations in this gene are linked to congenital myasthenic syndrome (CMS17), a severe neuromuscular disorder associated to respiratory failure [43]. Respiratory related symptoms (dyspnoea on effort, reduced physical capacity and recurrent infections) and respiratory failure are often evident in the later stages of LOPD [45]. The two rare and deleterious LRP4 variants identified only in this patient alter the highly conserved residues laying in important domains of the protein (Figure 4). In this context, we hypothesize that these variants, although reported as likely benign in the Global Variome shared LOVD database (https://databases.lovd.nl/shared/genes (accessed on 5 March 2021)), might have worsened the entire clinical picture of this patient. Moreover, she also shares with her brother II-4 a rare deleterious variant in TRAPPC11, whose mutations cause limb girdle muscular dystrophy 2S (LGMD2S), congenital muscular dystrophy (CMD) and dystroglycanopathy. Interestingly, this gene has a role upstream of autophagosome formation, by recruiting the ATG2B-WIPI4/WDR45 complex to preautophagosomal membranes, suggesting that the TRAPPC11 variants might also perturb autophagy flux [46,47].
Two rare and deleterious variants disturbing the lipid-binding domain of the PLIN2 gene segregated with obesity in the two families. PLIN2 encodes the most abundant lipid droplet (LD)-associated protein in non-adipose tissue, and its expression correlates with intracellular lipid accumulation [48,49]. Recent studies reported that mutations in PLIN2 affect lipolysis and are associated with reduced plasma triglyceride levels in humans [48], as well as with reduced insulin secretion and increased insulin sensitivity in obese Italian subjects [50]. Moreover, lipids in LDs were reportedly required for the induction of autophagy in yeast and served as the substrate of autophagosomes [51]. This aspect leads to the hypothesis of a link between obesity and autophagy impairment in LOPD patients. Here, we have demonstrated that RILP, FNIP2 and PLIN2 proteins localized at the sarcolemmal or sarcoplasmic level, with the same localization of markers of the autophagic-lysosomal system (LC3, LAMP2). Differences in the pattern of distribution in mutated versus wild type biopsies and even more of these versus external controls were evident. In particular, the different protein distribution between LOPD patients and the external control affected by vacuolar myopathy not linked to GAA mutations is very interesting.
Recent studies reported that bone involvement is a common finding in Pompe disease, particularly in wheel-chair users, primarily characterized by low BMD and scoliosis. LOPD patients have a high risk of fragility fractures, chronic immobilization and muscle weakness and BMD reduction, particularly at the femoral neck, seems to be related to lower mechanical loading resulting from structural and functional impairment of the primary muscles involved in walking, such as hip flexors and knee extensors muscles [21,52]. Vitamin D plays a key role in the bone-muscle crosstalk [53]. In this context, its effects consist mainly of the maintenance of intracellular Ca2 + homeostasis that is of crucial importance for muscle contraction and bone metabolism. In the present work, deleterious variants in IRS1, KL, RUNX1, TRPV5 and TRPV6 genes, involved in bone and mineral metabolism, were identified in family 1 [54][55][56][57][58][59][60][61][62][63][64][65][66].
In this context, we hypothesize that the load of polygenic rare and deleterious variants in genes closely related to the Ca2+ level and bone metabolism together with a perturbed autophagy, which plays an important role in bone remodelling, might explain the different spectrum of bone compromising in the family members.
Since 2006, ERT has been available for Pompe disease. ERT improved survival and motor function, normalized left ventricular cardiac hypertrophy and stabilized respiratory capacity [67]. However, such a treatment has brought to light the role of the immune response in abrogating the efficacy of the rhGAA treatment. Several patients are also cross-reactive immunological material negative (CRIM-) and develop high titre immune responses to ERT with rhGAA [68]. Thus, optimization of treatment for such patients includes development and utilization of strategies to prevent or eliminate immune responses, including modulating the immune system (prophylactic and therapeutic immune tolerance induction regimens) and engineering the enzyme to be less immunogenic and more effective.
In this study, we identified a deleterious variant in the FAM26F gene in the patient II-5 who manifested immuno-intolerance to rhGAA. In silico analysis showed that the splicing variant c.525 + 2T > G disrupted the donor splicing site of intron 2, producing a predicted truncated protein of 224 amino acids lacking the fourth transmembrane domain and of the cytosolic C-terminal tail of the protein [69]. The family with sequence similarity 26, member F (FAM26F), is a recently identified tetraspanin-like membrane glycoprotein involved in various immune modulating responses, including homophilic interactions and potential synapses between several immune cells, including CD4+, CD8+, NK, dendritic cells and macrophages. Several studies have demonstrated that FAM26F is highly crucial in determining/dictating the immune response of an individual, also in diseased conditions [70]. The ImmuNet bioinformatics tool predicts a high relationship confidence with the C1S gene that codifies a major constituent of the human complement subcomponent C1.
We hypothesize that this variant might affect the expression of the FAM26F protein, thus representing a good candidate to explain the adverse reaction to ERT associated with the C1 protein inhibitor in patient II-5. We suggest measuring the C1 inhibitor level before ERT treatment in LOPD patients.
Moreover, we demonstrated the important role of bone metabolism, in particular of vitamin D, as a key contributor in determining a worse phenotype and severe complications, including bone fractures. More studies are needed to confirm this association. However, we suggest that vitamin D serum levels and BMD should be monitored at regular intervals, and, if necessary, vitamin D should be supplemented.

Patients
Patients were recruited at the Reference Centre for Neurological and Neuromuscular Rare Disease, University of Campania "Luigi Vanvitelli" (Figure 1A,B) [18].
All clinical data were collected at baseline, excepted for immune response to ERT, for which data were referred to a period after one year from the start of therapy. Motor function was assessed by MMT using the Medical Research Council (MRC) grading scale for different muscle groups, and the 6MWT was used to evaluate functional endurance during prolonged ambulation. The GSGC scale was performed as a measure of global motor disability. Respiratory assessment included spirometry for evaluation of FVC, considering as normal value FVC > 80% of predicted, and measurement of ∆-FVC, defined as drop of FVC from stand up to supine position (∆-FVC > 20% is suggestive of restrictive pattern), according to the American Thoracic Society standards [20]. BMI was measured to provide an estimate of body fat in our patients. We defined a status of normal weight (19.5-24.9 kg/m 2 ), overweight (25-29.9 kg/m 2 ), obesity class I (30-35.9 kg/m 2 ), obesity class II (36-40 kg/m 2 ) and obesity class III (>40 kg/m 2 ), according to WHO recommendation for Caucasian subjects [71]. Vitamin D serum levels (normal values: >30 ng/mL) were measured and the BMD of the lumbar-sacral tract, femoral neck and total-body (T-score normal value > −1) was evaluated by computerized bone mineralometry (MOC); both parameters were used for assessment of bone metabolism. History of bone fractures was considered. Basilar artery dolichoectasia was detected by brain MRI, and structural cardiac valvular defects were assessed by M-MODE/2D/Color-doppler cardiac ultrasound. Immunity state was determined using IgG and IgE antibodies against glucosidase alpha (anti-rhGAA); an adverse reaction to ERT was taken in account. After informed consent, muscle biopsy was obtained from all patients. Specimens from each biopsy were processed according to standard procedures for histology and biochemistry studies. Measurement of percentage of vacuolated fibres on muscle samplings, performed by optical microscopy on 3 fields at 20× on ATP pH 9.4 staining, and quantification of lymphocytes containing PAS-positive granules on blood smears [72] were carried out as indicators of impaired autophagy, together with assessment of GAA activity on DBS. All DBS were performed at the diagnostic service of Metabolic Laboratory, Centre of Diagnostics (Hamburg, Germany) by Dr. Zoltan Lukacs. The activity value of GAA was measured at pH 3.8, pH 7.0 and after specific inhibition. Normal GAA activity was considered between 1.86 and 21.9 mol/h/L. Beta-galactosidase activity was used as the internal control. Finally, venous blood samples were taken from patients for genomic DNA extraction and genetic analysis.
The alignment of the 100-bp paired-end reads to the human reference genome was performed by using the Burrows Wheeler Aligner (BWA) MEM, v0.7.542. After removal of duplicate reads through the Picard Mark Duplicates command (with standard options), we called the single nucleotide variants (SNVs) and insertions/deletions (indels) for all samples using HaplotypeCaller (BP_RESOULTION option) and GenotypeGVCFs in Genome Analysis Toolkit (GATK), v3.5-0-g36282e4, following the manufacturer's best practice guidelines (available at https://software.broadinstitute.org/gatk/best-practices/ (accessed on 5 March 2021)) [73]. Variants with a Minor Allele Count (MAC) = 0, number of alternative alleles =2 and call rate <95% were also filtered out, as well as samples with identical-bydescent (IBD) sharing and sex mismatches, and samples with a call rate <90%. Variants passing quality control were annotated to genes (within 10 kb from transcription start/stop site) through ANNOVAR [74]. Functional prioritization of the candidate genes was performed by STRING database analysis. The entire list of 589 candidate genes was analysed with the STRING database to search for protein-protein interaction (PPI) networks and for metabolic pathways enrichment. We found a significant PPI enrichment (p value = 1.7 × 10 −6 ), which means that selected proteins have more interactions among themselves than what would be expected for a random set of proteins of similar size, drawn from the genome. In particular, we found a significant enrichment for genes related to skeletal muscle development and disease (p = 0.04) and in family 1 we identified a cluster of 18 gene products involved in calcium and bone metabolism (p = 0.03). Moreover, we analysed the metabolic pathways annotated by STRING to select the genes involved in the autophagy/lysosomal pathways as well as in immunity. Information were verified in OMIM and PubMed and selected variants were confirmed by direct sequencing, as already published [18].

Immunofluorescence Staining of Skeletal Muscle
For immunofluorescence (IF) analysis, skeletal muscle cryosections (7 µm-thick) were fixed in pre-cooled acetone (−20 • C) for 5 min at room temperature (RT). Endogenous peroxidase activity was blocked by means of a 0.3% hydrogen peroxide and methanol solution for 5 min at RT. After washing, incubation with a blocking buffer (1 × PBS/5% horse serum) was performed for 1 h in a humidified chamber. Primary and then secondary antibody incubations were performed at 4 • C overnight and at RT for 2 h, respectively. Anti-LYAG Muscle biopsies from a non-GAA vacuolar myopathy patient (CNT1) and from a non-myopathic orthopaedic patient (CNT2) were used as external controls, being available in the tissue bank of our department.

Plasmids Preparation, Cell Cultures and Transfection
The wild-type (wt) expression vectors of RILP and PLIN2 were purchased from the Addgene repository (EGFP-RILP catalogue: 110498; pEGFP-C1-ADRP (PLIN2) catalogue: 87161) and were used as the templates to introduce the RILP (c.G850A, p.G284S) and PLIN2 (c.G616A, p.E206K and c.C925T, p.R309C) variants. Mutagenesis were carried out by the Quick Change II XL site-directed mutagenesis technique according to the manufacturer's instructions (Agilent Technologies, Santa Clara, CA, USA). The whole coding sequences of wild-type and mutant plasmids were sequenced to confirm the mutagenesis and to exclude undesired mutations.
In total, 500 nanograms of each vector were incubated with 1 microgram of Lipofectamine 2000 (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) in Opti-MEM (Gibco, Thermo Fisher Scientific, Waltham, MA, USA) for 20 min; the mixture was added to HeLa cells (80,000 cells in 24 wells plate on coverslip) that were grown in complete medium for 20 h. Cells to be stained with anti-p62 and anti-LC3B antibodies were starved for 1 h in Hanks' Balanced Salt solution (HBSS) before transfection. Twenty hours after transfection, cells were fixed with paraformaldehyde (PFA) 4% and stained with anti-LAMP2, anti-p62 and anti-LC3B antibodies. Briefly, incubation with blocking buffer (1×PBS/0.1% triton/5% serum) was performed for 15 min in a humidified chamber. Primary and then secondary antibody incubations (1 × PBS/0.1% triton/antibody) were performed at 4 • C overnight and at RT for 1 h, respectively. Anti-LAMP2 (ab25631, Abcam, Cambridge, UK To assess the effect on splicing of the variant c.525 + 2T > G, located into the intron 2 of the FAM26F gene, we performed bioinformatics analysis by using the NetGene2 server (http://www.cbs.dtu.dk/services/NetGene2/ (accessed on 5 March 2021)), which produces neural network predictions of splice sites in human. We found that the c.525 + 2T > G variant completely disrupted (100% of probability) the canonical donor splicing site in intron 2.
To examine the expression pattern of the alternative transcript generated by the c.525 + 2T > G variant, an RT-PCR assay was performed using specific couples of oligonucleotides designed on exons 2 and 3 and on intron 2 of the gene. Blood RNA was isolated by using Tempus Spin RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA). A total of 1 µg of total RNA, digested with DNasi RNasi free (Thermo Fisher Scientific, Waltham, MA, USA), was reverse transcribed with the Superscript III-First strand kit (Thermo Fisher Scientific, Waltham, MA, USA). Quantitative PCR (qPCR) reactions were performed in triplicate, using human FAM26F transcript-specific primers and ITaq Universal Sybr Green Supermix (Bio-Rad, Hercules, CA, USA) following the manufacturer's directions. Results were normalized versus the expression of the beta-actin gene. The SD was calculated by using data of three different experiments. PCR products where directly sequenced on both strands using the Big Dye Terminator Ready Reaction Kit (Applied Biosystems, Foster City, CA, USA).

Statistical Analysis for Clinical Variables
Characteristics of the study population were presented using descriptive statistics. The mean and SD values were calculated for continuous variables, while frequencies were reported for categorical variables. Patients were classified as having a higher severity score (if moderate/severe or very severe score) or lower severity score (if mild or very mild score). Pearson correlation was used to explore the relationship between the severity score and clinical variables. A non-parametric test for independent samples (Mann-Whitney U) was used to compare the quantitative clinical variables (age, AAO, MMT-MRC, 6MWT, FVC, ∆FVC, BMI, VDD, BMD total body, BMD femoral neck, BMD lumbar-sacral and GSGC) in the two populations. Chi-squared tests were used to compare the qualitative variables (sex, BF, BAD, MVP, IgG-rhGAA and ERT-AE) between the two subgroups. A p value < 0.05 was considered statistically significant.

Conclusions
The findings of this study strongly support the hypothesis that disease severity and progression in LOPD patients is caused by the concomitant mutations of GAA and other autophagy genes. Of course, at this stage of advancement of our study, we cannot offer any improvement in practicing medicine; however, we believe it is important to report the two main results that emerge from this study. First is that the majority of the identified variants were located in autophagy-related genes, and second, some of these mutated genes were shared between the two LOPD families. Undoubtedly, the identification of potential modifier genes not shared by the majority of LOPD patients could cause great uncertainty on the interpretation of the clinical phenotype and of the ERT results collected in different series so far. However, at the molecular level, these findings offer a rational support to address additional genetic studies to identify in a larger cohort of LOPD patients the functional variants that could affect the autophagy pathway. A detailed molecular definition of Pompe disease would be essential to improve the treatment of the patients, in line with a personalized medicine approach.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The datasets analyzed during the current study are available from the corresponding author on reasonable request.