BBS is an extremely clinically and genetically heterogeneous ciliopathy characterized by intra- and inter-familial variability2,20. Generally, BBS presents with an autosomal recessive inheritance, but in some cases, an oligogenic inheritance has been proposed, in the form of triallelism or second-site modifiers8,11,12,13,21. In these families, some unaffected or less severely affected individuals have two pathogenic variants in a BBS-related gene, whereas the BBS-diagnosed or more severely affected relatives carry three alleles in two different BBS-related genes. Although the involvement of triallelism in BBS families is controversial8,9,10,13, there is further evidence for the possible existence of third modifier alleles11,12,13,14. Hence, our study presents a comprehensive study of the mutational burden in BBS and highlights the importance of considering non-Mendelian inheritance to improve the clinical management of BBS.

First, we recruited a large cohort of 99 cases from 77 families with genetic suspicion of BBS after molecular testing. Although other studies only consider those cases that met the diagnostic criteria for BBS10,12,13, we have included also patients who did not fulfill a priori the clinical criteria of BBS described by ref. 1. In fact, half of the syndromic cases with available clinical information (41/83) were not first classified as BBS but after genetic testing. Thus, for some of the clinical features, a slight bias was observed in our cohort compared to the distribution in other reports1,2,22, which can be explained by the inclusion of cases with a diagnosis of “BBS-like” or even more unspecific systemic findings (“RD-OTHERS”), or a poor clinical acquisition of some of the features.

Nowadays, NGS is the technology of choice for the study of BBS23,24,25. It allows the identification of new causative variants and further reanalysis to assess new BBS genes that may have been identified after a primary analysis or had not been covered by any of the classical methods initially used26,27. This statement is consistent with the fact that 78% of the suspected oligogenic families in our cohort were discovered through NGS approaches.

It has been estimated that oligogenic inheritance is present in less than 10% of the BBS families28. This value is confirmed by our triallelic distribution, with a rate of 13% in the informative triallelic families. However, it increases to 51% when all families with suspected modifiers alleles were included. We cannot elucidate if the third allele triggers a modifier effect in most families, because both siblings presented the same genotype-phenotype, the index case was the only affected in the family, or clinical information was unavailable.

BBS1 can cause both BBS and non-syndromic inherited RD. This gene is the most frequent source of BBS, accounting for 23–51% of characterized families12,29,30, which agrees with 42% of molecularly characterized families with BBS1 as the major primary gene in our cohort. The variant p.(Met390Arg) has been reported in up to 80% of BBS1-related alleles across different worldwide populations29, being mostly associated with BBS as only 21% of BBS1-positive patients showed non-syndromic presentations31. In our cohort, this variant is accordingly the most represented with an allelic frequency of 70%. The reason for this slight decrease may be that we only included BBS-associated families but excluded those with only non-syndromic affected individuals.

The implication of BBS1 in oligogenic inheritance remains unclear. Some reports claim that BBS1 is rarely involved in complex inheritance29, but BBS1 has been reported in a triallelic inheritance in 15% of families, being the primary gene instead of acting as a modifier in only 4% of cases12. However, our results go further, showing that only 65% of our BBS1-characterized families fit in an autosomal recessive inheritance, while in the remaining families, BBS1 might participate in oligogenic inheritance as the primary gene in 11 families, together with MKKS, IFT172, or other six BBS-related genes. This same behavior was observed for BBS12 in three families from our cohort, but in this case, there are no reports related to the likelihood of its participation in complex inheritance. In those families, biallelic BBS12 variants were found together with the missense variants p.(Gly352Val) in PDE6B, p.(Arg82Leu) in TRIM32 (MIM *602290), and p.(Arg400Cys) in SDCCAG8 (MIM *613524), respectively. While the combinations BBS12_SDCCAG8 and BBS12_TRIM32 were predicted only with DiGePred, only the combination BBS12_PDE6B was predicted as oligogenic with ORVAL platform with a 99% of confidence.

