A combined evolutionary and structural approach to disclose the primary structural determinants essential for proneurotrophins biological functions.

The neurotrophins, i.e., Nerve Growth Factor (NGF), Brain-Derived Neurotrophic Factor (BDNF), Neurotrophin 3 (NT3) and Neurotrophin 4 (NT4), are known to play a range of crucial functions in the developing and adult peripheral and central nervous systems. Initially synthesized as precursors, i.e., proneurotrophins (proNTs), that are cleaved to release C-terminal mature forms, they act through two types of receptors, the specific Trk receptors (Tropomyosin-related kinases) and the pan-neurotrophin receptor p75NTR, to initiate survival and differentiative responses. Recently, all the proNTs but proNT4 have been demonstrated to be not just inactive precursors, but signaling ligands that mediate opposing actions in fundamental aspects of the nervous system with respect to the mature counterparts through dual-receptor complexes formation with a member of the VPS10 family and p75NTR. Despite the functional relevance, the molecular determinants underpinning the interactions between the pro-domains and their receptors are still elusive probably due to their intrinsically disordered nature. Here we present an evolutionary approach coupled to an experimental study aiming to uncover the structural and dynamical basis of the biological function displayed by proNGF, proBDNF and proNT3 but missing in proNT4. A bioinformatic analysis allowed to elucidate the functional adaptability of the proNTs family in vertebrates, identifying conserved key structural features. The combined biochemical and SAXS experiments shed lights on the structure and dynamic behavior of the human proNTs in solution, giving insights on the evolutionary conserved structural motifs, essential for the multifaceted roles of proNTs in physiological as well as in pathological contexts.

Introduction

Neurotrophins (NTs) are known to play differential roles in the developing and adult nervous system and in non-neuronal organs by promoting differentiation, neuronal survival, synaptogenesis and modulating synaptic plasticity [1]. Only four NT genes are known in mammals i.e., Nerve Growth Factor (NGF), Brain-Derived Neurotrophic Factor (BDNF), Neurotrophin 3 (NT3) and Neurotrophin 4 (NT4). They exert their biological functions through two types of receptors, the specific Trk tyrosine kinases (Tropomyosin-related kinases) and the low-affinity pan-neurotrophin receptor p75NTR, which can bind all NTs with similar affinity [2], [3].

NTs are initially synthesized as larger precursor proteins, proneurotrophins (proNTs), consisting of a N-terminal pro-domain and a C-terminal mature domain [4]. These precursor proteins, rapidly associate as noncovalent homodimers and can enzymatically be cleaved by means of intracellular proteases (furin and proconvertases, plasmin, tryptase and specific matrix metallo proteases [5]), in the Golgi and endoplasmatic reticulum to produce C-terminal mature NTs [6].

ProNTs were historically considered to be just inactive precursors, with the pro-domains only acting as an intramolecular chaperone and facilitating the folding of the mature domain and regulating the respective secretory pathways. Recent findings demonstrated that all proNTs but proNT4 can be secreted to the extracellular space and are biologically active ligands [7], by exerting opposing actions in fundamental aspects of the nervous system (such as neuronal survival and synaptic plasticity) with respect to the mature counterparts [8]. ProNTs induce activation, at sub-nanomolar concentrations, of the apoptotic machinery [9], [10], [11] with subsequent cell death of specific neuronal populations; regulate axon guidance in neuronal development [12] and mediate synaptic plasticity, namely Long-Term Depression (LTD), in hippocampal neurons [13]. Besides their physiological apoptotic function, proNTs also play an important role in neuropathology. Indeed, apoptotic pathways activated by proNTs have been proven to be induced after injury and in several neurodegenerative disorders by promoting proNTs secretion to the extracellular space [14]. Likewise NTs, proNTs utilize dual-receptor complexes to mediate their functions. They form a ternary assembly composed of either a member of the VPS10 family, i.e., Sortilin, Sortilin related VPS10 Domain Containing Receptor 2 (SorCS2) or Receptor 3 (SorCS3) and p75NTR [9], [15], [16], to regulate synaptic plasticity and trigger apoptotic signaling. Therefore, the outcome of NT action depends on the kind of ligand that is secreted by the cells (proNT or mature NT), as well as on the interaction with distinct interplaying receptor complexes. Indeed, both proNTs and NTs are able to interact with VPS10p members, p75NTR and Trk, but in general proNTs preferentially bind to VPS10p members while mature NTs have higher affinity for p75NTR and Trk [9], [12]. The only exception is NGF that interacts with SorCS3 tighter than proNGF [17]. Co-expression of SorCS2 (or Sortilin) and p75NTR allows simultaneous binding to proNGF, substantially enhancing the binding affinity and leading to the ternary complex formation required for signaling [9], [10], [12], [15].

Crystallographic studies disclosed NGF and proNGF to be engaged with p75NTR via the same binding site on the mature domain of NGF with the only repositioning of a single loop [18], [19], [20]. Similarly, the interactions in the proNGF-SorC2 complex [21] occur via the mature domain. The binding site overlaps those present in the complexes with p75NTR and Trk, respectively. It is worth noting that the other member of VPS10p family, Sortilin, has been shown to specifically recognize the pro-domains of the three proNTs (proNGF, proBDNF, and proNT3) in the formation of the co-receptor complex with p75NTR. In this context, the structural and molecular determinants underpinning the interactions between the pro-domains and the proNTs receptors are still elusive. Indeed, in the previously discussed crystal structure of the complex between p75NTR and proNGF [18] the electron density corresponding to the entire pro-domain was poorly traceable, precluding to predict where it would interact on the receptor. This is likely ascribable to the intrinsically disordered nature of the pro-domain of proNTs thoroughly assessed for hproNGF and in part for hproBDNF by CD, FT-IR, fluorescence, chemical and thermal denaturation, H/D exchange, NMR and limited proteolysis studies [22], [23], [24], [25], [26], [27], [28], [29], [30], [31].

Pre-computed disorder predictions of the four hproNTs as well as their disorder-related functions are available by searching the curated IDP database D2P2 [32].

The main feature of intrinsically disordered proteins (IDPs) and intrinsically disordered regions (IDRs) is the partial or complete lack of a cooperative, folded structure under native conditions, due to the typical composition of their primary sequence enriched in charged, polar and structure-breaking amino acids and lacking bulky hydrophobic amino acids [33].

