Segregation Analysis of Rare NRP1 and NRP2 Variants in Families with Lymphedema

Neuropilins are transmembrane coreceptors expressed by endothelial cells and neurons. NRP1 and NRP2 bind a variety of ligands, by which they trigger cell signaling, and are important in the development of lymphatic valves and lymphatic capillaries, respectively. This study focuses on identifying rare variants in the NRP1 and NRP2 genes that could be linked to the development of lymphatic malformations in patients diagnosed with lymphedema. Two hundred and thirty-five Italian lymphedema patients, who tested negative for variants in known lymphedema genes, were screened for variants in NRP1 and NRP2. Two probands carried variants in NRP1 and four in NRP2. The variants of both genes segregated with lymphedema in familial cases. Although further functional and biochemical studies are needed to clarify their involvement with lymphedema and to associate NRP1 and NRP2 with lymphedema, we suggest that it is worthwhile also screening lymphedema patients for these two new candidate genes.


Introduction
Neuropilins 1 and 2 are transmembrane coreceptor proteins encoded by the NRP1 and NRP2 genes. They were first identified in the nervous system of Xenopus laevis embryos [1,2]. Although they map to different chromosomes (10p11.22 and 2q33.3, respectively; see Table 1), they have 45-50% amino acid sequence homology [3]. The two neuropilins are implicated in angiogenesis, and both are expressed by endothelial cells and neurons [4].
The two neuropilins have a very short intracellular domain (only about 40 amino acids long) and a longer extracellular region. Both neuropilins complex a large variety of ligands, by which they take part in cell signaling. These transmembrane proteins have a critical role in the regulation of vascular and neural development through binding vascular endothelial growth factors (VEGFs) and semaphorins [5]. The major neuropilin ligands are vascular endothelial growth factors VEGF-A, VEGF-C and VEGF-D [4,6]. While VEGF-A plays a primary role in angiogenesis, VEGF-C and VEGF-D are important for normal development of the lymphatic system. NRP1 and NRP2 also bind other ligands-for instance, hepatocyte growth factor (HGF) [7], EGF receptor [8] and an adhesion receptor L1-CAM [9]. Although they mediate cell signaling by ligand binding, it is unclear whether they also transduce signals on their own. In any case, they take part in many biological processes, and changes in their function may modulate tumor growth and metastasis, affecting vascular and lymphatic biology [10]. This can give rise to pathologies such as cancer and autoimmune disorders [11].
While NRP1 is predominantly expressed in the arterial endothelium and lymphatic valves, NRP2 mainly functions in venous and lymphatic endothelial cells (LECs) [12]. This indicates that NRP1 acts as a semaphorin receptor in the valve endothelium throughout lymphatic vessel development [12]. In contrast, NRP2 is implicated in VEGF-C-driven lymphatic vessel growth [13], while double-heterozygous nrp2 +/− vegfr3 +/− mice show a decrease in lymphatic vessel sprouting in adult organs (Table 1) [13].
Most information about neuropilin function comes from studies with mice. Early studies of Nrp1 −/− mice showed deficient neuronal vascularization, aortic arch malformations and impaired yolk sac vascularization [14,15]. Abnormal blood vessel formation due to overexpression of the Nrp1 gene was also observed [16], although the authors of the study did not investigate the lymphatic phenotype.
Nrp1 was recently implicated in lymphatic valve development via Sema3a/Nrp1 signaling [12]. Sema3a is a ligand bound by Nrp1. Deficiency or blockage of Nrp1 bound to Sema3a leads to malformations of a specific type of LEC that later migrates from the wall of lymphatic vessels to give rise to lymphatic valves. Defects in this process lead to abnormal lymphatic valve development [13,17,18].

