text:
Nanomechanics combined with HDX reveals allosteric drug binding sites of CFTR NBD1

Graphical abstract

Introduction
Cystic fibrosis is a lethal disease caused by the functional defect of CFTR (Cystic Fibrosis Transmembrane Conductance Regulator) chloride channel in the apical membrane of epithelial cells [1,2]. Biochemical, in silico, and structural studies of CFTR have contributed to understanding the effect of mutations ranging from misfolding to impaired regulation and channel gating [2][3][4]. CFTR is a member of the ABC (ATP Binding Cassette) protein superfamily that provides an ion conductance pathway through the cell membrane via two transmembrane domains (TMDs) each consisting of six TM helices [5][6][7]. CFTR possesses two nucleotide-binding domains (NBDs), which bind ATP, form a ''dimer" via transient interactions, and regulate channel gating. An NBD is composed of a b-subdomain that binds ATP, and an a-subdomain that contains the ABC signature motif. The subdomains are not formed by sequential sequence regions but are intertwined (Fig. 1a). The binding and hydrolysis events are communicated towards the TMDs by the so-called coupling helices, which are the regions of intracellular ''loops" that interact with the NBDs [5]. Most of these structural features have recently been confirmed by cryo-electron microscopy [6][7][8][9]. Mutations cover every region of the protein, and many of them are located in the N-terminal nucleotide binding domain, NBD1 [4].
Primary structure alterations of CFTR usually impair channel function, folding, processing and/or stability [4]. The most frequent mutation, the deletion of F508 residue (DF508 or F508del) located in NBD1, impairs cooperative domain-domain folding by affecting the interactions between the four structured domain, NBD1/2, and TMD1/TMD2 by a poorly understood mechanism [10][11][12]. We have shown earlier that the hydrophobic side chain of F508 is crucial for the NBD1/CL4 interface coupling and may be a drug target for corrector development [5]. Furthermore, the NBD1/TMD assembly problem affecting global CFTR folding originates from local effects of DF508 in NBD1, including decreased thermodynamic and kinetic stability. Detailed studies of isolated wild type and DF508-NBD1 demonstrated a substantially decreased melting temperature of the mutant domain [13,14] while suggesting a similar folding pathway compared to the wild type [15]. It was also proposed that a partially unfolded state is responsible for the aggregation propensity of NBD1. As primary objective in drug development the restoration of the DF508-CFTR folding was advocated by either the correction of the NBD1 stability [16] or targeting both the NBD1 and NBD1-CL4 interface instability [17,18]. Both correction methods are likely effective only on an already folded NBD1 subpopulation. Notably, the only approved drug that can target the isolated NBD1 was suggested to stabilize an (un)folding intermedi-ate and not the native fold [19]. Small molecules are categorized as type I, II, and III correctors based on their primary targets of putative structural defects, represented by NBD1-TMDs, NBD2 and its interfaces, as well as the NBD1 stability, respectively [17][18][19].
Lumacaftor (VX-809) has emerged as a type I corrector of the DF508-CFTR folding defect [20,21]. Although VX-809 has a low efficacy in restoring the DF508-CFTR folding [17], serious efforts have been devoted to identify its binding site [22]. Experiments suggested that VX-809 exerts its action on the TMD1 and TMD1/ NBD1 interaction and also binds to the cleft formed by the Cterminal helices and b-strands S3, S8, and S9 [22][23][24]. The recently approved treatment for DF508, Trikafta, is a combination of Tezacaftor (VX-661, type I corrector), Elexacaftor (VX-445, type III corrector), and the gating potentiator Ivacaftor (VX-770) [21]. Tezacaftor restores the NBD1/TMDs interface and Elexacaftor was shown to act on NBD1 [19,25], but none of them known to directly Organization of CFTR NBD1 structure and secondary structural element detachment. Topology (a) and structure (b) of CFTR NBD1. b-strands are colored, and ahelices are marked in shades of grey. Black circles label the S6-a-S8 core region. Green and yellow helices: a-subdomain. (c) Distribution of the detachment time points of each secondary structural unit calculated from all MD trajectories with the wild type S6-a-S8 core (pulling velocity was 1 m/s). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) affect NBD1 folding. It would be important to improve the efficiency of this combinational therapy, since similar to long-term Ivacaftor treatment, Trikafta administration also exhibits partial restoration, thus may not be sufficient to prevent a long-term decline of lung function [21,26,27]. BIA (5-bromoindole-3-acetic acid) has been demonstrated to promote DF508-CFTR maturation and modestly stabilize the DF508-NBD1 against thermal unfolding as a type III corrector acting on NBD1 [18]. However, BIA was effective only in 1-10 mM range in in vitro experiments [16].
Since the lack of a folding corrector is likely the bottleneck in developing more effective combination therapies, improving existing or developing novel folding correctors would be critical for more efficacious CF therapies. Therefore, our major objective was to elucidate the mechanism of BIA correction. As this small molecule can potentially act on NBD1 folding, we aimed to gain insights on NBD1 folding and its folding intermediates at a higher resolution as compared to earlier studies [28,29]. To identify unfolding pathways and intermediates, we performed both force-probe molecular dynamics simulations and force spectroscopy experiments on the wild type and DF508-NBD1. We also studied the effect of BIA on NBD1 unfolding applying pulling experiments. Since BIA stabilized NBD1 against thermal unfolding [18], it likely remains bound after the completion of NBD1 folding. Thus, learning its binding site possesses therapeutic relevance, helping in the rational modification of BIA to enhance its potency or the design of novel binders. As a complementary approach, we used hydrogen-deuterium exchange (HDX) to uncover NBD1 regions with decreased conformational dynamics upon BIA binding. However, binding may decrease the dynamics in regions distant from the binding pocket via allostery. In order to rule out false positive segments stabilized allosterically, we employed in silico methods including binding site prediction and docking. Our results also contribute to the understanding of important folding and unfolding steps of CFTR NBD1 and the impact of F508 side chain and backbone deletion on these processes in the DF508-NBD1.

Structural models
Wild type NBD1 structure based on an X-ray structure (PDBID: 2BBO) from an earlier study [5] was used as the starting point. In order to match the construct used in our experiments, regulatory insertion (a.a. 405-435) was removed and the gap was sealed by loop modeling of Modeller [30], setting residues 403, 404, 433, and 434 as a loop region. DF508 mutation was modeled similarly. The missense mutations were generated using the mutagenesis tool of VMD [31].

Molecular dynamics simulations
Conventional all-atom MD simulations were performed with the WT and DF508 S6-a-S8 core region of the NBD1 structure. The S6-a-S8 core (a.a. 487-604) consists of three b-sheets (S8, S7, S6) and five a-helices (H3, H4, H4b, H5, H6), including the F508 residue. MD simulations were run using GROMACS 2019 [32] with the CHARMM36m [33] force field and the TIP3P water model. A 150 mM KCl concentration was used. Hydrogen atoms were replaced with virtual interaction sites to speed up the calculations, such that a 4 fs time step could be used [34]. Electrostatic interactions were calculated using the fast smooth PME algorithm [35]; the LINCS algorithm [36] was used to constrain bond lengths.
All structures were energy minimized using the steepest descent integrator in the first step, then equilibration (NVT, NPT) procedure was performed prior to each pulling simulation to generate inputs for independent simulations with different starting velocities (T = 310 K, p = 1 bar). The Nose-Hoover thermostat and the Parrinello-Rahman barostat with isotropic coupling were employed for the production runs. Time constants for the thermostat and the barostat were set to 2 picoseconds and 5 picoseconds, respectively. The C-terminus of the S6-a-S8 core structure was restrained, and the N-terminus was subjected to constant velocity pulling. Two pulling velocities, 1 m/s and 0.1 m/s were used, requiring 40 ns and 400 ns long trajectories, respectively. 50 simulations were performed with both constructs and for each of the two pulling speeds. MD parameter files can be downloaded from https://resources.hegelab.org.

Protein expression and purification
We used a cysteine-less [37], His6-tagged SUMO-fusion NBD1 carrying a deletion of regulatory insertion (DRI), which improves the protein stability and solubility [38]. It was especially important for the DF508 mutant [38,14,39]. Cysteine-less construct was used to avoid interfering with immobilization via terminally introduced cysteines. For simplicity, we referred to Cys-less NBD1 DRI construct as NBD1. SUMO-NBD1 constructs were purified from E. coli Rosetta 2 (DE3) pLysS strain. His6-tagged proteins were purified using an Ni-NTA affinity column (Profinity IMAC Ni-Charged, BioRad).