Consequently, under physiological conditions, every IDP molecule samples a range of conformations over time. Some conformers of the same molecule might resemble a random coil structure, while others display structured elements, even if the dynamic nature of IDPs and the inherent ensemble averaging of in-bulk investigations can make sparsely populated conformations difficult to be detected by traditional biophysical techniques. These so-called transient secondary structures can act as seeds that promote the formation of local specific tertiary conformations upon binding to their partners in a coupled folding and binding process [33], [34], [35], [36]. Binding through unfolded or partially unfolded intermediates, though controversial, can provide a kinetic advantage through the ‘‘fly-casting’’ mechanism [37], [38]. According to this mechanism, a dimensionality reduction occurs when the folding of a disordered protein is coupled with binding, thereby speeding up the search for specific targets. These molecular interactions can be not only transient and dynamic, but also permanent, enabling IDPs to modify the activity of the partner or to act as protein hubs in intracellular networks. Furthermore, the folding of IDPs upon contact with partners can also be incomplete or even absent in some cases (according to the so-called “fuzziness” phenomenon) [39]. In this framework, disorder in protein can be static, if the IDP adopts a few or a multitude of distinct stable conformations (“polymorphic” model), or dynamic, in the case of the regions that are still unstructured in the complex either neighboring, i.e., “flanking” model, or linking, i.e., “clamp” model, ordered binding regions. If the complex formation does not induce folding of the IDPs, then the bound state is left entirely disordered (“random” model). These specific features determine the functional promiscuity of IDPs under different conditions.

Furthermore, the analysis of structured states of apparently disordered proteins from genomic sequences highlighted the co-occurrences of evolutionary constraints of functional importance [40].

Interestingly, IDPs generally appear to evolve more rapidly, to be more permissive to different point mutations and to have higher amino acids deletion and insertion rates compared to ordered proteins [41]. In this regard, proteins can display different modes of how sequence conservation can relate to disorder as a functional feature [42]. While in some cases the IDPs’ disorder itself is conserved, but not the sequence facilitating it (i.e., flexible disorder), other IDPs retain disorder together with a highly conserved sequence (i.e., constrained disorder). Thus, it is of primary interest to integrate structural and evolutionary analysis in order to reveal common mechanisms used by members of the proNTs family to deal with these fast rates of evolution and to gain insights into the structural and molecular determinants underpinning their functional adaptability.

Structural characterization of IDPs by conventional methods as X-ray crystallography and cryo-electron microscopy is rather challenging due to their extremely high flexibility in solution. Indeed, X-ray crystallography is not applicable because not only the intrinsic conformational heterogeneity of this class of proteins usually hampers crystallization, but also because, in the crystal structures of proteins containing IDRs, the corresponding electron density maps are poorly resolved. Therefore, high-resolution structural characterization of IDPs mainly relies on Nuclear Magnetic Resonance (NMR) spectroscopy that is especially well suited to probe structurally flexible systems and dynamics [43], [44]. Unfortunately, the size of the functional proNTs homodimers hampers to exploit this technique.

Small Angle X-ray Scattering (SAXS) represents instead the experimental technique of choice to deal with the structural flexibility of IDRs [45]. It has already been proven successful to structurally characterize a member of the proNTs family, i.e., mouse proNGF [26], [46]. This approach allows probing the size and shape of unstructured proteins and their complexes with receptors. Although its intrinsic low-resolution nature, SAXS analysis allows to quantitatively describe systems characterized by large-scale protein fluctuations and the presence of multiple species and/or conformations in solution [47], [48], [49], employing ensembles of conformers to describe the experimental data [50].

Here we present a computational biology approach coupled to an experimental study by SAXS aiming to uncover the structural and dynamical basis of the biological function displayed by proNGF, proBDNF and proNT3 but missing in proNT4. An in silico bioinformatic analysis, following an evolutionary approach recently implemented [51], allowed to elucidate the functional adaptability and to identify key structural features of the members of the proNTs family, by predicting their structural disorder and secondary structure propensities and computing the rates of evolution per site and the rates of structural properties transitions.

A combined experimental approach was then exploited to shed lights on the structure and the dynamic behavior of the human members of proNTs family in solution. We report on complementary biochemical and SAXS experiments aiming to uncover the molecular mechanisms underpinning the biological functions conserved in proNGF, proBDNF and proNT3, but lost in proNT4 during evolution.

Materials and methods

Bioinformatics

The intrinsic disorder property within the amino acid sequence of the pro-domains of the four hproNTs, was characterized by using a set of commonly per-residue disorder predictors, such as MFDp2 [52], the IUPred suite for the prediction of short and long intrinsically disordered regions, respectively [53], PONDR VLXT [54], PONDR VSL2 [55], PONDR VL3 [56] and PONDR FIT [57]. Similarly, the results of different secondary structure predictors (i.e., Psipred [58], Jpred4 [59], PredictProtein [60] and RaptorX SS3 [61]) were compared by analyzing the amino acid sequences of the pro-domains of the four hproNTs.