Lymphatic phenotype in mouse model
Nrp1 −/− malformation of valve precursor LECs [18] Nrp2 −/− reduction or loss of small lymphatic vessels and capillaries [19] Nrp2+/− and Vegfr3+/− abnormal lymphatic development and reduced lymphatic vessel branching [13] In comparison, the Nrp2 −/− mouse model provided evidence of the selective role of Nrp2 in the development of lymphatic vessels. Nrp2-deficient mice showed significantly fewer or a complete absence of small lymphatic vessels and capillaries in many organs. The small number of vessels that did develop were wrongly positioned, and in some cases, small vessels and capillaries were larger in Genes 2020, 11, 1361 3 of 17 size. A decrease in LEC proliferation was also observed, indicating a role of Nrp2 in the proliferation of LECs in selected types of lymphatic vessels and implying the control of positional guidance by Nrp2 [19].
Knocking out Nrp1 and Nrp2 separately has evident effects on the mouse phenotype, and the joint deletion of Nrp1 and Nrp2 has much graver consequences. Takashima et al. [20] reported that Nrp1 −/− /Nrp2 −/− mice died at an early stage of embryogenesis (E8.5) and displayed severely defective blood vessel development in the yolk sac. In the study, the authors observed that double-knockout mice had severer vascular malformations than Nrp1or Nrp2-deleted animals. Surprisingly, heterozygous mice with Nrp1 +/− /Nrp2 −/− or Nrp1 −/− /Nrp2 +/− genotypes did not survive embryo development either but lived longer than double-knockout mice (E10-E10.5) [20]. Although the authors did not study the development of lymphatics in these double-knockout mice, the premature death of the embryos and the fact that Nrp1 and Nrp2 deletions caused severe lymphatic system malformations suggested that a double knockout also results in impaired lymphatic system development.
The development of small lymphatic vessels, capillaries and valves is important, and its impairment may lead to an overaccumulation of lymph in the tissues and lymphedema. Lymphedema is a pathological condition characterized by swelling, usually at the extremities, accompanied by pain and inflammation [21]. Lymphatic capillaries, the development of which is regulated by NRP2, are indispensable for normal drainage of interstitial fluid into the lymphatic system. Lymphatic capillaries lack a basement membrane and feature endothelial cell junctions that grant the passage of interstitial fluid into the lymphatic system [22]. When lymphatic capillaries and small vessels are malformed, interstitial fluid cannot enter lymphatic collecting vessels and accumulates in the tissues. The proper functioning of valves is also crucial to ensure retrograde lymph flow. Valve defects cause lymph flow obstruction and fluid accumulation [22]. The blockage of lymph flow and lymph accumulation are the primary causes of lymphedema.
In our study, a cohort of 235 Italian lymphedema patients, who tested negative for variants in known lymphedema-associated genes, was screened by a next-generation sequencing (NGS)-targeted panel for variants the in NRP1 and NRP2 genes. A segregation analysis in familial lymphedema cases added evidence that NRP1 and NRP2 qualify as candidate genes to include in the genetic testing of lymphedema patients.

Clinical Evaluation
We retrospectively enrolled 235 Caucasian lymphedema patients in this study, which was conducted in accordance with the Declaration of Helsinki, and whose protocol was approved by the Ethics Committee of Azienda Sanitaria dell'Alto Adige, Italy (Approval No. 94-2016). No consanguinity was reported in any family. Clinical diagnosis of lymphedema was made in hospitals across Italy, according to generally accepted criteria. All subjects gave their informed consent for inclusion before they participated in the study.