Single-molecule force spectroscopy experiments and analysis
Freshly cleaved mica surface was functionalized with APTES. The terminal cysteine of NBD1 was cross-linked to the APTEScoated mica using Sulfo-SMCC. Force spectroscopy was carried out on a Cypher atomic force microscope (AFM) instrument (Asylum Research) using PNP-TR cantilevers (spring constant: 100-200 pN/nm, NanoWorld). Experiments were performed at 25°C in PBS buffer (pH 7.2) [40]. Unfolding of NBD1-SUMO was carried out by first attaching the protein to the tip non-specifically by applying a constant force of 1 nN for 1 s to the tip on the mica surface, then followed by constant speed retraction with a pulling velocity of (1 lm/s).
Data was fitted using the worm-like chain (WLC) model of polymer elasticity using Igor Pro (Wavemetrics) extended with the Asylum Research AFM driving software. Since AFM tip and protein adhesion occurred at random locations, most of the retraction curves did not show the unfolding of the full protein. Forceextension curves exhibiting an overall length compatible with a completely unfolded NBD1-SUMO protein (total length of $110-120 nm) and including the SUMO unfolding fingerprint as a terminal contour length (Lc) increase of $25 nm were selected. All calculated Lc values were collected and summarized in histograms.

Analysis of unfolding pathways
The analysis was completed using GROMACS tools [32], the MDAnalysis package [41] and in-house Python scripts. For identification of secondary structural units' detachment over the pulling trajectories, the native contacts (contacts in the initial folded structure) of every secondary structural unit were determined and the fraction of native contacts (Q) was calculated as a function of time over each pulling trajectory using the soft_cut metric from the MDAnalysis package [41,42]. We defined a detachment event if Q decreased below a threshold value, which was set 0.2 in the case of Gō simulations with NBD1 and 0.1 in the case of conventional MD simulations with the S6-a-S8 core region. The structures of the unfolding simulations were clustered to identify intermediate states during unfolding. We used a pairwise contact-based RMSD with a cutoff value of 0.8 nm as a distance metric for clustering as described by Mercadante et al. [43,44]. The pairwise residue distance matrix of all unfolding conformations along every pulling trajectory was calculated and all values above the 0.8 nm cutoff value were set to 0.8 nm. The pairwise contact-based root-meansquare deviation (RMSD) values were calculated from these distance matrices, as.
where N is the number of residues and f i is the pairwise a.a. distance matrix of a given frame with position index i. Using this type of RMSD values for clustering led to ignoring the changes between distant residues that would have masked changes in important contacts within the remaining folded part. We applied a densitybased clustering algorithm, DBSCAN for clustering [45]. The structures of each pulling simulation were clustered separately using the contact based RMSD as a distance measure and the centroids of these clusters were pooled. The centroid structures from all simulations were clustered again to yield all observable intermediates from every simulation. Clusters with fewer members than 5% of the clustered structures were omitted from the re-clustering step. The DBSCAN clustering parameter, Eps, was set based on the Elbow method [45] and considering the Silhouette Score [46]. The cluster centroids were selected based on a calculated similarity score [47] as.
where D: pairwise distances of the structures of the cluster (contact based RMSD values), S ij : pairwise similarity of two structures, D ij (RMSD ij ): pairwise distance of structures i and j from the conformational ensemble (cluster), D std : standard deviation of D, argmax: the structure most similar to all other structures within the cluster (centroid).

Proximity calculations
The number of native and non-native contacts during each unfolding trajectory were calculated and normalized. Amino acid residues were in contact if the distance between any atom of the two residues is smaller than the 0.45 nm cutoff value. Contacts between residues that were present in 75% of the initial structure of the 50 simulations were labeled as native. Residues that were not native contacts but got closer than the cutoff in the course of the unfolding simulations were assigned as non-native contacts. Proximity values [48] were used to identify important non-native contacts where their cumulative number was increased during the simulations (18-25 ns and 180-250 ns range at pulling velocities of 1 m/s and 0.1 m/s, respectively). Proximity is minimal (zero) for the cutoff value or larger distances and maximal (one) if the distance between two residues is zero: where cut: distance cutoff between amino acids (0.45 nm), a i : amino acid residue in position i, d ij : distance between amino acids in position i and j. Proximity values for amino acid pairs in each structure were summed for all the WT and the DF508 simulations (separately for both pulling velocities) in the time interval of interest ( Supplementary Fig. 9). The summed 2D proximity values for each amino acid residue were also calculated by summing the previously calculated proximity values of the contacts of the given residue ( Fig. 4 and Supplementary Fig. 10).

Calculation of rupture forces and contour length increments in simulations
Rupture forces were derived from the pulling trajectories. Peaks of the force curves were collected for each SSE around the detachment time point of the given SSE. Normalized frequency of rupture forces from all and from individual SSE detachment events and their Gaussian density estimates were calculated and visualized. Contour lengths and their increments (DL) were calculated using the simple polynomial worm-like chain (WLC) interpolation formula [49].