Vertebrate sequences were collected by searching the NCBI database, using NCBI BLAST [62] by the Blastp (protein–protein BLAST) algorithm. The protein sequences of the four human proNTs, i.e., hproNGF: P01138; hproBDNF: P23560; hproNT3: P20783; hproNT4: P34130 [63] were used as query. A filter for high sequence coverage was applied and for each of the species only the sequences with the highest BLAST Max score were selected to minimize redundancy. In the case of high sequences abundance of a single taxon, i.e., Mammalia and Aves, sequences have been further filtered in order to equally cover all the orders of each taxa. For Multiple Sequence Alignment (MSA), datasets of the whole vertebrate sequences and of the representative organism were constructed for each of the proNTs. The species included in the representative datasets were: Homo sapiens, Pongo abelii, Ornithorhynchus anatinus, Monodelphis domestica, Mus musculus, Columba livia and Gavia stellata (for proNGF, proBDNF and proNT3); Gekko japonicus, Podarcis muralis, Nanorana parkeri, Rhinatrema bivittatum, Amblyraja radiata and Rhincodon typus (for proBDNF and proNT3 and proNT4); Latimeria chalumnae and Paramormyrops kingsleyae (for proBDNF and proNT3 and proNT4); and Anarrhichthys ocellatus (for proNGF, proNT3 and proNT4). All datasets were aligned using Clustal Omega [64] and the tree building tool was employed to construct the corresponding phylogenetic trees. Structural disorder propensity of all unaligned vertebrate datasets was predicted using MFDp2 (Multilayered Fusion-based Disorder predictor v. 2.00), which allows to process simultaneously up to 100 protein sequences [52]. It combines per-residue disorder probabilities predicted by MFDp with per-sequence disorder content predicted by DisCon [65]. Novel post-processing filters are applied to provide disorder predictions with improved predictive quality. The order–disorder propensity scores were mapped onto the MSA and taxon position in the tree to be visualized as a heat map and were further converted into a binary matrix to analyze the evolutionary dynamics of structural order–disorder (i.e., 1 for order and 0 for disorder). Secondary structure propensity was predicted for all sequences in each dataset using Jpred 4 [59], the latest version of the popular JPred protein secondary structure prediction server, which provides predictions by the JNet algorithm, one of the most accurate methods for secondary structure prediction, being its secondary structure prediction accuracy reported to be 82%. As for the analysis of the order–disorder propensity, the results of the secondary structure propensity were mapped to be visualized in a heat map and then converted into a binary matrix (i.e., 1 for helix and strand, and 0 for loop), to examine the evolutionary dynamics of secondary structure-loop transitions. The substitution rates of amino acids for each vertebrate dataset alignment were analyzed using ConSurf [66], a bioinformatics tool that applies statistical inference methods, machine learning and multiple sequence alignment to attribute a conservation score to each residue based on the phylogenetic relationship between the protein and its homologous sequences.

The evolutionary dynamics of structural order–disorder and predicted secondary structures were analyzed by calculating the evolution rates based on the phylogenetic trees and the binary matrices. The Gain Loss Mapping Engine (GLOOME) [67] was used in order to study the gain/loss transition events of the structural property by adopting the phyletic patterns. The gaps were treated as missing data and the outputs of evolution rates were standardized as Z-score.

DNA constructs

The cDNA encoding for human BDNF (hBDNF), human NT3 (hNT3) and human NT4 (hNT4) were purchased from Harvard Medical School PlasmID Database (clone_Ids: HsCD00005889, HsCD00001730 and HsCD00000988) corresponding to Homo sapiens brain-derived neurotrophic factor (NCBI Gene ID: 627), Homo sapiens neurotrophin 3 (NCBI Gene ID: 4908) and Homo sapiens neurotrophin 4 (NCBI Gene ID: 4909), respectively.

hBDNF, hNT3 and hNT4 in fusion with the corresponding N-terminal pro-domains (hproBDNF, hproNT3 and hproNT4) without the secretion leader sequence, corresponding to codons 57–741 of hpreproBDNF cDNA, 57–771 of hpreproNT3 cDNA and 75–630 of hpreproNT4 cDNA were amplified by PCR using the following forward and reverse primers, respectively:

5′- CTTTAAGAAGGAGATATACATATGGCCCCCATGAAAGAAGCAAACATC

3′-GGCTTTGTTAGCAGCCGGATCTCATCTTCCCCTTTTAATGGTCAATGTAC

for hproBDNF

5′- CTTTAAGAAGGAGATATACATATGAACAACATGGATCAAAGGAGTTTG

3′- GGCTTTGTTAGCAGCCGGATCTCATGTTCTTCCGATTTTTCTCGAC

for hproNT3

5′- CTTTAAGAAGGAGATATACATATGCAACCCCCACCCTCAACATTGC

3′- GGCTTTGTTAGCAGCCGGATCTCAGGCCCGGCCAGTCCG

for hproNT4.

They were cloned in pET22b vector by Restriction free Cloning [68], removing the bacterial secretion leader, the polylinker and the C-terminal His Tag. After screening by colony PCR, the positive clones were verified by automated DNA sequencing.

All the sequenced clones for hproBDNF carried the V66M polymorphism that was indeed present in the purchased plasmid. The sequence encoding the wild type protein was generated from the hproBDNF V66M construct by retro mutation (M66V) according the Stratagene site-directed mutagenesis protocol.

The following forward and reverse primers were employed respectively:

5′-GCTGACACTTTCGAACACGTGATAGAAGAGCTGTTGG

3′- CCAACAGCTCTTCTATCACGTGTTCGAAAGTGTCAGC.

The construct was verified by automated DNA sequencing.

Protein expression, refolding and purification

Human proNTs expression and purification were performed based on the method previously described for hproNGF [69]. pET22b constructs encoding for the proNTs were transformed into Rosetta (DE3) strain (Stratagene) and protein expression was induced at OD(600 nm) = 1.0 by 1 mM IPTG overnight at 30 °C. Pulsed refolding of hproNT3 was performed as previously reported for hproNGF, while in the case of hproBDNF and hproNT4 the refolding buffer was further optimized and 750 mM Arginine, 100 mM Tris pH 8.5, 5 mM EDTA, 5 mM reduced glutathione, 2 mM oxidized glutathione, 0.05% (v/v) PMSF was employed. After extensive dialysis against 50 mM Sodium Phosphate pH 7.0, 1 mM EDTA and 0.05% (v/v) PMSF at 4 °C, the proteins were purified by cationic exchange chromatography on HiTrap SP column (GE Healthcare) and concentrated by stirred ultrafiltration cells (Millipore).

Differential Scanning Fluorimetry (DSF)

DSF experiments were performed using a CFX96 Touch Biorad Real-Time PCR system (Bio-Rad) with λex = 470–505 nm and λem = 540–700 nm. The proNTs stocks in 50 mM Sodium Phosphate buffer pH 7.0 were mixed at final concentration of 0.5 mg/mL with SYPRO Orange (Sigma) at final concentration of 90×. Fluorescence was measured in a 20°–90 °C temperature range with increments of 0.2 °C/min. The melting temperatures (Tm), represented by the inflection points of the transition curves, were calculated using Boltzmann Sigmoid fit. Experiments were performed in triplicate.

Small-angle X-ray scattering data collection and data processing