Genetic Analysis
Genetic testing was performed on genomic DNA extracted from saliva or peripheral blood of probands, as previously described [23]. Briefly, a custom-made oligonucleotide probe library was designed to capture all coding exons and flanking exon/intron boundaries (±15 bp) of genes known to be associated with lymphedema [24]. We added the candidate genes NRP1 and NRP2 to our panel. DNA from probands was analyzed for germline variants. Variants were confirmed by bidirectional Sanger sequencing on a CEQ8800 Sequencer (Beckman Coulter, Brea, CA, USA).
We searched the international Single Nucleotide Polymorphism Database (dbSNP) and the Human Gene Mutation Database professional for all nucleotide changes. In-silico evaluation of the pathogenicity of sequence changes in NRP1 and NRP2 was performed using the Variant Effect Predictor Genes 2020, 11, 1361 4 of 17 tool and MutationTaster. Minor allele frequencies were checked in the Genome Aggregation Database (gnomAD, https://gnomad.broadinstitute.org/). All variants were evaluated according to the American College of Medical Genetics and Genomics guidelines [25]. Detailed pretest genetic counseling was provided to all subjects, who were then invited to sign informed consent to use of their anonymized genetic results for research.

In-Silico Analysis
The primary amino acid sequences of NRP1 and NRP2 in FASTA format ( Table 2) were used as targets to search template libraries (i.e., the SWISS-MODEL Template Library (version 2019-10-24) and the Protein Data Bank (release 2019-10-18)) [26] for matching evolution-related structures by means of BLAST [27] and HHBlits [28], as previously described [23]. Briefly, models were based on target-template alignment using ProMod3 of the SWISS-MODEL server [29]. An alternative model was obtained with ProMod-II in those cases in which loop modeling failed with ProMod3 [30]. Coordinates conserved between the target and the template were copied from the template to the model. Insertions and deletions were remodeled using a fragment library. Side chains were then rebuilt. Finally, the geometry of the resulting model was regularized with the CHARMM27 force field [31]. Global and per-residue model quality was assessed using the QMEAN scoring function [32]. BioVia Discovery Studio Visualizer (version 17.2) [33] was used to visualize the modeled protein, to vary the targeted amino acids and to analyze interactions at the molecular level.

Clinical and Genetic Evaluation
DNA from 235 Caucasian patients diagnosed with lymphedema, from different Italian hospitals, previously testing negative for variants in known lymphedema and lymphatic malformation-causing genes, was screened for variants in new candidate genes. Here, we report the results for NRP1 and NRP2. No consanguinity was reported in their families.
NGS screening detected two different variants in NRP1 in two families and four different variants in NRP2 in four families. All the variants were heterozygous. The clinical features of the probands and family members are shown in Table 3. A segregation analysis performed on family members carrying variants in NRP1 are shown in Figure 1. The pedigrees of families with variants in NRP2 are shown in Figure 2.  Regarding variants in the NRP1 gene, the first proband, male, 68 years, has had lymphedema of the lower limbs since his youth. The case is sporadic, and no other family members of the proband were tested. The proband carries a missense single-nucleotide variation c.2557A > C (dbSNP rs760388137) with a GnomAD-reported frequency of 0.0001178.
The second proband with an NRP1 variant, female, 54 years, has stage 2 lymphedema of the left foot and ankle, diagnosed at age 27. The case is familial: the proband's daughter also has stage 2  Regarding variants in the NRP1 gene, the first proband, male, 68 years, has had lymphedema of the lower limbs since his youth. The case is sporadic, and no other family members of the proband were tested. The proband carries a missense single-nucleotide variation c.2557A > C (dbSNP rs760388137) with a GnomAD-reported frequency of 0.0001178.
The second proband with an NRP1 variant, female, 54 years, has stage 2 lymphedema of the left  Regarding variants in the NRP1 gene, the first proband, male, 68 years, has had lymphedema of the lower limbs since his youth. The case is sporadic, and no other family members of the proband were tested. The proband carries a missense single-nucleotide variation c.2557A > C (dbSNP rs760388137) with a GnomAD-reported frequency of 0.0001178.
The second proband with an NRP1 variant, female, 54 years, has stage 2 lymphedema of the left foot and ankle, diagnosed at age 27. The case is familial: the proband's daughter also has stage 2 lymphedema of the lower limbs ( Figure 1). Both carry the same heterozygous missense variant c.1655G > A (dbSNP ID rs757990959) with a frequency of 0.00004873 (GnomAD) classified variant of unknown significance (VUS) by the ACMG [34].
Regarding variants in the NRP2 gene, we identified four families with rare variants. In the first family, the proband, female, 18 years, has congenital edema of the left hand (lymphedema stage 2). Indeed, upper limb lymphoscintigraphy failed to visualize the left axillary lymph nodes. She carries a missense variant c.580T > C (rs755679361), which has an allele frequency in healthy control populations of 0.000008122, according to GnomAD. The segregation analysis in available family members showed that the proband inherited the variant from her father. The father was apparently healthy; however, lymphoscintigraphy showed clear bilateral lymph transport problems, especially on the left side ( Figure 3).  The second proband with an NRP2 variant, male, nine years, developed edema of the left foot, ankle and heel after vaccination at age three months. The case is sporadic, and no other family members were tested. The missense variant is classified as VUS (c.1748T > C) (dbSNP rs746130411), with a GnomAD frequency of 0.000004068.
The third proband, female, 45 years, was diagnosed with edema of the right foot at age 32. The The second proband with an NRP2 variant, male, nine years, developed edema of the left foot, ankle and heel after vaccination at age three months. The case is sporadic, and no other family members were tested. The missense variant is classified as VUS (c.1748T > C) (dbSNP rs746130411), with a GnomAD frequency of 0.000004068.
The third proband, female, 45 years, was diagnosed with edema of the right foot at age 32. The heterozygous missense variant is c.838C >T (dbSNP rs79750907) and has a frequency of 0.001783 (GnomAD). This is a familial case, since the proband's son was also diagnosed with left foot lymphedema at age nine. The son was found to carry the same variant as the proband.
The fourth proband, female, 50 years, was diagnosed with swelling of the lower limbs at age 42. She carries a heterozygous single-nucleotide missense variant c.1000C > T (dbSNP rs114144673), with a frequency of 0.001525 (GnomAD). This is a sporadic case without family history, and no other family members were tested.

In Silico Analysis, Template Selection and Building
In an attempt to study the possible role of the variants found in our cohort compared to the wildtype, we performed an in-silico study. A wildtype template (Table 2) search with BLAST and HHBlits against the SWISS-MODEL template library (last update: 2019-10-24; last included PDB release: 2019-10-18) produced a total of 267 and 254 templates that matched the NRP1 and NRP2 genes, respectively, with different sequence identity and quality percentages. Details of the top 10 templates are shown in Table 4.
Based on the percentage of sequence identity, similarity and best quality square, the 4gz9.1.A and 2qql.1.A chains were selected to align the template and query sequences in order to build models of NRP1 and NRP2 (Figures 4 and 5). Based on the percentage of sequence identity, similarity and best quality square, the 4gz9.1.A and 2qql.1.A chains were selected to align the template and query sequences in order to build models of NRP1 and NRP2 (Figures 4 and 5).      Then, we entered the model in the Discovery studio visualizer to generate versions of the NRP1 structure with Arg552Gln and versions of the NRP2 structure with Phe194Val, Pro280Ser, Arg334Cys and Ile583Thr. Unfortunately, the template to generate the model for Ile853 was not available, so we could not study its interaction with adjacent residues. A molecular-level interaction analysis between native/variant residues was performed ( Figures 6-10 show snapshots). Details of the residues involved in the interactions, the types of bond they formed and bond lengths in angstrom units are shown in Tables 5 and 6. Then, we entered the model in the Discovery studio visualizer to generate versions of the NRP1 structure with Arg552Gln and versions of the NRP2 structure with Phe194Val, Pro280Ser, Arg334Cys and Ile583Thr. Unfortunately, the template to generate the model for Ile853 was not available, so we could not study its interaction with adjacent residues. A molecular-level interaction analysis between native/variant residues was performed (Figures 6-10 show snapshots). Details of the residues involved in the interactions, the types of bond they formed and bond lengths in angstrom units are shown in Tables 5 and 6.
Briefly, for all the residues analyzed, there is a change in the number or length of the bonds that are formed in the variant with respect to the modeled wildtype structure. Indeed, the in-silico analysis showed that NRP1 encoded with Arg552 has slightly different stability from the Gln552 variant, because it has seven interactions with adjacent residues, whereas the variant has only three that exactly match three of the Arg552 interactions (Table 5 and Figure 6).  By contrast, NRP2 encoded with Phe194 has slightly different stability from the Val194 variant due to four interactions with adjacent residues, whereas the variant has only two, similar to Phe194 but slightly different in bond length (Table 6 and Figure 7).  In the case of Pro280Ser, there is a shift from a nonpolar residue (Pro) to a polar residue (Ser). Differences in the amino acid side chain greatly alter the interaction of the residues with the phenylalanine at residue 425, leading to a loss of Pi interaction. These results suggest that the overall protein structure is altered by these different interactions with nearby residues and is functionally defective ( Table 6 and Figure 8). In the case of Arg334Cys, Arg334 has six interactions with adjacent residues, whereas Cys334 has only three similar to those of Arg334, and the interaction with Ile226 has a different hydrophobic interaction bond length (Table 6 and Figure 9). In the case of Pro280Ser, there is a shift from a nonpolar residue (Pro) to a polar residue (Ser). Differences in the amino acid side chain greatly alter the interaction of the residues with the phenylalanine at residue 425, leading to a loss of Pi interaction. These results suggest that the overall protein structure is altered by these different interactions with nearby residues and is functionally defective ( Table 6 and Figure 8). In the case of Arg334Cys, Arg334 has six interactions with adjacent residues, whereas Cys334 has only three similar to those of Arg334, and the interaction with Ile226 has a different hydrophobic interaction bond length (Table 6 and Figure 9).  In the case of Pro280Ser, there is a shift from a nonpolar residue (Pro) to a polar residue (Ser). Differences in the amino acid side chain greatly alter the interaction of the residues with the phenylalanine at residue 425, leading to a loss of Pi interaction. These results suggest that the overall protein structure is altered by these different interactions with nearby residues and is functionally defective ( Table 6 and Figure 8). In the case of Arg334Cys, Arg334 has six interactions with adjacent residues, whereas Cys334 has only three similar to those of Arg334, and the interaction with Ile226 has a different hydrophobic interaction bond length (Table 6 and Figure 9).  Likewise, in the case of Ile583Thr, Ile583 has six interactions with adjacent residues, whereas the Thr583 variant has only four that are the same, but the Pi interaction with Trp578 varies slightly in bond length (Table 6 and Figure 10).

Discussion
Lymphedema is a progressive disease caused by lymphatic system malformations and impaired lymph flow. A panel consisting of 29 known lymphedema genes is currently used for the genetic testing of patients [24] but only detects relevant variants in 25-30% of patients [35,36]. This highlights the need for new target genes to screen. Here, we report our findings from screening a cohort of 235 patients, negative for variants in the said 29 genes, for variants in the neuropilin-coding genes NRP1 and NRP2.
Neuropilins are highly conserved transmembrane receptors that perform various functions in vessels, neurons and tumors, because they bind structurally to different ligands and coreceptors. Although NRP1 and NRP2 share basic structural features, they are implicated in different biological processes, such as the development of blood and lymphatic endothelial cells [37]. They have a role in the development of lymphatic valves [13,17,18], lymphatic vessels and capillaries [12,19,22].
Overall, our results showed rare variants in NRP1 and NRP2 in patients who tested negative for variants in the 29 genes already associated with lymphedema. In three out of six cases, family members were screened, highlighting segregation of the variant with lymphedema or with subclinical signs of abnormal lymph flow. Bioinformatic modeling added evidence of a possible impact of the variant on the stability of the wildtype protein. Further biochemical studying or in-vitro studying is necessary to verify the altered function of coreceptors carrying the variant. In addition, other studies aiming at evaluating the effect of the interaction of the altered NRP1 with its principal ligand SEMA3A, or of the altered NRP2 with its coreceptor VEGFR3, could be important to define whether the variants also have a role in such signaling. Moreover, the loss of SEMA3A function by the specific inhibition of its interaction with NRP1 leads to impairment of lymphatic valve development [12], while variants in VEGFR3, a tyrosine kinase receptor important in the development and establishment of the lymphatic system, are known to cause Milroy disease, a form of primary lymphedema [38].
Additional studies could shed light on altered pathways, leading to a better understanding of the mechanism of action of the alteration; in particular, it could help determine whether the phenotype seen in our patients is due to haploinsufficiency or to a dominant negative effect of the coreceptor carrying the variant.  Briefly, for all the residues analyzed, there is a change in the number or length of the bonds that are formed in the variant with respect to the modeled wildtype structure. Indeed, the in-silico analysis showed that NRP1 encoded with Arg552 has slightly different stability from the Gln552 variant, because it has seven interactions with adjacent residues, whereas the variant has only three that exactly match three of the Arg552 interactions (Table 5 and Figure 6).
By contrast, NRP2 encoded with Phe194 has slightly different stability from the Val194 variant due to four interactions with adjacent residues, whereas the variant has only two, similar to Phe194 but slightly different in bond length (Table 6 and Figure 7).
In the case of Pro280Ser, there is a shift from a nonpolar residue (Pro) to a polar residue (Ser). Differences in the amino acid side chain greatly alter the interaction of the residues with the phenylalanine at residue 425, leading to a loss of Pi interaction. These results suggest that the overall protein structure is altered by these different interactions with nearby residues and is functionally defective ( Table 6 and Figure 8).
In the case of Arg334Cys, Arg334 has six interactions with adjacent residues, whereas Cys334 has only three similar to those of Arg334, and the interaction with Ile226 has a different hydrophobic interaction bond length (Table 6 and Figure 9).
Likewise, in the case of Ile583Thr, Ile583 has six interactions with adjacent residues, whereas the Thr583 variant has only four that are the same, but the Pi interaction with Trp578 varies slightly in bond length (Table 6 and Figure 10).

Discussion
Lymphedema is a progressive disease caused by lymphatic system malformations and impaired lymph flow. A panel consisting of 29 known lymphedema genes is currently used for the genetic testing of patients [24] but only detects relevant variants in 25-30% of patients [35,36]. This highlights the need for new target genes to screen. Here, we report our findings from screening a cohort of 235 patients, negative for variants in the said 29 genes, for variants in the neuropilin-coding genes NRP1 and NRP2.
Neuropilins are highly conserved transmembrane receptors that perform various functions in vessels, neurons and tumors, because they bind structurally to different ligands and coreceptors. Although NRP1 and NRP2 share basic structural features, they are implicated in different biological processes, such as the development of blood and lymphatic endothelial cells [37]. They have a role in the development of lymphatic valves [13,17,18], lymphatic vessels and capillaries [12,19,22].
Overall, our results showed rare variants in NRP1 and NRP2 in patients who tested negative for variants in the 29 genes already associated with lymphedema. In three out of six cases, family members were screened, highlighting segregation of the variant with lymphedema or with subclinical signs of abnormal lymph flow. Bioinformatic modeling added evidence of a possible impact of the variant on the stability of the wildtype protein. Further biochemical studying or in-vitro studying is necessary to verify the altered function of coreceptors carrying the variant. In addition, other studies aiming at evaluating the effect of the interaction of the altered NRP1 with its principal ligand SEMA3A, or of the altered NRP2 with its coreceptor VEGFR3, could be important to define whether the variants also have a role in such signaling. Moreover, the loss of SEMA3A function by the specific inhibition of its interaction with NRP1 leads to impairment of lymphatic valve development [12], while variants in VEGFR3, a tyrosine kinase receptor important in the development and establishment of the lymphatic system, are known to cause Milroy disease, a form of primary lymphedema [38].
Additional studies could shed light on altered pathways, leading to a better understanding of the mechanism of action of the alteration; in particular, it could help determine whether the phenotype seen in our patients is due to haploinsufficiency or to a dominant negative effect of the coreceptor carrying the variant.
Although these results are not sufficient to associate the two genes, known to be involved in the development of lymphatic capillaries and valves, with lymphedema, they suggest that it would be worthwhile screening a larger cohort of patients for variants in NRP1 and NRP2, especially in large families with more than one affected subject, and performing a further functional analysis of the variants identified by us in in vitro and in vivo models to confirm their involvement in the development of lymphedema.

Conclusions
A comprehensive review of the literature strongly suggests that neuropilins NRP1 and NRP2 play a role in the development of lymphatic capillaries and lymphatic valves, crucial components of the lymphatic system. If these components do not function normally, the lymph does not drain, and fluid may accumulate in tissues, eventually leading to lymphedema. Given their functional role and our findings, we suggest that NRP1 and NRP2 should be considered candidate genes for inclusion in the gene panel for genetic testing of lymphedema patients.