Generally, BBS-associated proteins are located at the base of the cilium and participate in ciliary biogenesis and in cilia function1,32, but a variety of specific locations and functions have been described33,34. In our work, the involvement of BBS proteins in a specific complex, structure, or process does not seem to be related to the level of involvement of each gene in the triallelic inheritance. However, we have seen that the genes that encode for the chaperonin-like complex (BBS10, BBS12, and MKKS) are mostly implied in a recessive inheritance. Moreover, BBS12 and MKKS are normally the primary genes when involved in oligogenic inheritance. Therefore, these genes are usually the principal gene causing the disease regardless of the type of inheritance. It has been reported that the activity of genes encoding BBSome components (ARL6, BBIP1, BBS1, BBS2, BBS5, BBS7, BBS9, and TTC8) may be dependent on the chaperonin-like genes35 and families with variants in the chaperonin-like complex present a more severe phenotype36, so these chaperonin-like genes may not normally require a second-site modifier. Nevertheless, we found three families from our cohort carrying the variant p.(Ala242Ser) in MKKS as a possible modifier allele in combination with biallelic BBS1 variants. This non-synonymous change, which has been previously described as a dominant-negative allele17,37, disrupts the protein conformation of the BBSome, thus preventing them from doing their proper function17,37. Intrafamilial variable severity has been seen in one of these three families carrying the p.(Ala242Ser) variant in combination with homozygous p.(Met390Arg) in BBS1. Our result suggests that the increase in the detected mutational load may correlate with a more severe phenotype, which could be explained by its chaperone function38. Variable expressivity between siblings involving other genetic combinations can be also found in the other five families in our cohort.

The presence of modifier alleles can determine the phenotype, since they may influence the presentation of the BBS phenotype17. This can be the scenario for family RP-2634, in which the syndromic index case has biallelic variants in BBS1 and the heterozygous missense p.(Arg207His) in BBS5, whereas her non-syndromic sister is just biallelic for the BBS1 variants. This BBS5 variant, presenting a minor-allele frequency of 0.9% in Europeans, has been predicted as a null mutation17.

Specific heterozygous variants acting as modifiers has previously been associated with the existence of particular findings (e.g., ocular, neurological, or renal features)39,40,41. A sibling from family GBB23 carries the homozygous pathogenic variant p.(Met390Arg) in BBS1, and additionally, the index case suffering from T2DM also carries the heterozygous variant p.(His3880Tyr) in ALMS1. Mutations in ALMS1 are the cause of Alström syndrome (ALMS; MIM #203800), an ultra-rare metabolic ciliopathy associated with severe visual impairment, sensorineural deafness, obesity, insulin resistance, T2DM, and hypogonadism, among other features42. One of the explanations of glucose metabolism alterations in this syndrome are defects in the ALMS1 protein, which participates in the insulin-regulated glucose transport43. In our family, the variant identified in ALMS1 may be acting as a second-site modifier altering the possibility of suffering T2DM. However, the high frequency of diabetes mellitus in the general population44 might also be a plausible explanation for its presence in this case.

Some of the possible modifier alleles in the BBS-causing genes found in our cohort could be good candidates for functional studies to analyse a possible modifying effect on the BBS phenotype, e.g., IFT172, TRIM32, or WDPCP. Furthermore, we identified third alleles in other genes previously reported as possible candidates or modifiers of BBS, e.g., ALMS1, CORO2B, NPHP4, or PDE6B36,45,46,47. Therefore, our findings could also support a possible involvement of these genes in the pathogenesis of BBS. For example, PDE6B is a gene associated with non-syndromic RD but also reported to BBS phenotype in a consanguineous family with homozygous pathogenic mutations in BBS10 and PDE6B47. Here, three different heterozygous variants in PDE6B were found in three families (family IDs RP-1378, RP-0496, and B64), accompanying biallelic variants in BBS1, MKKS, and BBS12, respectively. It is unknown how the effect of the mutational load detected in these families may be affecting their phenotype, but PDE6B and other genes, which are involved in the phototransduction and visual transduction pathways, are downregulated in BBS and ALMS zebrafish models and may be drivers of the retinal degeneration48. Nonetheless, the hypothetical role of PDE6B and other unclear modifiers in BBS should be further studied functionally.

The effect of modifier alleles on clinical manifestation needs to be assessed usually by means of in vitro and in vivo strategies16,17. However, due to technical limitations, it is not always possible to perform these analyses in a clinical setting. Alternatively, in silico tools can help to discover and predict combinations of variants that can be affecting the patients’ phenotype18,19. According to our data, ORVAL and DiGePred together support an oligogenic inheritance in 44% of our families with more than two alleles, reaching 91% with only one positive method.

To understand the phenotypic variability in BBS families, both genetic and environmental factors should be considered. Recently, multi-omics analyses are being considered to understand and elucidate the role of the mutational load in BBS-associated mechanisms through the integration of multiple analyses (mutational load, differential gene and/or protein expression, epigenetic, and/or metabolome-based signatures). Therefore, obtaining larger data sets would help to clarify the role of possible modifiers in BBS and related ciliopathies23,49.

In summary, this work deepens into the controversial topic of oligogenic inheritance in BBS, finding new evidence for the existence of second-site genetic modifiers as a cause of intrafamilial variability in this disease. Besides, it highlights the importance of the use of NGS in the genetic diagnosis of BBS.