SAXS measurements on 5 different concentrations (the ranges are reported in Table S1) of the four hproNTs in 50 mM Sodium Phosphate pH 7.0 were performed on the P12 beamline EMBL SAXS-WAXS (Hamburg, Germany at PETRAIII/DESY) [70]. Data were collected with 20x 0.05 sec exposures using a Pilatus 2 M pixel X-ray detector, sample-detector distance 3.00 m, wavelength 1.24 Å. No radiation damage was detected, comparing scattering profiles for the collected frames.

Data were processed employing the ATSAS 2.6.0 program package [71]. In detail, normalization to the intensity of the transmitted beam, averaging of the frames for each sample and buffer subtraction were performed by PRIMUS [72]. Since for the hproBDNF sample some inter-particle repulsion effect was observed at increased protein concentration, the low s-data of diluted samples (0.41 mg/ml), where inter-particle interactions did not occur, were merged with the high s-data of the concentrated samples.

Guinier approximation [73] was applied to calculate the radius of gyration Rg and the forward scattering I(0), assuming that at very small angles (0 < s < 1.3/Rg), the intensity I(s) = I(0) exp(−1/3(Rg·s)2). GNOM [74] was used to compute the maximum sizes Dmax and the pair distance distribution functions p(r). Molecular masses (MM) were estimated by comparison of the calculated forward scattering I(0) of the samples with that of the standard solution of bovine serum albumin (MM 66 kDa). The Porod approximation [75]:was used to calculate the excluded volume of the hydrated protein molecule (Vp).

The web-based server tool SWISS-MODEL [76] was exploited to build a 3D homology model of mature hBDNF dimer using the high-resolution X-ray crystal structure of the hBDNF/hNT3 heterodimer (PDB_ID 1BND) [77] as a template. The homodimeric assembly was subjected to 4000 steps of steepest-descent in vacuo energy minimization, followed by 10,000 steps of conjugate gradient in vacuo energy minimization using GROMACS 4.5 [78].

Rigid-body modeling was performed by using CORAL [79] in order to generate the approximate conformations of the missing regions of the propeptides. The crystallographic structures of hNT3 (PDB_ID 1NT3) [80] and hNT4 (PDB_ID 1HCF) [81] and of the obtained high resolution 3D models of hNGF [82] and hBDNF, have been kept fixed and either P1 or P2 symmetry have been imposed. Ten independent CORAL runs were performed for each data set. The resulting models were superimposed by using the program SUPCOMB [83] and averaged by using DAMAVER [84] to identify the most typical models representing the human proNTs conformation in solution. The similarity of the resulting models was estimated by DAMAVER calculating the normalized spatial discrepancy parameter (NSD) [83]. NSD values ≤ 1.0 are expected for similar models.

The Ensemble Optimization Method (EOM) [85] allowed to quantitatively assess propeptide flexibility and size distribution of possible conformers. Randomized 10,000 conformations were initially generated for the hproNTs based on amino acid sequences and the crystallographic structures of hNT3 and hNT4 and the 3D-structural models of hNGF and hBDNF, used as rigid bodies, imposing the native conformation option and applying P2 symmetry. The scattering profiles of these pools of randomly generated conformations, computed by CRYSOL [86], were compared and representative structures, whose scattering curve fits to the experimental scattering curve, were selected by a genetic algorithm.

Results

The bioinformatic analysis of the evolution of structural order-disorder and secondary structures propensities highlights a divergent behavior of proNT4 with respect to the other members of the proNTs family

The comparison of the pro-domains among the human members of the proNTs family highlights a high variability both in their length and in their sequences (Fig. 1a), being the pro-domain of hproNT4 the most divergent.

Fig. 1

Bioinformatic comparison of the primary sequences of the hproNTs: (a) multiple alignment of the primary sequences of the corresponding pro-domains. Each residue is colored according the conservation scores based on the phylogenetic relationship between the pro-domains of hNTs and their homologous sequences calculated by ConSurf [66]. Above each sequence, the secondary structure propensity based on the Jpred4 (Drozdetskiy et al. 2015) prediction score (i.e., red for helix and blue for strand) is reported. (b) MFDp2 predictions of intrinsically disordered regions [59]; below each graph the schematic representation of the corresponding amino-acid sequence is shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Indeed, its length is almost half of the other members and it shares a lower sequence similarity. It displays a lower mean percentage of pair sequence identity (i.e., 17.7 ± 2.0% instead of 23.7 ± 3.6%) and almost halved mean percentage of pair sequence homology (i.e., 24.2 ± 1.4% instead of 43.0 ± 3.1%) with respect to the others hproNTs. The mature domains of hproNTs, instead, show a lower degree of variability. hNT4 displays very similar and much higher mean percentage of pair sequence identity (i.e., 52.2 ± 2.6%) and homology (i.e., 73.2 ± 3.7%) if compared to the others hNTs (55.0 ± 2.5% and 71.8 ± 4.4%, respectively). In analyzing the sequence similarity within this family, it deserves attention a main structural feature, at present reported only for proNGF and proBDNF [26], [29], i.e., the intrinsic disorder of the pro-domains. It is known that IDRs evolve more rapidly than globular folded domains. The significant differences in sequence conservation between the pro-domains and the mature regions within hproNT family prompt to the intrinsic disorder to be a common feature among all the pro-domains. Thus, by taking into account that disorder is a fundamental property of the proteome, robustly predictable from primary sequence relying on characteristic patterns of amino acids distribution and overall amino acids content [87], we extended the analysis of the disorder content to the other two members of the proNTs family (Fig. 1b). We assessed a set of commonly used per-residue disorder predictors, selected on their specific sensitivities to different features associated with intrinsic disorder (Supp. Fig. S1). Almost all the predictors concur with the reported intrinsically disordered nature of the pro-domain of hproBDNF and point to a very low disorder content in the pro-domain of hproNT4, mainly confined at the first 10 residues. On the contrary, a wide variability in the degree of disorder has been observed for the pro-domains of hproNGF and hproNT3 (Supp. Fig. S1a and c). Considering that the intrisic disorder of the pro-domain of hproNGF and hproBDNF has been experimentally characterized [26], [29], [31], the most accurate predictor resulted to be MFDp2 [52], (Fig. 1B), which allows to quantitatively predict the disorder content by means of DisCon [65], a predictor that takes advantage of a set of descriptors that aggregate and hybridize information originated from sequence, evolutionary profiles, and secondary structure prediction, solvent accessibility, flexibility, and annotation of globular domains. Indeed, while the pro-domain of hproNT3 is estimated to be mainly disordered having a higher disorder content than the other two biologically active hproNTs (being the DisCon index of 0.379, 0.297 and 0.414 for hproNGF, hproBDNF and hproNT3, respectively), hproNT4 is characterized by a significantly lower disorder content (being the DisCon index of 0.104). Moreover, the comparison of the evolutionary conservation of the amino acid sequences (Fig. 1a) suggests that the pro-domains of hproNGF, hproBDNF and hproNT3 are characterized by constrained disorder, especially at the level of the regions that are predicted to be involved in secondary structure formation, according to the consensus obtained employing three different secondary structure predictors (Supp. Fig. S2).