Identification of drug binding sites in NBD1
FTMap webserver [50] (https://ftmap.bu.edu/) was used with default options to identify potential binding pockets on the NBD1 surface. Docking of BIA was performed with Autodock Vina [51]. Default options were applied except exhaustiveness, which was increased to 128. The search space was defined as described and shown in Supplementary Fig. 14.

Hydrogen-deuterium exchange (HDX) experiments
The DF508-NBD1 of human CFTR was purified as described [13]. Deuteriation time course of the DF508-NBD1 was measured by HDX coupled with mass spectrometry (HDX-MS) technique [53]. The sample concentration was 5 lM in buffer containing: 10 mM HEPES, 150 mM NaCl, 1 mM ATP, 2 mM MgCl 2 and 1 mM TCEP at pH 7.5. The deuterium uptake was performed in D 2 Obased buffer in the presence and absence of 2 mM BIA. For each deuteration time, NBD1 was mixed with 1:14 dilution ratio into D 2 O-based buffer, resulting more than 90% D 2 O contents, and incubated for 10 s, 40 s, and 120 s. HDX reaction was quenched by adding chilled quenching buffer (300 mM glycine and 8 M urea at pH 2.4) with 1:2 ratio. Quenched solution was flash frozen in MeOH containing dry ice and stored at À80°C until use. 10 lL of quenched sample was injected into the sample loop, followed by in an on-line immobilized pepsin column prepared in house. Online pepsin digestion was carried out at a flow rate of 50 lL/min for 1.5 min, and resulting peptides were trapped on a C18 trapping column (Optimized technologies, Oregon City, OR). Following desalting for 1.5 min at a flow rate of 180 lL/min, the peptides were loaded onto a C8 analytical column (1 mm i.d. Â 50 mm length, Thermo Fisher Scientific) and separated with Agilent 1290 Infinity II UHPLC system. Separated peptides were detected by LTQ Orbitrap XL (Thermo Fisher Scientific) in positive-ion mode for m/z 200 -2000 using electrospray ionization. For peptide identification, tandem MS (MS/MS) analysis in data-dependent acquisition mode with collision-induced dissociation was performed in separate measurements. All MS/MS data were analyzed in Proteome Discoverer 1.4 (Thermo Fisher Scientific). The deuteration were determined from triplicate measurements and the collected data were analyzed using HDExaminer 2.3 (Sierra Analytics). The relative deuterium uptake (%D) for each peptide was calculated by comparing the centroids of the isotope envelopes of the deuterated samples against the undeuterated controls. Deuterium uptake plots were generated using Prism 6 (Graphpad).

Differences in unfolding MD simulation with DF508 and WT S6-a-S8 cores
The structure of the CFTR NBD1 domain was investigated using a combination of experimental and simulation methods. To help the interpretation of the experimental AFM results, we first obtained a set of possible unfolding pathways of the whole wild type CFTR NBD1 that can be correlated with the experiments. Since pulling a large system such as the 250 amino acid (a.a.) long CFTR NBD1 in regular steered molecular dynamics simulations (MD) is highly limited, we performed simulations using an all-atom Gō model (SMOG2, Supplementary Text) [54]. To describe a pathway, we identified the unfolding steps of NBD1 by analyzing the hierarchy of the detachment of secondary structural elements (SSE) in each pulling simulation. This was determined by calculating the fraction of native contacts (Q) as a function of pulling time (Supplementary Fig. S1). We calculated the unfolding pathways of the NBD1 in the case of all simulations (n = 100) to attain an unfolding pathway set ( Supplementary Fig. 2). Analyzing these data obtained from the Gō simulations made the assignment of the unfolding steps of structural units to the AFM force-extension data possible.
Investigation of the unfolding pathways of Gō simulations revealed that the unfolding of NBD1 consisted of two parts. The unfolding always began with the unfolding of the b-sheet subdomain followed by the a-helical subdomain. The last steps constitute of the unfolding of the center of the a-helical subdomain, the S6-a-S8 core (a.a. 487-603, Fig. 1). Since several studies [15,29,55] have indicated that the S6-a-S8 core is crucial for NBD1 folding, we focused on this region and studied its unfolding in fully solvated atomistic force field pulling MD simulations for high resolution and accuracy. In the MD simulations of the WT and DF508 S6-a-S8 core we analyzed the order and the timing of SSE detachments and the rupture forces at which these detachments occur. The order of SSE detachments was determined by monitoring the fraction of the native contacts as was done for the Gō simulations ( Supplementary  Fig. 1). Detachment sequences of the secondary structure units were determined in each simulation, and the frequencies of the resolved pathways were calculated (Fig. 2a-c). In order to simplify the comparison between the WT and mutant pathways, we separated the S6-a-S8 core unfolding into two stages. In the first stage (Fig. 2b), the b-strands S6, S7 and S8 and the helix H6 unfold. Most frequently S8 decoupled in the first step, followed by H6 unfolding and the concurrent detachment of S6 and S7. The frequency of this pathway is much lower in DF508 S6-a-S8 core than in WT (54%, 27 out of 50 versus 78%, 39 out of 50; p < 0.05, v 2 test). The main divergence between the WT and DF508 protein was the different timing of the H6 detachment. An increased frequency of the pathways, in which the detachment of H6 occurs after the decoupling of S6 and S7, can be observed in DF508 S6-a-S8 core (DF508: 40%, WT: 22%; p = 0.0517, v 2 test). This could be caused either by weaker binding of S6 or by stronger binding of H6 to the folded part of the core when compared to the WT. Comparing the S6 rupture forces of WT and mutant shows that in the case of DF508 a greater proportion of S6 detached at lower forces, while the distribution of H6 rupture forces were unchanged (Fig. 2d), indicating that the mutation weakened the S6 interactions. At the second stage of the S6-a-S8 core unfolding, the a-helices H5, H4, H4b and H3 detached (Fig. 2c). H5 and H3 in the WT core showed a synchronized unfolding slightly more frequently than in the DF508 core (20%, 10 out of 50 versus 6%, 3 out of 50; p = 0.074, v 2 test). Comparing the time points of the SSE detachments of the WT and the mutant core, we observed faster detachment of secondary structural units (S6, H4 and H5) in the mutant compared to WT (Supplementary Figs. 3 and 4).
To assess the effect of speed on the unfolding forces, we performed pulling simulations also at a lower speed, which resulted in smaller forces and loading rates ( Supplementary Figs. 5 and 6). Nevertheless, the main pathways remained the same with somewhat shifted ratios; and some SSE detached earlier in the mutant than in the wild type ( Supplementary Figs. 4 and 7).

An increased number of non-native contacts are characteristic for late WT NBD1 intermediates
To characterize the details of different aspects of the unfolding steps, intermediate structures of the S6-a-S8 core from our MD simulations were determined and analyzed by using a special contact-based metric [43], which is efficient to compare highly different structures developing along pulling trajectories. Accordingly, the conformations from each pulling simulation were clustered separately based on contact RMSD. Then, the centroids of these clusters were pooled as intermediates, which were clustered again to yield all observable intermediates from every simulation. We identified four intermediate clusters describing the structural changes during unfolding in the case of both constructs (Fig. 3). However, there is an intermediate structure forming a well-defined wild type cluster (WT cluster #3), which was not observed for the DF508 S6-a-S8 core. Instead, the mutant conformations in the corresponding period (18-25 ns) changed their structure continuously into the totally unfolded state and were clustered into the last unfolding group.
To analyze the intramolecular interactions contributing to the formation of the well-defined intermediate structure observed in the WT cluster #3, the native and non-native contacts were calculated during unfolding. As expected, the number of native contacts decreased monotonically during unfolding (Supplementary Fig. 8c) and appeared very similar for both constructs. In contrast, the number of non-native contacts increased at the beginning of pulling simulations ( Fig. 4a and Supplementary Fig. 8e), probably reflecting the equilibration of the structure under pulling conditions, then exhibited a decreasing trend. Interestingly, a secondary increase was observed in conformations between 18 and 25 ns and it was more pronounced in the WT than in the mutant, coinciding with an unfolding intermediate state detected only in the WT (WT cluster #3).
In the next step, we analyzed which residues were involved in the non-native contacts observed in these conformations between 18 and 25 ns (Fig. 4, Supplementary Figs. 9-11). Per residue proximity values were exploited to determine those residues which participated in non-native contact formation (Fig. 4b, and Supplementary Figs. 10 and 11). Residues exhibiting pronounced differences between the WT and mutant proteins were located in the H3-H4 loop around the F508 residue and in the H4-H5 loop. Most of the residues in non-native contacts were hydrophobic and located in the gap between helices H3, H4 and H5 (Fig. 4c). The above studied non-native intermediate was the last one during unfolding, thus it is likely formed as the first intermediate in the reverse, folding process. This suggests that the F508 residue may play an important role in the interaction network of a non-native intermediate during the early folding that was also observed in our folding simulations ( Supplementary Fig. 12). Clustering and contact analysis were also performed on simulations using 0.1 m/s pulling speed (Supplementary Results and Supplementary Figs. 8-10) and indicated some differences from simulations at higher pulling speed.

AFM experiments revealed distinct unfolding steps of NBD1
For the pulling experiments, the NBD1 N-terminus was tagged with SUMO (small ubiquitin-like modifier) protein (Fig. 5a). We applied SUMO as a fusion protein, since it has been reported to enhance protein expression and it possesses a well-characterized unfolding pattern in force-extension curves obtained by AFM, aiding selection of successful experiments and their analysis [39]. The C-terminus of NBD1 was cross-linked to the mica-surface and the N-terminal SUMO tag was grabbed and pulled via non-specific binding. A successful pulling event is characterized by typical saw tooth profiles with peaks indicating rupture events, by a total length of $110-130 nm (length of the SUMO-NBD1 protein) and a terminal contour length increment of $25 nm (the SUMO fingerprint) (Fig. 5b).
In order to quantitatively compare the pulling events of WT and DF508-NBD1 constructs, contour length increments (DL) derived from the NBD1 unfolding, detected before the unfolding of SUMO, were collected from force-extension curves. For WT, the DL histogram exhibited five peaks (Fig. 5d). In order to correlate these peaks to molecular events, we used the unfolding pathway set of natively folded NBD1 obtained by Gō simulations and NBD1 structural information. In simulations, the last event was the breaking of the S6-a-S8 core. This core consists of two SSE groups, gS8 and gS6. The SSEs were grouped when they unfolded together in Gō simulations ( Supplementary Fig. 1). For example, the b-strand S8 and the a-helix H6 together form the group gS8. Since groups gS8 and gS6 are regions of 32 and 83 amino acids in length, respectively, their detachments result in contour length increments of $12 nm (32Â0.35 nm) and $30 nm (83Â0.35 nm). Depending on the order in which the gS8 and gS6 regions of a given NBD1 unfold, and whether they unfold separately or together, contour length increments of 12, 30 and 42 (12+30) nm resulted during the unfolding.
These data corresponded to well-defined peaks observed in the last part of the force-extension curve (Fig. 5). Breaking of group gS3 resulted in the unfolding of a $13 nm (36 a.a.) long segment. If this unfolding happened together with the breaking of gS8 and gS6, then the produced contour length increment was 55 (42+13) nm. By similar reasoning, breaking of gS10 and gS9 at the start of unfolding produced a 74 or 72 nm peaks in the force-extension curves, respectively. The first two peaks in the histogram of the contour length increments (Fig. 5d, upper histogram) were around 55 and 70 nm. Albeit these data were noisy, since they included tip positions close to the mica surface with adverse electrostatic inter- actions, the changes in contour length corresponded quite well to these specific breaking events of WT NBD1.
3.4. DF508 mutation decreased the proportion of natively folded S6-a-S8 core of NBD1 DF508-NBD1 pulling experiments were analyzed as that of WT. The histogram built from DF508-NBD1 data did not contain welldefined peaks as pulling WT NBD1 in AFM experiments, but the distribution of contour length changes was homogeneous. Additionally, in a large number of pulling experiments with the DF508 mutant, unexpected contour length increments were observed as peaks at approximately 20 and 35 nm compared to WT (Fig. 5d). These peaks suggested that a significant number of DF508-NBD1 exhibited modified mechanical resistance or incorrect folding.
In order to resolve these two mechanisms, our MD simulations and experiments were correlated. We calculated contour length increments from force-extension curves of our fully solvated atomistic simulations with the S6-a-S8 core using the WLC model as in experiments. We found that the histogram peaks at 42 nm and 30 nm from MD simulation with WT NBD1 match the experimental data of WT NBD1 suggesting that we were able to detect the gS8 and gS6 unfolding by AFM (Fig. 5d). Since we pulled a natively folded S6a-S8 core region in our simulations, the similarity of the DL in silico and in vitro peaks confirmed that these were the peaks characterizing the mechanical resistance of the correctly folded structure.
Unfolding of a-helices H3, H4 and H5 as force peaks can be observed in the simulations, albeit they did not emerge in experiments likely because of the buffered unfolding of a-helices in in vitro experiments [56] (see below, Supplementary Text, and Fig. 5).
The same analysis performed for the DF508 simulations resulted in a histogram, which was highly similar to the one from the simulations with the WT S6-a-S8 core (Fig. 5d). This was not unexpected, since only the natively folded and not misfolded NBD1 structures were known and applied in simulations. Importantly, the WT-like peaks were observed in spite of decreased forces ( Supplementary Fig. 5) in simulations, suggesting that the difference between WT and DF508 experimental pulling curves was mainly not caused by a decreased mechanical resistance of the mutant, but by its decreased folding yield. Counting the num- ber of AFM unfolding events corresponding to the breaking of the natively folded S6-a-S8 core (any combination of groups gS8 with 12 nm and gS6 with 30 nm DL) showed that the unfolding signature of the S6-a-S8 core was native-like only in 28% of the DF508 experimental curves that is exactly a twofold decrease when compared to the 56% of WT curves with these peaks (Table S1).

The corrector molecule BIA acts on the a-helical subdomain
We investigated the effect of BIA on the unfolding of the mutant NBD1 in AFM experiments. We compared the contour length increments of DF508-NBD1 in the presence and absence of BIA and found the appearance of a pronounced peak around 24 nm in the presence of this compound (Fig. 6a and Supplementary Fig. 13). The typical WT peaks of 42 and 30 nm were not observed in the presence of BIA, suggesting that this molecule did not restore the native, WT-like conformation or mechanical properties of NBD1, but it bound to and stabilized an NBD1 region, resulting in a peak around 26 nm. This peak is close to the H5-H4-H3 peak observed in simulations, indicating that BIA binds to this part of the NBD1 asubdomain. This observation also suggests that BIA binding provides improved stabilization of this region than a-helices would alone in the native structure, since with our AFM setup the helix unfolding was not detected in the case of WT (Fig. 5d).
Two peaks around 70 nm and 55 nm were also detected in the presence of BIA, which likely correspond to the b-subdomain and were observed in WT but not in DF508. Thus, we aimed to narrow the potential BIA binding sites using in silico methods. FTMap 44 identified several putative drug binding sites in NBD1 (Fig. 6b), but sites at NBD1/CL4 and NBD1/NBD2 interfaces were excluded, since binding to those locations would interfere with CFTR assembly and maturation. Two out of three sites in the a-subdomain were reinforced by Autodock Vina [51] docking ( Fig. 6c and Supplementary Fig. 14). One includes ends of H4, H5, H6 and loops, and another involves H6, the loop between S2 and S3, and H7.
Additional evidence for the BIA ligand-binding site was sought for by determining the backbone amide hydrogen deuteration of the isolated DF508-NBD1 in the absence and presence of BIA using the hydrogen deuterium exchange and mass spectrometry (HDX-MS) technique. The time course of the domain deuteration was monitored by continuous labeling for 15, 40 and 120 sec at 37°C. The accelerated deuteration of the DF508-NBD1 relative to the WT-NBD1 was consistent with thermal destabilization of domain that was reported previously [13,14,28]. The DF508-NBD1 HDX was partially suppressed at four a-subdomain peptides, including of H4 (a.a. 512-525), H5 (a.a. 559-568) and H6 (a.a. 591-594), and at a group with several secondary structure elements (end of H6, S8, H7, S9 and S10; a.a. 595-624) (Fig. 6c). Importantly, BIA was able to attenuate the HDX kinetics of the mutant below the WT-NBD1 in a single peptide of the H6 (a.a. 591-594) at all incubation times. According to in in silico docking calculations, the last residue of the H6 (a.a. 591-594) peptide, L594 sidechain serves as a docking surface to BIA (Fig. 6d). These observations strongly suggest that BIA bound this region of H6 and its stabilizing effect is communicated allosterically to distant parts of NBD1, also explaining the native-like AFM unfolding signature of DF508-NBD1 b-subdomain in the presence of BIA.

Discussion
In this study, we explored the mechanical unfolding of CFTR NBD1 domain using both atomistic molecular dynamics (MD) simulations and single molecule force spectroscopy (AFM) experiments. We found that deletion of F508 has a significant effect on NBD1 unfolding. Our results suggest that these effects were likely exerted by hindering the formation of crucial non-native intermediate states in the late stages of unfolding, thus potentially in the early stage of folding. Additionally, by combining HDX and AFM experiments with simulations we identified drug binding sites. HDX and AFM results indicated that BIA delays NBD1 unfolding and potentiates the formation of NBD1 intermediates, respectively.
Although the first step of CFTR biogenesis affected by the DF508 is the folding of NBD1, there are only a handful of studies investigating this process. Qu et al. demonstrated that DF508 increased the probability of off-pathway intermediates and affected an early folding step before the formation of the ATP binding site [57,58]. They used isolated NBD1 and measured its folding yield by light scattering and intrinsic Trp fluorescence, showing that the F508 residue affects the rate of maturation and suggested that the side chain of F508 makes crucial contacts during the folding process. Interestingly, it has also been demonstrated by Bartoszewski et al. that deletion of the three nucleotides resulting in DF508 mutation causes alteration in mRNA structure, leading to a decreased rate of translation [59]. Bali et al. also described that the lower translation speed resulted in altered CFTR conformations in metabolic pulse chase experiments [60]. A set of other experiments targeted NBD1 folding on ribosomes using truncation constructs and FRET [15,29,61]. These studies revealed that the folding of the N-terminal part involving b-strands S1-S6 took place while the F508 residue was in the ribosome tunnel and the deletion affected a later stage of domain folding [15]. Kim et al. found that the ribosome delayed the a-subdomain folding which was essential in the correct insertion of S7/S8 strands into the b-sheet core [29]. These results seem to be contradictory to that of Qu et al. [57,58]. However, the N-terminal b-subdomain may fold rapidly and independently according to FRET experiments by Kim et al. [29], but the state of this subdomain is likely not sufficiently mature to provide an environment for forming a correct ATP binding site [57,58]. Most likely the different levels or timescales of experiments provided data on different aspects of the folding process. Although the above studies on NBD1 folding were fundamental and agreed in the involvement of altered intermediate states, their resolution was low.
Due to the success of using simplified models in pulling simulations [62,63], we first used a native structure based Gō model that overcomes computational limitations associated with pulling simulations with NBD1-sized proteins. This simplified model was essential for the analysis of our AFM experiments, since it is highly challenging to correlate the peaks of the force-extension curves with molecular events and to deduce the exact order of the unfolding of protein segments without molecular modeling. The unfolding patterns observed in our simplified simulations matched the experimental data obtained from AFM experiments, thus we could identify the distribution of the unfolding pathways of the secondary structural elements. Considering that the accuracy of the simulations using Gō model is limited by the lack of explicit water molecules and non-native contacts, we also performed regular, fully solvated atomistic force field pulling simulations with a smaller part of NBD1, the S6-a-S8 core. A detailed analysis of the unfolding was performed on the results of these atomistic simulations of the S6-a-S8 core. We detected altered pathway frequencies and faster detachment of certain secondary structure elements in the mutant core ( Fig. 2 and Supplementary Fig. 3), suggesting differences in the interaction network around the b-strand S6 and in the final unfolding intermediate unit containing ahelices H4, H5 and H3. Importantly, our results showed that the F508 residue remained a component of the folded part of the wild type NBD1 almost until the end of unfolding (Figs. 3 and 4). The WT core exhibited a higher number of non-native interactions at a late stage of unfolding compared to DF508, suggesting that non-native interactions contribute to the stability of the late unfolding inter-mediate detected in the wild type S6-a-S8 core. These interactions included positions with known CF-causing mutations (DI507, V520F, L558S and A559T) that have been shown to affect the acore compaction of the nascent NBD1 during a critical window of folding [61].
Taken together, the F508 residue supports the development and persistence of non-native interactions that may be an important factor for off-pathway avoidance and self-chaperoning. The non-native contacts, which have been described to influence the folding free-energy barrier by Clementi et al. [64] and can become the rate-limiting step of protein folding according to Shao et al. [65]. These non-native contacts likely serve as a deceleration mechanism to provide time for the NBD1 polypeptide to acquire the right intermediate state before engaging the next step of folding. This was also supported by the in vitro translation experiments of Kim et al., showing that faster codons inhibited folding [29]. The same residues, which were involved in the formation of non-native contacts during unfolding ( Fig. 4 and Supplementary  Fig. 10), were also in contact in our folding simulations (Supple- Fig. 6. BIA partially restores the WT unfolding pattern and binds to the a-subdomain. (a) The histogram shows the experimental contour length increments measured from the SUMO-tag, in the presence of BIA (n = 126). Increments corresponding to the expected length increases based on WT pulling are indicated with dashed lines. Kernel density estimates show the distribution of contour lengths increases from WT S6-a-S8 core pulling MD simulations (blue). (b) Potential drug binding pockets in NBD1 were identified by FTMap. Some small molecule fragments bind to CL4 binding pocket (pink spheres) or NBD1/NBD2 interface (orange spheres). Green: NBD1, yellow: CL4, gray: other parts of CFTR. Binding of fragments (blue spheres) to rational pockets does not interfere with CFTR assembly. (c) BIA (red sticks) was docked to NBD1 using Autodock Vina. The locations exhibiting higher or lower decrease in HDX are colored pink or magenta, respectively. Corresponding HDX rates for these peptides in WT-(blue) and in DF508-NBD1 in the absence (orange) or presence (green) of BIA are shown as insets. The significant BIA inhibition on the DF508-NBD1 a.a.591-594 peptide deuteration is depicted by the bar plot. (d) The docked pose with the best binding score indicates BIA interaction with L594 beside R516 and S519 of the DF508-NBD1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) mentary Fig. 12), confirming that the residue F508 and its surroundings may serve as a folding nucleus. Earlier, a decreased folding time was observed for the DF508-NBD1 in Gō folding simulations, suggesting that the self-chaperoning of NBD1 was diminished [66].
Although the quantitative comparison of the S6-a-S8 core unfolding simulations with the AFM experiments by analyzing the contour length increments showed correlation, it also draws attention to limitations and differences of these methods. First, the histogram of contour length increments from experimental DF508-NBD1 data did not contain well-defined peaks, whereas the data from the DF508 simulations showed a WT-like pattern (Fig. 5d). Considering that we started the pulling simulations from the natively folded S6-a-S8 core region, and the distribution of the contour length increments of the experiments and simulations matched those obtained from WT, we presumed that these were the peaks characterizing the mechanical resistance of the correctly folded structure. Thus, in our simulations we did not pull misfolded conformations, which were present in experiments. An additional marginal difference is the lack of peaks in experiments corresponding to the in silico unfolding of a-helices. The difference likely results from the higher pulling speed used in simulations. Helices have been described to act as force buffers under mechanical stress [56], thus the mechanical resistance of helices is negligible at lower pulling velocity and does not result in considerable peaks in experimental force-extension curves.
Importantly, by combining experiments and computer simulations, we identified drug binding sites that are located on the surface of NBD1 away from the protein axis and exposed to the solvent. Therefore, drug binding to these regions is unlikely to interfere with CFTR domain-domain assembly and maturation ( Fig. 6b and c). In addition, drugs targeting these regions may not only rescue the volatile folding and stability of the a-subdomain, but potentially allosterically stabilize the b-subdomain, as confirmed experimentally by BIA binding (Fig. 6a and d). We also demonstrated by computational methods that secondary site mutations either in the b-subdomain or in the a-subdomain restored the WT-like allosteric network in the absence of F508 ( Supplementary Fig. 16). Because of this allosteric subdomain coupling, we propose that a drug rationally designed to bind the asubdomain, not only corrects DF508 and other mutations in the a-subdomain, but also has the potential to rescue CF mutations localized in the b-subdomain.
In summary, we found that the deletion of F508 has a significant effect on the unfolding pathways of NBD1 and accelerated the detachments of certain secondary structure elements compared to the wild type. Our results suggest that these effects were likely exerted by hindering the formation of crucial non-native intermediate states in the late stages of unfolding, thus potentially in the early stage of folding (Figs. 3 and 4 and Supplementary  Fig. 8-11). Furthermore, our results confirm that F508 residue maintains a network of non-native contacts and suggest a role in slowing down the translation, thereby aiding self-chaperoning. The experimental results suggest that the S6-a-S8 core is folded incorrectly in a significant portion of the wild type NBD1 population, and that misfolding is greatly enhanced by the F508 deletion.
We conclude that the a-subdomain is prone to folding instability. To overcome the effect of destabilizing NBD1 mutations, the asubdomain could be a potential target for designing small molecules that allosterically stabilize the full domain.

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.


In vitro Inhibition of HIV-1 by Cyclotide-Enriched Extracts of Viola tricolor

Since viral infectious diseases continue to be a global health threat, new antiviral drugs are urgently needed. A unique class of therapeutic compounds are antimicrobial peptides (AMPs). They can be found in humans, bacteria and plants. Plants express a wide variety of such defense peptides as part of their innate immune system to protect from invading pathogens. Cyclotides are non-classical AMPs that share a similar structure. Their unique topology consists of a circular peptide backbone and disulfide bonds. In previous studies they have been attributed to a wide range of biological activities. To identify novel cyclotides with antiviral activity, we established a library of plant extracts largely consisting of cyclotide-rich species and screened them as inhibitors of HIV-1 infection. Subsequent extraction and fractionation revealed four cyclotide-containing subfractions from Viola tricolor with antiviral activity. These subfractions inhibited HIV-1 infection with IC50 values between 0.6 and 11.2 μg/ml, and selectivity indices of up to 8.1. The identification and characterization of antiviral cyclotides and the determination of the antiviral mechanisms may allow to develop novel agents to combat viral infections. Therefore, cyclotides represent a natural source of bioactive molecules with prospects for development as therapeutics.

INTRODUCTION
Viral infectious diseases continue to be a global threat to human health (WHO, 2020), as best illustrated by the ongoing SARS-CoV-2 pandemic. The development of preventive or therapeutic antiviral drugs allowed to effectively treat several viral diseases, but most approved antivirals target only a limited number of viral pathogens such as HIV-1, HCV, or herpes viruses. In comparison, the arsenal of drugs available to combat the more than 200 known human pathogenic viruses is very limited with no specific antiviral therapies against most of them (Nii-Trebi, 2017;WHO, 2022). Therefore, new drugs are needed, especially those that are less prone to resistance induction (Mahlapuu et al., 2016;Browne et al., 2020) and that have a broad activity which prepares for reemerging or emerging pandemic viruses (Kuroki et al., 2021). One unique class of therapeutic compounds are peptides which attained much attention over the course of the 20th century (Lau and Dunn, 2018). Peptides are recognized for their high selectivity and efficacy and are relatively safe and well tolerated (Mahlapuu et al., 2016;Wang et al., 2022). Peptide drugs now account for~7% of the total number of approved pharmaceuticals (Al Musaimi et al., 2021) and this number is expected to increase in the near future because of the tremendous progresses made in peptide synthesis and purification (Wang et al., 2022).
Antimicrobial peptides (AMPs) can also be found in humans, bacteria and plants (Huan et al., 2020). Plants have developed mechanisms to protect themselves against invading pathogens such as viruses, bacteria, fungi, nematodes and insects (Mott et al., 2014;Nishad et al., 2020). Since they only have the innate immune response, a wide variety of peptides is expressed (Mott et al., 2014;Mammari et al., 2021), many of which are rich in cysteines. Plant cysteine-rich peptides have been isolated from thousands of species (Liu et al., 2017); they are extremely diverse with hundreds of different peptides produced by a single plant species. These peptides are classified into different families mainly based on their cysteine spacing, disulfide bond network and three-dimensional fold (Silverstein et al., 2007;Tam et al., 2015;Hellinger and Gruber, 2019;Li et al., 2021). One of these families are cyclotides, which are circular, disulfide-rich peptides and typically consist of 28-37 amino acids (de Veer et al., 2019). They comprise a typical structural architecture, termed the cyclic cystine knot (CKK) motif. This is characterized by a head-to-tail cyclized peptide backbone and interlocking arrangement of three disulfide bonds (Figure 1). Due to this unique topology, they exhibit exceptional resistance to thermal, chemical and enzymatic degradation (Colgrave and Craik, 2004). Cyclotides are nonclassical AMPs since they do not exhibit the typical AMP structure comprising a cationic amphipathic α-helical or βsheet structure (Mahlapuu et al., 2016).
Cyclotides do not occur in all plants, but were isolated from flowering plants-including many medicinal plants with ethnopharmacological relevance-belonging to Violaceae, Rubiaceae, Cucurbitaceae, Fabaceae, and Solanaceae (Gruber et al., 2008;Poth et al., 2011Poth et al., , 2012Nguyen et al., 2013;Burman et al., 2015;de Veer et al., 2019). Typically, cyclotides do not appear as singletons, but a single plant produces dozens of cyclotides (Hellinger et al., 2015;de Veer et al., 2019). Most plants have their own individual suite of cyclotides and most cyclotides are unique to a particular plant meaning that it is uncommon for one cyclotide to occur in multiple plant species (Gruber et al., 2008). However, variations in cyclotide content between tissues as well as over different seasons or in different geographic locations have been reported (de Veer et al., 2019). For example, in Viola tricolor ( Figure 1C), a member of the Violaceae family, we previously isolated and characterized 164 different cyclotides (Hellinger et al., 2015). This plant is a common horticultural plant, and has been used in traditional medicine for heat-clearing, detoxification, and relieving coughs (Tang et al., 2010). Viola tricolor is also widely used in Russian official medicines, and a monograph on this plant is even included in the Russian Pharmacopoeia (Shikov et al., 2014(Shikov et al., , 2017(Shikov et al., , 2021. Since violaceous plants at a global scale may be the source to as many as 150,000 individual cyclotides, this commercially available medicinal herb may be a suitable starting point for bioactivity-guided screening studies (Hellinger et al., 2015). As a part of the plant host defense, cyclotides were shown to be insecticidal (Jennings et al., 2001), hemolytic , anti-bacterial (Strömstedt et al., 2017;Slazak et al., 2018), anti-fungal (Slazak et al., 2021), anti-cancer (Pinto et al., 2018), antifouling (Göransson et al., 2004), and antiviral (Daly et al., 2004;Ireland et al., 2008;Wang et al., 2008;Gerlach et al., 2013). It is hypothesized that the diversity of cyclotides in a single plant may reflect a strategy to overcome resistance development by pathogens or pests (de Veer et al., 2019).
Therefore, cyclotides represent a rich natural source of bioactive molecules with prospects for development as therapeutics against cancer and infectious diseases (Mammari et al., 2021). Making use of the natural evolution, a huge variety of already optimized peptides can be assessed (Lazzaro et al., 2020;Atanasov et al., 2021). Peptides have a good efficacy compared to small molecules (Wang et al., 2022) and a good safety profile because their degradation products are natural amino acids and they are generally less immunogenic than recombinant proteins or antibodies (Mahlapuu et al., 2016).
In contrary, first-generation antivirals have severe side effects due to pour specificity (Mammari et al., 2021). Cyclotides are expected to be broadly active due to their direct mode of action destabilizing lipid membranes (Henriques and Craik, 2012; which will most likely prevent resistance development (Mahlapuu et al., 2016;Browne et al., 2020;Huan et al., 2020). Finally, peptides can be chemically modified and structure-activity-relationship studies can help to improve therapeutic features such as bioavailability, specificity and bioactivity (Browne et al., 2020;Huan et al., 2020;Dijksteel et al., 2021;Wang et al., 2022). Typical AMPs have a low oral and systemic bioavailability due to enzymatic degradation (Mahlapuu et al., 2016), but cyclotides have an exceptional resistance to degradation because of their unique topology (Colgrave and Craik, 2004), making them a promising class of antivirals (Lau and Dunn, 2018). And indeed, anti-HIV activity was shown for the prototypic cyclotide kalata B1, but also other cyclotides like kalata B8, varv E and circulins A and B (Gustafson et al., 1994;Daly et al., 2004Daly et al., , 2006Wang et al., 2008). Several classical AMPs are already on the market and further AMPs are currently evaluated in clinical trials (Mahlapuu et al., 2016). For example, the antiviral peptide T20 (enfuvirtide) is FDA-approved and commercially available (Huan et al., 2020). Other plant-derived peptides are also in clinical development. For example, griffithsin, a lectin isolated from the red alga Griffithsia sp., is investigated in clinical studies (NCT02875119; NCT04032717) as an anti-HIV microbicide for prevention of sexual transmission (Lee, 2019). Since it recognizes mannose, it has a broad-spectrum antiviral activity and was shown to inhibit SARS-CoV-2 recently (Cai et al., 2020). Additionally, the first cyclotide (called T20K) is currently in clinical trials for multiple sclerosis (Gründemann et al., 2019).
As nature-derived products are attractive candidates for translational application (Atanasov et al., 2021), we aimed to identify novel antiviral cyclotides. To this end, we established a library of plant extracts largely consisting of species with ethnopharmacological relevance that have previously been identified as a rich source of cyclotides or cyclotide-like knottin peptides. This encompasses ten plant species: Viola odorata, Viola tricolor, Psychotria solitudinum, Palicourea tomentosa, Carapichea ipecacuanha, Bryonia alba, Beta vulgaris, Sambucus nigra as well as Momordica charantia (He et al., 2013;Koehbach et al., 2013;Karpyuk et al., 2015;Fahradpour et al., 2017;Álvarez et al., 2018;de Veer et al., 2019;Retzl et al., 2020). Using solvent and solid-phase extraction yielded the plant extracts enriched in peptides. These extracts were analyzed for cytotoxic and anti-HIV-1 activity in cell culture-based infection assays. The resulting inhibitory concentration (IC 50 ) and cell toxic concentration (CC 50 ) were used to calculate selectivity indices and select the plants species for more detailed analysis. Using these criteria, the particular extract from Viola tricolor-a species used in traditional medicine and listed in the European Pharmacopoeia-had the most selective antiviral activity and was analyzed more thoroughly. To confirm the presence of cyclotides as bioactive molecule class in the extract, we applied (sub) fractionation by HPLC and mass spectrometry and identified several mixtures of cyclotides. These mixtures of cyclotides were subsequently analyzed again for IC 50 , CC 50 and selectivity indices. Lastly, in an approach to determine the antiviral mechanism, we performed pre-treatment experiments of the virions with the cyclotides. (leaves) and Palicourea tomentosa (Aubl.) Borhidi (formerly known as Psychotria poeppigiana) (leaves) were collected in La Gamba (Costa Rica) . Dry material was stored at room temperature in the dark and was powdered just before extraction.

Plant Extraction
The plant extracts and peptide-enriched fractions were prepared using previously established protocols (Muratspahić et al., 2021). Briefly, the peptides were extracted from dried plant material (50 g) using a mixture of methanol/dichloromethane, 1:1 (v/v), overnight and under continuous agitation at room temperature. Following removal of the plant material and filtration with 0.5 volume of ddH 2 O, the peptide-containing methanol/water phase was separated from the organic phase. The aqueous phase was subsequently concentrated using a rotary evaporator, lyophilized, and subjected to C 18 solid-phase extraction. The dried, crude extract was dissolved in 5% acetonitrile/95% ddH 2 O (v/v) and loaded onto the C 18 material ZEOprep 60 Å, irregular 40-64 µm (Zeochem, Uetikon, Switzerland). The column was equilibrated with solvent A (99.9% ddH 2 O/0.1% trifluoroacetic acid, v/v) and washed with 10-30% solvent B (90% acetonitrile/9.9% ddH 2 O/0.1% trifluoroacetic acid, v/v/v). The peptide-containing fractions were separated from hydrophilic components by elution with 50-80% solvent B, depending on the plant species. A mass fingerprint of peptideenriched fractions was determined by MALDI-TOF mass spectrometry (MS) using a MALDI-TOF/TOF 4800 analyzer (AB Sciex, Framingham, MA, United States) in a reflector positive ion mode with 2,000 to 10,000 total shots per spectrum and a laser intensity of 3,500. 0.5 µl of each sample was mixed with 3 µl of a matrix solution and spotted directly onto the MALDI target plate. An α-cyano-hydroxyl-cinnamic acid (α-CHCA) matrix (Sigma-Aldrich, St. Louis, MO, United States), dissolved in ddH 2 O/acetonitrile/ trifluoroacetic acid, 50/49.9/0.1% (v/v/v) with a final concentration of 5 mg/ml was used. Spectra were acquired, processed, and analyzed using the Data Explorer Software (AB Sciex). Most abundant peptides of extracts were labelled by molecular weight and, if applicable, identified by comparison to databases, such as Cybase (www.cybase.org. au; Wang et al., 2007) or previously published literature (Hellinger et al., 2015).

Bioactivity-Guided Fractionation of Viola tricolor
The bioactivity-guided fractionation of the cyclotide-enriched V. tricolor extract was carried as previously described (Hellinger et al., 2014;Muratspahić et al., 2021). The extract was dissolved in 5% solvent B and manually loaded onto the preparative Phenomenex Jupiter C 18 column (250 mm × 21.2 mm, 10 μm, 300 Å; Phenomenex, Aschaffenburg, Germany). Following preparative fractionation, a fractionation round 2 of V. tricolor was performed on a semipreparative Kromasil C 18 column (250 × 10 mm, 5 μm, 100 Å). The mobile phase consisted of solvent A (99.9% ddH 2 O/0.1% trifluoroacetic acid, v/v) and solvent B (90% acetonitrile/9.9% ddH 2 O/0.1% trifluoroacetic acid, v/v/ v). The preparative RP-HPLC fractions were collected automatically on a Dionex 3000 LC unit (Dionex, Amsterdam, Netherlands) machine using a linear gradient of solvent B between 5% and 65% at a flow rate of 8 ml/min while the semipreparative RP-HPLC were collected manually with a linear gradient of solvent B between 5% and 65% at a flow rate of 3 ml/min. Analytical RP-HPLC was performed on a Kromasil C 18 column (250 mm × 4.6 mm, 5 μm, 100 Å; dichrom GmbH, Marl, Germany) using a linear gradient of solvent B between 5% and 65% at a flow rate of 1 ml/min. The eluting peptides were monitored by UV absorbance at 214, 254, and 280 nm wavelengths. Cyclotides in V. tricolor were identified by molecular weight and retention time by comparison to CyBase (www.cybase.org.au; Wang et al., 2007) and data published in Hellinger et al. (Hellinger et al. , 2015).

HIV-1 Production
Virus stock of the CCR5-tropic HIV-1 NL4-3 92TH014-12 (Papkalla et al., 2002) was generated by transient transfection of 293T cells with proviral constructs as described (Münch et al., 2007). Transfection mixture was replaced by 2 ml DMEM supplemented with 2 mM L-glutamine, 100 units/ml penicillin, and 100 mg/ml streptomycin and 2.5% heatinactivated FCS after overnight incubation. 2 days later, virus was collected by centrifuging the cell supernatant for 3 min at 330 x g to remove cell debris. Virus stocks were stored at −80°C. The amount of virus to be used in infectivity assays was determined by titration.

HIV-1 Infection Assay
HIV-1 infection was quantified using TZM-bl reporter cells stably transfected with an LTR-lacZ cassette (Wei et al., 2002). Upon infection with HIV-1 the viral protein Tat is expressed, which activates the LTR promotor and results in βgalactosidase expression. For infection assays, 10,000 TZMbl cells were seeded the day before into 96-well plates. The next day, the medium was replaced with serumfree X-vivo 15 medium (Lonza, BE02-060F). For cell treatment, cells were treated with the titrated compounds and afterward infected. For virus treatment, virus and titrated compounds were mixed and preincubated and subsequently added onto the cells. 2 days later, infection rates were determined by detecting the β-galactosidase activity in cellular lysates using the galscreen β-galactosidase reporter gene assay system for mammalian cells (Thermo Fisher Scientific) and the Orion II microplate luminometer (Titertek Berthold). Measured values represent reporter gene activity (RLU/s) and were corrected for the background signal derived from uninfected cells. Untreated controls were set to 100% infection.

Cell Viability Assay
To assess toxicity of tested compounds, assays were performed in parallel to infection assays using medium instead of virus. 2 days later, metabolic activity was quantified via an MTT-based assay. To this end, medium was removed and 100 µl of 1:10 diluted MTT solution (Sigma Aldrich) was added. After 2.5 h at 37°C, supernatant was removed and 100 µl of DMSO-ethanol (1:1) solution was added to dissolve formazan crystals. Then, absorption was measured at 570 nm and baseline corrected for 650 nm using the Vmax Kinetic ELISA microplate reader (Molecular Devices, LLC). Untreated controls were set to 100%.

Bright-Field Microscopy
Cells were treated like for cell viability assessment. 2 days later, images of the cells were taking using a Cytation ™ 3 cell imaging system and processed with Gen5 (BioTek).

Nonlinear Regression
The determination of the inhibitory concentration 50 (IC 50 ) and cell toxic concentration 50 (CC 50

Screening of Cysteine-Rich Plant Extracts for Anti-HIV Activity
Driven by previous findings that cysteine-rich peptides exhibit anti-HIV properties (Bokesch et al., 2001;Ireland et al., 2008;Wang et al., 2008), we established a new library of plant extracts comprising cysteine-rich peptides, as previously described (Muratspahić et al., 2021). This library largely consisted of plants with ethnopharmacological relevance that have previously been identified as a rich source of cyclotides such as Viola odorata, Viola tricolor, Psychotria solitudinum, Palicourea tomentosa and Carapichea ipecacuanha Fahradpour et al., 2017;de Veer et al., 2019). Furthermore, the plant extract library included species that are rich in other cyclotide-like or knottin peptides such as Bryonia alba, Beta vulgaris, Sambucus nigra as well as Momordica charantia (He et al., 2013;Karpyuk et al., 2015;Álvarez et al., 2018;Retzl et al., 2020). Using solvent and solid phase extraction yielded the plant extracts enriched in cysteine-rich peptides. As quality control of the molecular content of peptideenriched extracts MALDI-TOF MS analysis was used, which confirmed the presence of peptide mass signals in the range of 2,500-4,000 Da (Supplementary Figure S1). Subsequently, these extracts were tested for antiviral FIGURE 2 | Effect of plant extracts on HIV-1 infection of TZM-bl cells. Ten peptide-enriched plant extracts were obtained using solvent and solid phase extraction, and their molecular content was analyzed by MALDI-TOF MS ( Figure S1). Extracts were dissolved in 5% DMSO. Serial dilutions of extracts #1 to #10 and DMSO control were added to TZM-bl cells, which were then infected with HIV-1 (infection, blue) or medium was added (metabolic activity, grey). After 2 days, infection and metabolic activity was quantified via β-galactosidase or MTT-based assay, respectively. Untreated controls were set to 100%. Shown are means ± SD of one (#10, metabolic activity), two (#10, infection) or three (all others) independent experiments performed in triplicates. *CC50 and IC50 were determined by nonlinear regression from concentration-response curves (Figure 2). Values represent means ± SD of three independent experiments in triplicates with the exception of extract #10 (only one (metabolic activity) or two (infection) independent experiments). #SI, was calculated by dividing CC50 by IC50.
Frontiers in Pharmacology | www.frontiersin.org May 2022 | Volume 13 | Article 888961 5 activity against HIV-1 by adding serial dilutions of the extracts to TZM-bl cells, which were then infected with HIV-1. All extracts reduced HIV-1 infection rates dosedependently, however they also reduced metabolic activity ( Figure 2). Hence, the 50% inhibitory concentration (IC 50 ) as well as the 50% cytotoxic concentration (CC 50 ) was calculated and used to determine the selectivity index (SI) which is the ratio of the toxic concentration of a sample against its effective bioactive concentration ( Table 1). This revealed IC 50 values between 0.4 ± 0.1 μg/ml and 278.7 ± 31.8 μg/ ml, but SIs for extracts #3 to #10 were low. In contrary, extracts #1 (Viola odorata) and #2 (Viola tricolor) had low IC 50 values of 3.3 ± 0.1 μg/ml and 1.8 ± 0.7 μg/ml, which yielded good SIs of 22.0 ± 3.9 and 31.8 ± 17.3, respectively.

Cyclotide-Enriched Fractions of V. tricolor Exhibit Anti-HIV Activity
Considering that extract #2 (V. tricolor) exhibited the most pronounced anti-HIV effect and had the best SI, we next conducted a multistep bioactivity-guided fractionation of this extract to identify the peptide(s) with anti-HIV properties. Accordingly, we performed preparative RP-HPLC thereby generating five distinct cyclotide-enriched fractions, referred to as fractions 1-5 (Supplementary Figure S2). These fractions were characterized by analytical RP-HPLC and MALDI-TOF MS to confirm the presence of cyclotides (Supplementary Table S1) and subsequently tested for anti-HIV-1 activity. Fraction 1 and 2 had no effect, whereas fraction 3-5 inhibited HIV-1 in concentration-dependent manner with IC 50 values of 111.1 ± 114.2 μg/ml, 21.0 ± 18.3 μg/ml and 20.5 ± 28.2 μg/ml, respectively ( Figure 3; Table 2). Since they also reduced the metabolic activity, the SI was calculated again, showing highest selectivity FIGURE 3 | Anti-HIV-1 activity of Viola tricolor peptide-enriched fractions. Extract #2 (V. tricolor) was fractionated and fractions dissolved in PBS. Then, fractions 1-5 and PBS control were titrated and added to TZM-bl cells. Cells were infected with HIV-1 (infection, blue) or medium was added (metabolic activity, grey). 2 days later, infection and metabolic activity was quantified via β-galactosidase or MTT-based assay, respectively. Data are normalized to untreated controls and represented as means ± SD of three independent experiments in triplicates.  Cyclotides in each subfraction were identified by molecular weight obtained by MALDI-TOF MS, as compared to CyBase entries for Viola tricolor (Wang et al., 2007). were titrated and added onto seeded TZM-bl cells in serumfree X-vivo medium, incubated for 30 min and subsequently infected with HIV-1 (cell treatment, blue). Alternatively, subfractions were first incubated with virus for 30 min. Then, these mixtures were used to infect the cells (virus treatment, red). 2 dpi, infection was quantified by β-galactosidase assay. In parallel, cell viability was quantified using MTT-based assay (metabolic activity, grey). (B) Experiments were performed like in A described. Additionally, 3 hpi medium was replaced by fresh compound-free medium. Data in A and B are normalized to untreated controls and represented as means ± SD of three independent experiments in triplicates. 0.7 ± 0.1 3.0 ± 0.5 12.9 ± 1.0 4.4 ± 0.4 0.6 ± 0.2 2.9 ± 0.9 20.6 ± 0.9 7.5 ± 1.8 mara.
3.8 ± 1.0 4.3 ± 1.8 --3.7 ± 1.4 3.3 ± 1.2 --CC50 and IC50 were determined by nonlinear regression from concentration-response curves ( Figure 4) and are given in µg/ml units for fractions and nM units for maraviroc. Values represent means ± SD of three independent experiments in triplicates. SI was calculated by dividing CC50 by IC50. Since 100 μg/ml was the highest concentrations used, no higher halfmaximal concentrations could be calculated. CT, cell treatment; VT, virus treatment; mara, maraviroc.
Since the cyclotide-enriched fractions 3-5 showed antiviral activity, all three will most likely contain antiviral cytlotides and could be used for future studies. For this study, we opted for fraction 4 and further generated four sub-fractions by semipreparative RP-HPLC fractionation, namely fractions 4.1-4.4 (Supplementary Figure S3). These fractions contained several cyclotides, which were identified and characterized by MALDI-TOF MS and RP-HPLC (Table 3, Supplementary Table  S2) including vigno 5, acyclic vitri E, acyclic cO22, cO28 (fraction 4.1), vitri peptide 2, vigno 9, cO2 (fraction 4.2), kalata S, varv C/D, acyclic vigno 5 (fraction 4.3) and kalata B1, varv E, vigno 3, vigno 4, chacur 1, and cO22 (fraction 4.4). The four fractions were then analyzed using the same experimental conditions as described above ("cell treatment"). In parallel, the virus was pre-incubated with the titrated subfractions and then used to infect the cells ("virus treatment") reaching the same final concentrations as under cell treatment conditions. In line with the observed activity of mother fraction 4, semipure V. tricolor cyclotide fractions 4.1-4.4 displayed anti-HIV activity and also reduced metabolic activity of the treated cells under both experimental conditions ( Figure 4A, Supplementary Figure S4). The IC 50 values were in the range between 0.7 ± 0.1 μg/ml and 11.2 ± 9.0 μg/ml. The highest activity was observed for virus pretreatment with fraction 4.4 with an IC 50 of 0.7 ± 0.1 μg/ml ( Table 4). The best selectivity had fraction 4.1 with a SI of 8.1 ± 5.4. All four subfractions where more active in virus than in cell treatment, suggesting that the inhibitory effect is directed against the viral particle. In contrary, the control inhibitor maraviroc, an entry inhibitor blocking the HIV receptor CCR5, had similar antiviral activities under both experimental conditions ( Figure 4A; Table 4).
As the antiviral activity is directed against the virus particles, we were wondering whether a reduced incubation time of the fractions with the cells may result in similar antiviral effects. To test this, we repeated the experiments using both, cell and virustreatment conditions, but removed the viral inoculum and the compounds 3 hours post infection. Cells were then further cultivated in pure medium, and cytotoxicity and HIV-1 infection rates were determined 2 days later as before. As shown in Figure 4B, the antiviral activities were in the same range (IC 50 s between 0.6 ± 0.2 μg/ml and 5.3 ± 2.3 μg/ml; Table 4) as compared to Figure 4A, further indicating that the inhibitory activity is rapidly directed against the virus. In contrary, cell toxicity was reduced ( Figure 4B, Supplementary Figure S5).

DISCUSSION
Plant defense peptides such as cyclotides are thought to be part of the innate plant host defense. The unique disulfide-knotted and head-to-tail cyclized topology may be one reason for the previously observed bioactivity of cyclotides (Daly et al., 2004;Ireland et al., 2008;Wang et al., 2008;Gerlach et al., 2019). To identify novel antiviral peptides, we established a library of plant extracts largely consisting of cyclotide-rich species with ethnopharmacological relevance and screened them against HIV-1, a representative and important member of the retrovirus family. Following the initial screen, the most active extract of Viola tricolor was selected to detailed analysis. Through two HPLC purification steps, four Viola tricolor-derived peptide subfractions were obtained which inhibited HIV-1 concentration-dependently with IC 50 values between 0.6 and 9.8 μg/ml. The antiviral fractions contained several cyclotides such as vigno 2/5, vitri 2, vigno 8/9, kalata S, vigno 3/4, kalata B1, and cycloviolacin O12/O22, belonging to the Möbius or the bracelet subfamily (Table 3, Supplementary Table S2). Cyclotides have been found in every member of the Violaceae family examined so far and they occur in all tissues of these plants, including roots, stems, leaves and flowers (Burman et al., 2015). These cyclotides have previously been identified in different Viola species (Ireland et al., 2006;Hashempour et al., 2013;Hellinger et al., 2015;Esmaeili et al., 2016), which is not surprising since more than 200 cyclotides are known within the Violaceae family (de Veer et al., 2019). However, only for some anti-HIV activity was shown to date (Gustafson et al., 1994;Daly et al., 2004Daly et al., , 2006Wang et al., 2008).
All subfractions tested were more active when the virus was pre-incubated, suggesting that antiviral activity is not mediated indirectly via a cell factor, but by the direct action of cyclotides on viruses. When the compounds were washed out 3 hours after infection, the antiviral activities were in the same range as when treated for 2 days, which is further evidence that the antiviral activity is mediated by the actual effect on the virions. This is in line with previously published data showing that cyclotides are able to permeabilize model membranes (Henriques and Craik, 2012;. It has been published that cyclotides target membranes through specific interactions with phospholipids containing phosphatidylethanolamine (PE) headgroups and then insert through nonspecific hydrophobic peptide-lipid interactions (Henriques et al., 2011Henriques and Craik, 2012;Cranfield et al., 2017). Subsequently, this promotes outward movement of PE phospholipids, exposing more PE in the outer leaflet. This self-promotes the binding of more cyclotides until a threshold concentration is achieved leading to self-aggregation, pore formation and eventual membrane disruption (Henriques and Craik, 2012). Whether the cyclotides identified in this work may exert similar effects needs to be further examined.
The analyzed cyclotide fractions also reduced metabolic activity, which is in line with previous studies Wang et al., 2008). This effect on the cells can probably be explained by the mechanism of action of the cyclotides acting on lipid membranes (Henriques and Craik, 2012). However, the activity of cyclotides is dependent on their affinity for lipid bilayers and the lipid composition. Although the total lipid in mammalian cells comprises~20% PE phospholipids (Patel and Witt, 2017), model membrane studies showed that kalata B1 targets HIV particles with high preference, since raft-like HIV membranes contain PE phospholipids and cholesterol/ sphingomyelin domains, which enhances the interaction (Brügger et al., 2006;Henriques et al., 2011;Mücksch et al., Frontiers in Pharmacology | www.frontiersin.org May 2022 | Volume 13 | Article 888961 8 2019). The preference for HIV particles over eukaryotic cells is also supported by the calculated selectivity index (SI) of up to 8.1 ± 5.4 herein. Such therapeutic windows have also been observed by others Wang et al., 2008), e.g. Wang et al. determined a CC 50 of 5.7 µM and an IC 50 of 0.66 µM resulting in an SI of 8.6. The ideal compound should have a relatively high toxic but a very low active concentration. The higher the SI ratio, the theoretically more effective and safer a drug would be during in vivo treatment (Cushnie et al., 2020). Besides, the hemolytic activity of several cyclotides isolated from V. tricolor has been determined previously with 50% hemolytic concentrations (HD 50 values) ranging from 4 to 226 µM suggesting a similar or even higher SI for red blood cell lysis (Tang et al., 2010).
In the future, it would be interesting to investigate if cyclotides have preferences for specific vesicles sizes. The so-called curvature-sensing peptides such as the AH-peptide preferentially target smaller vesicles and are antiviral (Jackman et al., 2018;Park et al., 2019). Membrane bending is enlarged in smaller vesicles, leads to lipid packaging defects between neighboring lipid molecules and subsequently to insertion of the amphipathic α-helical peptides (Kawano et al., 2019;Park et al., 2019). It cannot be expected that cyclotides have the same mechanism since they have a different tertiary structure, however, stronger membrane bending in smaller vesicles leading to packaging defects may also facilitate hydrophobic peptide-lipid interactions, pore formation and membrane disruption. This could lead to the preferential destruction of viral over cellular membranes improving the therapeutic window.
It is noteworthy that the antiviral fractions contained several cyclotides, which raises the question of synergistic effects. A wellknown advantage of plant substances is their complex composition, consisting of many compounds with multiple activities that together give a greater overall activity (Schmidt et al., 2008). Also, typical AMPs can have synergistic activity (Lazzaro et al., 2020). Therefore, mixtures of cyclotides could also have additive antiviral effects that reduce the tendency to develop resistance. Moreover, Gerlach et al. showed that the cyclotide cycloviolacin O2 increases the activity of the HIV protease inhibitors saquinavir and nelfinavir by pore-formation in HIVinfected cells and viral membranes at non-hemolytic concentrations (Gerlach et al., 2013(Gerlach et al., , 2019. It is also suggested to use AMPs in combination therapies to have higher activity but also to reduce resistance development (Dijksteel et al., 2021). Hence, cyclotides may augment antiretroviral therapy efficacy and be useful for combination therapies.
Since cyclotides act on lipid membranes it is reasonable to assume that they could be broadly effective antiviral agents. Indeed, the antiviral activity of a cyclotide against influenza A virus has been reported (Liu et al., 2014) and it would be interesting to investigate the antiviral activity of the extracts or purified cyclotides against other viral pathogens such as SARS-CoV-2. In addition, cyclotides have been successfully used for molecular grafting, i.e. inserting small bioactive epitopes into the stable cyclotide scaffold . Due to the plasticity of cyclotides, this usually preserves their structural integrity and activity  and can enable the development of dual-function antiviral agents by conjugating the cyclotides with specific virus entry inhibitors.
In sum, cyclotides and extracts described herein may offer a promising starting point for innovative therapeutic antiviral agents: the peptides have been structurally optimized by evolution to serve particular biological functions, and can be optimized for activity and stability for therapeutic applications (Atanasov et al., 2021). Moreover, cyclotides have an exceptional resistance to thermal, chemical and enzymatic degradation (Colgrave and Craik, 2004) which might increase bioavailability, and their mode of action has limited propensity for resistance development (Mahlapuu et al., 2016;Browne et al., 2020;Huan et al., 2020). Above all, V. tricolor is a medicinal plant, listed in the European Pharmacopoeia (Hellinger et al., 2015; EMA/HMPC/131734/2009), and its herein described properties may be explored in the future for applications of herbal preparations with antiviral activity.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CC performed HIV experiments, arranged figures, and drafted and revised the manuscript. EM and NT performed plant extraction and (sub) fractionation and supported the manuscript draft; JM and CG are responsible for the study, supervised the work and revised the manuscript; all authors critically reviewed the manuscript.

FUNDING
JM acknowledges funding through the DFG (CRC1279). Work in the laboratory of CG has been supported by a grant of the Austrian Science Foundation (FWF) through project P32190.