To this regard, despite their designation as “disordered”, IDPs do not generally behave as structure-less random polymer chains and frequently display residual secondary structures. Indeed, although the high degree of intrinsic disorder of the pro-domains, hproNGF, hproBDNF and hproNT3 are characterized by both the presence of helical content, that is absent in hproNT4, and by a conserved short beta strand that, in turn, iscommon to the pro-domains of all the members of the family.

In order to further investigate the relationship between rates of amino acid substitution in the pro-domains and their structural and functional properties across evolution, we calculated the conservation Z-scores for each position, being positives the scores for those sites evolving faster than average and negative the scores for those sites evolving slower than average (upper graphs, Fig. 2).

Fig. 2

Rates of evolution per site of hproNTs and evolutionary dynamics of structural properties: a) hproNGF, b) hproBDNF, c) hproNT3 and d) hproNT4. Two graphs are shown in each panel: the rate of evolution per site (top) and the evolutionary dynamics of predicted order–disorder and secondary structures (bottom). The secondary structure propensity bar (middle) indicates conserved helix and strand structures as red and blue box, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Moreover, to elucidate the changes in protein dynamics/structure of the pro-domains that are likely to be highly related to their biological function, we focused on the evolutionary patterns of structural properties of the proNTs members in vertebrates. The heat maps of secondary structure propensity and order–disorder structure propensity (Supp. Figs. S3 and S4) have been mapped to the MSA, ordered by the taxa position in the phylogenetic trees. The dynamic evolution of those structural properties has been analyzed by calculating the corresponding transition rates (lower graphs, Fig. 2). Despite each of the pro-domains is characterized by a specific evolutionary pattern of the structural properties, common features as well as diverging behaviors can be clearly identified.

Although the four pro-domains are very divergent both in length and in sequence, the secondary structure prediction results showed that all the pro-domains share a highly conserved beta strand (Supp. Fig. S3). The fact that this short secondary structure element, ranging between 2 and 4 residues, is extremely conserved across evolution with very similar transition patterns and rates (Fig. 2) in the proNTs superfamily, might suggest that it could act as an essential seed, mediating pro-domain biological functions associated with NT folding and maturation. An additional conserved secondary structure helix element is found in all the pro-domains but not in proNT4. Besides showing a highly conserved predicted long helical structural element (ranging 12 to 25 residues), located close to the amino-terminal, proNT3 shares with proNGF an homolog region (ranging 3 to 10 amino acids) having helix propensity, which is highly conserved across evolution. A small, but distinct, difference in the transition patterns and rates (Supp. Fig. S3a and c and Fig. 2a and c, respectively) has been detected. In addition, proBDNF shows conserved helix secondary structure propensity mapped in two regions (Supp. Fig. S3b), whose transition patterns and rates (Fig. 2b) are very similar to the helical structure shared by proNGF and proNT3. Since Sortilin recognizes a binding site encompassing amino acids 44–103 of proBDNF [88], the amino-terminal alpha helix element might then be regarded as the structural seed involved in this interaction, despite differences in the amino acid composition if compared to the corresponding alpha helical segment shared by proNGF and proNT3. Nevertheless, a typical hallmark of IDPs is their tolerance to mutations across evolution of paralogs allowing the biological function to be preserved. It is known that proBDNF derives from a different intermediate neurotrophin ancestral gene with respect to the one that originated proNGF and proNT3 [89], so explaining the basis for the divergence of the primary sequence with respect to the other biologically active proNTs.

Interestingly, the pro-domain of proNT4 is characterized by an apparent distinct behavior with respect to all the other members. Indeed, no conserved pattern of helical secondary structure could be identified in proNT4. A helical structure propensity is only present in a number of mammalian orders, but is almost completely absent in primates (Supp. Fig. S3d). This divergent behavior of proNT4 from the other members of the family might suggest that the very conserved helical structure could be involved in proNT specific biological functions associated with cell death and LTD, which on the contrary are not mediated by proNT4.

The comparison of the order–disorder propensity in the proNT family seems to corroborate this divergent behavior of proNT4. Indeed, while the proNGF pro-domain shows a conserved propensity of structural disorder, excepting very few order propensities being observed in some mammalians (Supp. Fig. S4a), opposed to proNT4, proBDNF and proNT3, it shares a specific order–disorder propensity pattern through evolution in this domain. Indeed, both proBDNF and proNT3 pro-domains are characterized by a major transition to disorder from amphibia and reptiles to birds, non-placental mammals, mammals and primates (Supp. Fig. S4b and c). On the contrary, proNT4 shows an opposed behavior with a major transition to order increasing from amphibia and reptiles to mammals and primates (Supp. Fig. S4d).

This finding is further supported by comparing the distributions of the rates of evolution in each of the pro-domains (upper graphs, Fig. 2). It is worth noting that conserved shared features that are related to the biological function are considered more reliable for slowing amino acid substitution rates in disordered regions, which instead generally have a faster rate of evolution. To this regard, all the proNTs pro-domains except proNT4 are characterized by a disorder propensity. Moreover, besides proNT3 that shows the lower rate of amino acid substitution (Fig. 2c) along almost the whole pro-domain sequence, the regions of the lowest rate of evolution in the pro-domains of proNGF and proBDNF are coincident with amino acid sequences with secondary structure propensity (Fig. 2a and b). The fact that, most of the sequence of the pro-domain of proNT4 (but the amino and the carboxy terminals) is characterized by a low rate of evolution and by the absence of helical structure propensity prompts to the alpha helical region present in proNGF, proNT3 and proBDNF as a key structural determinant in mediating pro-domains biological functions.

hproNT4 shows a higher conformational stability compared to the other members of hproNTs family

The preliminary part of the biophysical characterization dealt with the expression and purification of the hproNTs. The obtained yields for hproBDNF were similar to those for hproNGF (10 mg/L culture), but despite any further optimization the yields for hproNT3 and hproNT4 were ten times lower (1–1.5 mg/L culture). Nevertheless, the purity of the resulting samples was > 95% (Supp. Fig. S5). The absence of aggregation was checked by dynamic light scattering experiments (data not shown), being the % of polydispersity 15, 13 and 16 for hproBDNF, hproNT3 and hproNT4, respectively. All these data confirmed the suitability of the samples for the subsequent biochemical and SAXS experiments.

In details, Differential Scanning Fluorimetry (DSF) spectroscopy was exploited in order to compare the thermal stability of the hproNTs. Indeed, monitoring protein thermal denaturation by using a high quantum yield hydrophobic fluorescent dye (SYPRO Orange) allowed to measure the melting curves, whose first derivative can be used to calculate Tm, a parameter that well describes thermal and conformational stability of proteins with low levels of surface hydrophobicity.

The designed experimental protocol resulted to be not suitable for the hproNT3 sample due to the strong interactions with the hydrophobic fluorophore occurring in the pre-transition region (data not shown). The resulting high initial signal did not allow monitoring protein melting by the fluorescence increase. This behavior is likely ascribable to the higher hydrophobicity profile of hproNT3 compared to the other human members of the family, as estimated by their respective GRAVY indexes.

Indeed, the thermo-melting profiles of the other hproNTs (Fig. 3a) show low and flat background fluorescence in the pre-transition region, and a bimodal profile that can be ascribed to the presence of the two domains. In details, the second sharp transition at high T is due to the unfolding of the globular folded domain of the mature proNT (confirmed by measuring the melting curve of the mature hNGF domain, data not shown). Instead, the first transition that occurs at low T, is likely the result of the denaturation of the pro-domains. In order to better compare pro-domains denaturation of the hproNTs, the respective melting curves have been normalized on the first transition (Fig. 3b).

Fig. 3

The pro-domain of hproNT4 shows an increased thermal stability compared to the ones of hproNGF and hproBDNF. Thermal denaturation assay by DSF with SYPRO© dye showing the overall melting curves in the sampled temperature range (a) and the melting profiles relative to the first transition (b), which describes the denaturation of the pro-domain. Experiments were performed in triplicate: the mean curve is shown for clarity.

Interestingly, the melting curves of the pro-domain of hproNGF and hproBDNF are almost superimposable (with Tm of 37.8 ± 0.5 °C and 38.3 ± 0.4 °C for hproNGF and hproBDNF, respectively), whereas the pro-domain of hproNT4 shows a significant positive shift in its melting curve (Tm 44.3 ± 0.3 °C). This reflects the higher conformational stability in solution of the pro-domain of hproNT4 with respect to the more flexible conformations of the pro-domains of hproNGF and hproBDNF. This experimental result corroborates the bioinformatic analysis of the structural order–disorder propensity in the pro-domain of proNTs, which highlighted a lower disorder content in the pro-domain of hproNT4 if compared to the other members of proNTs family.

hproNT4 shows a lower conformational plasticity in solution with respect to the other members of the hproNTs family

Solution techniques such as SAXS play a central role in the structural studies of IDPs, because they allow to gain insights not only on protein conformations in solution, but also on protein plasticity that is the main hallmark of the physiological functions of IDPs.

Thus, we performed batch SAXS measurements on serial dilutions of the four hproNTs samples (see Supp. Table S1, that summarizes the overall parameters from the SAXS data), without observing any concentration dependent association state changes. The resulting scattering curves of the four hproNTs are shown in Fig. 4 with the corresponding Guinier plots reported in the inserts.

Fig. 4

Experimental SAXS patterns and the obtained fits of the hproNTs: hproNGF (a), hproBDNF (b), hproNT3 (c) and hproNT4 (d). Experimental SAXS patterns, scattering calculated from the crystallographic models of the respective mature NTs by CRYSOL (fit crystal) and scattering calculated from models obtained by EOM (fit EOM) are plotted as indicated. The plots display the logarithm of the scattering intensity as a function of momentum transfer (s). The corresponding Guinier plots are shown in the inserts.

In all cases, the linearity of the Guinier plots at 0 < s < 1.3/Rg confirmed the homogeneity of the samples, although a small level of aggregation could be detected in both hproBDNF and hproNT4 samples. The MM values, derived from the calculation of the corresponding excluded Porod volume Vp, were in agreement with values expected for homodimeric species at all the measured concentrations, taking into account that for globular proteins the Vp in Å3

should numerically be about 1.5–1.7 times larger than the MM in Da, according to an empirical finding by Petoukhov et al. 2012 (Table S1). These findings were also confirmed by the estimation of the MM from the relative forward scattering intensity (s = 0, with s the scattering vector). The increased dimensions of the proNTs with respect to those of the corresponding mature NTs, have been highlighted by comparing the experimental scattering patterns to those computed using the high-resolution crystallographic structures and the models of the mature NTs (dot line in Fig. 4). This caused a poor fit to the experimental scattering curves (i.e., the discrepancy χ2 values are 5.233, 9.853, 4.241 and 2.754 for hproNGF, hproBDNF, hproNT3 and hproNT4, respectively).

The computed distance distribution functions p(r) (Fig. 5) of the hproNTs are in agreement with the typical features of IDPs, being very asymmetric if compared to the highly symmetric p(r) of globular proteins and displaying an extended tail, due to the variety of extended conformations that co-exist in solution.

Fig. 5

Comparison of the distance distribution functions of the hproNTs.

Interestingly, the structural parameters derived from the scattering curves depict a scenario regarding the proNTs conformations in solution that was hard to predict based just on the relative lengths of their primary structures. Indeed, the radii of gyration (Rg) and consistently the maximum dimensions (Dmax) derived from the p(r)s (Table S1) do not directly correlate with the MM of the hproNTs. In details, the hproNGF shows the most compact conformation with an Rg of 28.8 ± 0.01 Å and Dmax of 95 ± 3 Å. Consistently hproNT3, whose pro-domain is 17 residues longer, is slightly less compact with an Rg of 31.9 ± 0.05 Å and Dmax of 110 ± 5 Å. On the contrary, hproBDNF, whose pro-domain is intermediate in length among the former hproNTs, is by far the most elongated with an Rg of 37.2 ± 0.06 Å and Dmax of 125 ± 6 Å. The most striking result concerns hproNT4, whose extended conformation is comparable to that of hproBDNF, with an Rg of 35.0 ± 0.05 Å and Dmax of 115 ± 6 Å, despite its pro-domain is half in length.

Besides deciphering the molecular dimensions of a protein in solution, SAXS data are extremely useful when used in combination with other methods (in particular with high-resolution techniques such as X-ray crystallography), especially in the case of proteins containing IDRs. Indeed, SAXS measurements allow to infer the low-resolution shape of a protein in solution. Thus, to get insights of the mean conformation of the hproNTs, multiple runs of the program CORAL [79] were employed for the 3D reconstruction of the missing pro-peptide regions, by assuming either an asymmetrical conformation or by applying P2 symmetry (Supp. Figs. S6 and S7, respectively). The resulting models, which provide a representation of the average conformations in solution, corroborate the calculated structural parameters that point to larger dimensions of hproNT4 in comparison to hproNGF and hproNT3, despite the differences in the length of the respective pro-domains (for details, see Supplemetary Material).

As discussed, IDPs and multi-domain proteins that include large IDRs, like the proNTs family, are extremely dynamical structures in solution that are very hard to characterize at an atomic resolution by conventional methods especially in the case of proteins whose MM exceeds the limits of NMR. To this regard SAXS is not only suitable to determine averaged structural parameters and low resolution structures of flexible proteins in solution, but allows to analyze protein dynamic behavior providing the basis for the structural characterization of IDRs. Although SAXS typically yields low-resolution information, the combination of SAXS and atomistic simulations can provide insight into conformational changes in proteins and protein flexibility [45], [85], [90], [91]. In particular, the Ensemble Optimization Method (EOM) is a powerful tool that allows to facilitate the construction of conformational ensembles for which the ensemble-averaged theoretical SAXS profile is in agreement with experimentally determined SAXS profiles [85]. The resulting conformational ensemble provides a rich dataset that can be used to study the role of flexibility in protein function. Therefore, the development of this tool has deeply improved the characterization of flexible systems, and in particular of IDPs. Indeed, instead of only a mean Rg value, obtained by the Guinier plot and the p(r) profile just depicting an averaged behavior in solution, the EOM generates the statistical distributions of Rg and the representative conformational models that reflect the conformational heterogeneity of the macromolecule in solution.

Therefore, besides qualitatively confirming the intrinsic disorder porperties of human proNTs by Kratky plot analysis (data not shown), EOM was employed to obtain a quantitative characterization of plasticity and size distribution of possible multiple configurations of the four proNTs in solution. EOM searches for optimal conformational ensemble models of flexible proteins to match the experimental scattering profiles (i.e χ2 of 1.007, 0.595, 0.564 and 1.066 for hproNGF, hproBDNF, hproNT3 and hproNT4, respectively; see fit EOM in Fig. 4). As shown in Fig. 6, the behavior of the four proNTs in solution is significantly different.

Fig. 6

hproNT4 displays a reduced conformational heterogeneity in comparison with the other members of neurotrophin family. Rg distributions obtained by EOM for hproNGF (a), hproBDNF (b), hproNT3 (c) and hproNT4 (d). The distributions for the initial random pools of models and for the selected ensembles are shown by grey dot lines and black solid lines, respectively. The representative conformations (semitransparent surfaces) are shown near the distributions: compact on the left and extended conformation on the right. Folded homodimers of mature NTs are depicted as cartoons. Figure produced by Pymol [92].

Indeed, the profile of the Rg distribution of hproNGF closely resembles that of hproNT3 (Fig. 6a and c, respectively). Overall, both curves show a bimodal distribution with a sharp and symmetric peak in the range of 20–30 Å and an extended tail, consistent with their disordered nature that is more pronounced in the case of hproNT3, where the distribution of the selected ensemble is as wide as the distribution of randomly generated models. Although the peak accounts for the great majority of the hproNGF conformers population (75%), the corresponding peak of hproNT3 is less populated, being only 46.2% of the total. Accordingly, the tail of the hproNGF profile spans a shorter range (30–50 Å), which is less populated (25%) compared to the one of hproNT3 that covers a range of 30–62 Å (accounting for the 53.8% of the population with a significant residual of 5% at high range, i.e., 50–62 Å). Thus, these two hproNTs in solution are characterized by a high degree of conformational plasticity that is the typical hallmark of IDPs. This feature is even more dramatic in the case of hproBDNF (Fig. 6b). In fact, an almost continuous distribution of the selected ensembles is observed that covers a wider range with respect to the distribution of randomly generated models. The profile is still bimodal and it is characterized by two peaks. The main one (62% of the total population) is sharp and symmetric with Rg ranging 22–42 Å whereas the second peak is broader and asymmetric (38% of the total population), with a tail accounting for very extended conformations (0.5% of the total population), with Rg spanning between 42 and 70 Å and having a maximum around 49 Å.

Intriguingly, the profile obtained for hproNT4 is very divergent and points to a significantly lower conformational heterogeneity in solution. The Rg values of the individual members of the conformational ensemble (Fig. 6d) spanned a range shifted to higher dimensions (ranging from 28 to 46 Å), characterized by a roughly unimodal distribution with the main population (91.8%) spanning between 32 and 42 Å. Although the shift of the peak position towards higher Rg values accounts for the overall more extended conformation with respect to hproNGF and hproNT3 (in agreement with the parameters calculated for the mean conformation in solution), the distribution of hproNT4 is significantly narrower if compared to the other members of the hproNT family, confirming the predicted lower degree of disorder. Its reduced plasticity might be regarded as a key feature that prevents hproNT4 to fulfill proNTs biological functions, related to receptor binding.

Discussion

The present study aimed to understand the molecular determinants essential to mediate physiological functions of the intrinsically disordered proNTs. To this end, widely adopted biophysical and bioinformatics approaches have been exploited for their structural characterization. The results allowed us to shed lights on the dynamics and conformational landscapes of biologically active proNTs and identify specific and conserved structural motifs that are likely to be involved in the interactions with physiological binding partners.

Concerning the evolution of the proNTs, their molecular phylogeny suggests that this family is the result of two gene duplications [89]. The two different intermediate neurotrophin ancestral genes, resulting from the first duplication event, gave separately rise to two couple of paralogues (i.e., proNGF/proNT3 and proBDNF/proNT4, respectively). These duplications events took place early in vertebrate history, before the appearance of cartilaginous fish, but after the split of lampreys and hagfish from the common vertebrate lineage. The consequence of this early evolution is a high variability at the level of the pro-domains both in sequence length and in sequence composition, being the pro-domain of hproNT4 the most divergent. To this regard, the analysis of the molecular evolution of three avian NTs genes revealed that the high conservation exhibited by the BDNF pro-domain results from intense functional constraints that are relaxed in NGF and somewhat relaxed in the NT3 pro-domains respectively [93]. It is worth noting that the NT4 gene has apparently been lost during evolution in avian taxon, prompting BDNF to likely play a compensatory role in the development of the avian nervous system.

Previous structural studies on members of proNTs family revealed that the pro-domains of proNGF and proBDNF lack of folded structure [26], [29]. In agreement with these findings, the prediction of order–disorder propensity of the pro-domain of proNTs presented in this study, experimentally validated by DSF and SAXS measurements, also showed that the three members of this family but proNT4 are highly disordered in structure. To this regard, the present analysis on the evolution of the structural order–disorder propensity reveals an opposite transition of the order–disorder propensity for proNT4 with respect to other proNTs during the evolution of vertebrates that is most likely related to the loss of its biological functions (i.e., apoptosis and LTD). The EOM analysis of the SAXS patterns allowed to characterize the wide variety of conformational ensembles as result of the hproNTs high flexibility. hproNT4 shows a very low degree of disorder comparing to the others hproNTs. Indeed, it explores a rather restricted conformational space with a global low plasticity that in turn prevents from adapting to the physiological binding partners or to the changing environments of the cell, which on the contrary are known conditions to rapidly modulate IDPs flexibility and conformation [94]. Moreover, comparing the EOM profiles, it is tempting to speculate that the active conformers, which are able to fulfill the proNTs biological functions associated with receptors binding, are clustered in the main peak of the three bimodal distributions (ranging from 20 to 30 Å). These conformational ensembles can be considered in dynamic equilibrium with all the other population of conformers due to the high degree of disorder observed for hproNGF, hproBDNF and hproNT3. This crucial behavior might be considered as the structural determinant that provides flexible and rapid adaptation networks triggering proNT biological activities. On the contrary, hproNT4 shows a very limited plasticity in solution, sampling a very limited range of conformations that are more extended respect to the putatively active conformation ensemble of the other members of the family.

Even if IDPs are highly dynamic and do not form tertiary structures, they are known to contain variable amounts of transient secondary structures. The presence of these transiently formed structural elements has been reported to facilitate partner recognition and to tune the thermodynamics and kinetics of the binding interactions [95], [96], [97]. Interestingly, the present analysis on the evolution of secondary structure propensity highlighted conserved secondary structures in all vertebrates.

All proNTs share a conserved small beta strand secondary structure that might act as a pre-formed structural element involved in the chaperonin function of the pro-domains that are known to assist proper folding of the mature part.

Most importantly, the three biologically active proNTs are characterized by the predicted presence of an alpha helix secondary structure, which exhibits lower rates of evolution than average. This structural determinant is absent in the pro-domain of proNT4 and therefore might be envisaged as the structural basis mediating receptor recognition in a highly specific manner. Regarding this feature, previous studies both on hproBDNF by CD spectroscopy [98] and on hproNGF by hydrogen/deuterium (H/D) exchange and by combined NMR-SAXS measurements [46], [99] confirmed the presence of a localized higher-order structure motif in their respective pro-domains. The observation that the reported conserved alpha-helix structure is predicted to be proximal to the observed H/D exchange protection, can be ascribed to the fact that the first four backbone amide hydrogens at the N terminus of an alpha-helix are unprotected because they are not engaged in hydrogen bonding.

The findings that the disorder and the alpha-helix secondary structure propensities appear to be highly conserved throughout evolution of the three biologically active proNTs, prompt these features to likely play in concert to mediate receptor binding. Indeed, these structural determinants could act as functional modules within the frames of disordered regions in furtherance of the interaction with the globular domains of the binding partners, assisted by flexible structures around them [100], and thus be constrained to maintain the binding affinity to the specific receptors of the VPS10p family. Moreover, being predicted to be conserved in disordered structures across evolution, these elements that have lower rates of evolution are generally considered to be relevant to mediate physiological protein–protein interactions. Accordingly, the three proNTs, which play similar critical function in their interaction with specific receptors of the VPS10p family, are characterized by constrained disordered structures in vertebrates. These observations are in line with the well-known hallmarks of IDPs, by which essential intermolecular interactions and conserved binding partners could constrain disordered regions through evolution.

Conclusion

In conclusion, the disordered and dynamic nature of the pro-domain of biologically active proNTs have been addressed, for the first time, in the present study, by a combined evolutionary and biophysical approach. The present study represents so far, an important advance in knowledge with respect to the only previous reported phylogenetic analyses limited to the molecular evolution of three avian NTs genes, i.e., proNGF, proBDNF and proNT3 [93] and to BDNF in 36 vertebrates [101].

The gained deeper insights on the evolutionary conserved structural motifs appears to be the main determinant essential for the multifaceted roles of proNTs in physiological as well as in pathological contexts.

CRediT authorship contribution statement

S. Covaceuszach: Conceptualization, Formal analysis, Investigation, Visualization, Writing. L.Y. Peche: Formal analysis, Investigation. P.V. Konarev: Formal analysis, Investigation. D. Lamba: Conceptualization, Formal analysis, Investigation, Visualization.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.