The transient receptor potential (TRP) channel superfamily plays a central role in transducing diverse sensory stimuli in eukaryotes. Although dissimilar in sequence and domain organization, all known TRP channels act as polymodal cellular sensors and form tetrameric assemblies similar to those of their distant relatives, the voltage-gated potassium (Kv) channels. Here, we investigated the related questions of whether the allosteric mechanism underlying polymodal gating is common to all TRP channels, and how this mechanism differs from that underpinning Kv channel voltage sensitivity. To provide insight into these questions, we performed comparative sequence analysis on large, comprehensive ensembles of TRP and Kv channel sequences, contextualizing the patterns of conservation and correlation observed in the TRP channel sequences in light of the well-studied Kv channels. We report sequence features that are specific to TRP channels and, based on insight from recent TRPV1 structures, we suggest a model of TRP channel gating that differs substantially from the one mediating voltage sensitivity in Kv channels. The common mechanism underlying polymodal gating involves the displacement of a defect in the H-bond network of S6 that changes the orientation of the pore-lining residues at the hydrophobic gate.
A fit cell must perceive and comprehend the conditions of its inner and outer worlds, integrating diverse and transitory physicochemical stimuli into concerted cellular decisions. For this reason, the membranes of even the simplest bacteria are studded with ion channel proteins that detect cellular conditions and translate them into electrochemical information via gated ionic conduction (Hille, 2001).
In eukaryotes, the complexity of cellular life has taken this requirement to its apex. Accordingly, natural selection has elaborated on the ion channel, producing an impressive array of polymodal cellular sensors, the transient receptor potential (TRP) channels (Clapham, 2003). All TRP channels detect multiple physicochemical stimuli, with some overlap among the eight extant TRP subfamilies. However, the response to each stimulus varies substantially from channel to channel, presumably dictated by heterogeneous and subfamily-specific intracellular and extracellular domains (Owsianik et al., 2006). Indeed, since their divergence from the voltage-gated potassium (Kv) channel superfamily over a billion years ago, TRP channel proteins have maintained a tetrameric six-transmembrane (6-TM) architecture and little else (Fig. 1, A–C).
Even so, there is reason to suspect that a common logic underlies TRP channel structure and function (Su et al., 2007). TRP subunits often form heteromeric channels in vivo, even if those subunits come from long-diverged TRP subfamilies (Kobori et al., 2009; Fischer et al., 2014). Also, functional cross-family TRP channel chimeras have been reported (Brauchi et al., 2006). For these heteromers and chimeras to assemble and function, one would expect a common mechanism of pore opening, even if the stimuli that trigger this mechanism differ substantially. These observations suggest that despite heterogeneity in domain organization and gating modality, selection has maintained a set of features in TRP channels underlying their aptitude for polymodal signal transduction.
Discovery of such features necessitates a comparative framework for TRP channels, but high resolution structural and functional data are scant for all but a few known members of the TRP channel superfamily (Liao et al., 2013; Winter et al., 2013). Protein sequences of thousands of TRP channels, on the other hand, are available and have been useful in subfamily-specific analyses of TRP channels (Matsuura et al., 2009; Doñate-Macián and Perálvarez-Marín, 2014). If a large ensemble of TRP channel sequences could be gathered and aligned, statistical features in the resulting multiple sequence alignment (MSA) such as conservation or covariation could be related to selective pressures common to all TRP channels (Halabi et al., 2009; Palovcak et al., 2014).
In this spirit, we searched for sequence–structure–function relationships common to the TRP channel superfamily. We scoured protein sequence databases for an exhaustive set of TRP channel sequences and developed a phylogenetically informed protocol to align the low similarity 6-TM region. Using approaches from information theory and machine learning on a large TRP channel MSA, we characterized the evolutionary pressure exerted on positions throughout the TRP TM region and inferred networks of contacting, coevolving residues running from helix to helix. These observations were contextualized with a parallel analysis on the distantly related Kv channels.
We report sequence features throughout the TRP channel TM core that are not shared with Kv channels, and offer molecular interpretations reasoned on the recently resolved TRPV1 structures (Cao et al., 2013; Liao et al., 2013). Our results establish an evolutionarily informed framework for investigating TRP channel structure and function and suggest a coherent, unified model of TRP channel gating.
MATERIALS AND METHODS
Collection and curation of exhaustive TRP and Kv channel MSAs
Annotated TRP channel sequences were collected from the expert-curated SwissProt sequence database and aligned with the fast statistical alignment (FSA) algorithm in the pore domain where most similarity among subfamilies is observed (Bradley et al., 2009). Included sequences spanned the TRPV, TRPM, TRPA, TRPC, TRPP, TRPML, TRPY, and TRPN subfamilies. A profile hidden Markov model (HMM) was constructed on the pore-only TRP MSA and searched against the uncurated UniProt sequence database with HMMER (Eddy, 1998). Thousands of putative TRP channel sequences were collected (E-value of <0.001). Voltage-sensor domain (VSD)-containing Nav, Cav, and BacNav channels were excluded by screening each sequence against an HMM previously constructed to identify VSDs (Palovcak et al., 2014). An approximately maximum-likelihood phylogenetic tree was inferred with FastTree2.1 using the JTT model of amino acid substitution (Price et al., 2010). Major branches of this tree were identified as known TRP channel subfamilies. Two pore and Catsper channels were also identified as major branches and excluded.
Each identified subfamily was realigned throughout the TM region (S1–S6) with FSA. HMM profiles of each subfamily were constructed and aligned with HHalign (Söding, 2005). Pairs of HMMs were aligned and progressively combined according to the following Newick-formatted phylogenetic tree: (((TRPV,(TRPA,TRPY)),((TRPC,TRPN),TRPM)),TRPP). The TRPML (mucolipins) could not be confidently aligned with our approach (the posterior probability of the profile–profile alignment was <0.50), likely because their sequences were highly similar and provided little information for the profile alignment algorithm. TRPML channels were subsequently excluded. The final TRP channel MSA contained 2,851 sequences. Histograms of pairwise sequence identity were calculated for each TRP subfamily by considering only the aligned 6-TM region, counting the matches for each pair and dividing by the length of the MSA. Sequence logos were generated with Weblogo (Crooks et al., 2004).
A Kv channel MSA was generated similarly. The initial seed MSA for the initial HMM was constructed from rat Kv1-9 sequences in the TM region. Sequences lacking 6-TM helices or the “T[VL]G[FY]GD” potassium selectivity filter sequence were manually removed. Small branches associated with calcium-activated potassium and Kv10-12 channels were also identified and excluded (Yu and Catterall, 2004). A tree inferred from the final Kv MSA contained branches representing Kv1-9 and two large clades of microbial Kv channels. The final Kv channel MSA contained 3,746 sequences.
Relative entropy calculations for the TRP and Kv superfamilies
To compare structurally equivalent positions, a structure-based alignment of apo TRPV1 to the paddle chimera (Protein Data Bank [PDB] accession nos. 3J5P and 2R9R, respectively) was generated with MultiSeq (Roberts et al., 2006; Long et al., 2007; Liao et al., 2013). We used the apo “closed” structure of TRPV1, as we found alternative conformations of TRPV1 produced identical structure-based sequence alignments. Because of differences in relative orientation between the two structures, the sensor and pore domains were aligned separately. A discrepancy in the structural alignment of the S6 helix was resolved by considering the pore domain of the bacterial voltage-gated sodium channel NavAb structure (PDB accession no. 4EKW), as discussed in the Results (Payandeh et al., 2011).
Amino acid frequency distributions were calculated as described in Palovcak et al. (2014), with sequences with >90% identity down weighted as a single sequence to avoid oversampling of highly similar sequences. Relative entropies for TRP and Kv MSAs were calculated according to:
where X represents the TRP or Kv MSA, the subscript i denotes the position in the MSA, and a is each amino acid in the 20–amino acid alphabet q. Qref is a previously calculated amino acid reference distribution for membrane proteins (Palovcak et al., 2014).
Average polarity calculations for the TRP- and Kv-sensor domains
For each position in the sensor domains of TRP and Kv channels, the average polarity was calculated as the expected value of the amino acid frequency distributions for that position (normalized to exclude gaps) and a biological amino acid hydrophobicity scale, as described in Chowdhury et al. (2014). In this scale, greater values indicate a larger free energy associated with inserting a monosubstituted TM helix into a membrane; more polar residues have greater values, given in kilocalorie per mole.
Inferring networks of evolutionary coupled residues and estimating their significance
The asymmetric pseudolikelihood-maximization direct coupling analysis algorithm (aplmDCA) was used to estimate the parameters of the maximum entropy probabilistic model consistent with selected MSA statistics (univariate and bivariate frequency distributions). Default parameters were used for field and coupling regularization and sequence reweighting (lamba_h = 0.01, lambda_j = 0.01, theta = 0.1). The optimal parameters are used to calculate an “evolutionary coupling score” (EC score) for each pair of positions given by the quadratic mean of all coupling parameters involving those positions (Ekeberg et al., 2013). Hereafter, we will call the matrix of all ECs between all pairs of positions the “EC map.”
In a fully connected model, the number of parameters increases quadratically with sequence length. Specifically, the number is 400((N(N − 1))/2)+20N, where N is the MSA length. To reduce the number of parameters to be estimated, we calculated the EC maps separately for the sensor domain (S1–S4), pore domain (S5–pore–S6), and the interface (S1, S5, S5–pore–S6) in both the TRP and Kv superfamily. Still, with so many parameters, overfitting is a serious concern. To assess the robustness of EC maps with respect to a different choice of the training dataset used in the statistical inference, we performed an error analysis to assign statistical significance to EC scores.
Our goal is to calculate the standard error of the maximum-likelihood estimator used in the statistical inference step. Ultimately, to gauge the “strength” of a coupling, we would like to take into account not only its value but also the statistical fluctuations resulting from the finiteness of the sample. Expressing each coupling in units of its standard error (in the form of a Z-score) allows us to assess its statistical significance and provides an expedient metric to appreciate, simultaneously, the strength and degree of robustness of each evolutionary-coupled pair. Moreover, the statistical significance enables us to associate a p-value to each threshold chosen to select the top ranking couplings. Calculation of the standard error of a maximum-likelihood estimator is, in general, a complex task. A conceptually simple scheme to achieve this goal, the so-called jackknife approach, is based on the use of appropriately chosen subsets of the original dataset (Reeds, 1978). Specifically, the jackknife estimator of a parameter θ is obtained by calculating the former for all subsamples resulting from the removal of a single observation without replacement, i.e., for each so-called pseudo-observation (Quenouille, 1949; Shao and Wu, 1989). Here, we have used a variation of this approach, the delete-d jackknife, which is obtained by omitting d observations each time, because it is statistically consistent for a wider class of estimators compared with the delete-one version (Efron and Tibshirani, 1993). In particular, the variance of the estimator is:
where N is the sample size, and d is the number of omitted observations. The sum runs over all the distinct subsets of size d that can be extracted without replacement from the original dataset. For each of them, a pseudo-observation is built, and the parameter is recalculated. Finally, is the average of all the terms; namely,
Note how the prefactor multiplying the sum in Eq. 3 compensates for the increased mean-square fluctuation as d increases. This estimator of the variance is consistent, i.e., it converges in probability to the true variance as the number points increase, as long as d is comprised between and N (Shao and Tu, 1995). An important byproduct of this analysis is a bias-corrected estimator of the parameter (i.e., with bias terms proportional to 1/N removed):
which enables comparisons between datasets of different size (here, is the parameter calculated on the entire dataset) (Narsky and Porter, 2014). A convenient choice for d is N/2, which makes the prefactor of Eq. 3 equal to 1. In practice, exhaustive enumeration of all the possible pseudo-observations is not computationally feasible; thus, we approximated the mean-square fluctuation term of Eq. 3 by stochastically sampling K pseudo-observations (Gentle, 2009).
The parameters, in the equations above, that we are interested in here are the EC scores provided by the aplmDCA estimator. Following the approach outlined above, we computed the EC scores 100 times by omitting half of the sequences from the MSA each time. Thus, the ECZ or Z-score of the EC is defined, for each (i,j) pair as:
where denotes the average over the K resampled MSAs (K = 100 in this case), and is the standardized EC score:
In Eq. 7, indicates the average over the distinct residue pairs. The standardization step is performed to homogenize all the pseudo-observations by removing constant terms and scale factors. In the following analysis, we will consider ECZ scores above 4 as corresponding to a one-tailed p-value of 3.2 × 10−5.
Online supplemental material
Fig. S1 shows a small MSA of representative TRP channel sequences aligned in the TM region, with the residues considered for the analysis highlighted by red shading. Fig. S2 shows the relative entropy of the TRP channel TM region with respect to Kv channels. Mapping of the relative entropy on the structure of TRPV1 highlights the regions where the distributions of residues are most similar: positions with low relative entropy primarily putatively face the lipids, whereas the core-facing regions of both channel families have distinct residue distributions except for the indicated lower portions of S3 and S5.
Profiles of TRP subfamilies enable alignment in the TM region
A reliable alignment of the TRP channel TM region is needed to reflect the evolutionary events of amino acid substitution, insertion, or deletion that have shaped TRP sequences. Additionally, insights obtained from such a dataset can only be generalized to the superfamily if the MSA is comprehensive. However, in the majority of the TM region, TRP channel sequences are highly heterogeneous, precluding sequence collection and curation by familiar methods. To overcome this obstacle, we developed a custom alignment protocol for the TRP superfamily.
Most simply, we leverage knowledge of TRP channel phylogeny to construct statistical profiles (profile hidden Markov models, or HMMs) of each TRP subfamily (Eddy, 1998; Price et al., 2010). These subfamily profiles are then aligned in pairs, reconstructing the full superfamily branch by branch (Söding, 2005). In this way, we were able to generate a plausible MSA of thousands of TRP channel sequences throughout their TM region. This method works because the statistical profiles of two protein families (for instance, the TRPV and TRPA subfamilies) contain much more information for a pairwise alignment than two representative sequences (such as human TRPV1 and TRPA1) (Marti-Renom et al., 2004). A detailed description of this workflow can be found in the top row of Fig. 2 A and in Materials and methods.
To compare the sequence features of the TRP channel superfamily with those of the well-studied Kv channels, we used the same protocol to construct a Kv channel MSA of similar size (described in Materials and methods). The final TRP MSA contained 2,851 sequences, whereas the Kv channel MSA contained 3,746 sequences.
The TRP channel MSA captures the expected diversity of the TRP superfamily
From the phylogenetic tree inferred from the TRP MSA, we identified eight subfamilies varying in size and diversity. However, the small and homogeneous TRPML (mucolipin) subfamily could not be confidently aligned with our protocol and was subsequently excluded from the analysis (see Materials and methods). For the remaining seven subfamilies, histograms of sequence identity showed major modes between 14 and 44% identity, with the major mode of the full MSA showing 16% identity (Fig. 2 B). The multimodal nature of the subfamily histograms mirrors the finer phylogenetic structure observed in the inferred tree (Fig. 2 C). Consistent with recent phylogenetic analyses on smaller trees of TRP channels, we find homologues of all TRP channel subfamilies (except TRPN) in single-celled eukaryotes and a highly diversified and expanded TRPA subfamily in insects (Cai, 2008; Peng et al., 2015).
Although the TRPY and TRPN subfamilies were relatively small (127 and 48 sequences, respectively), TRPCs, TRPVs, TRPAs, TRPPs, and TRPMs each contained hundreds of putative members, implying that the TRP channel MSA is heterogeneous and does not over-represent one or two TRP subfamilies. We therefore expect that insights generated from this dataset should be general to the TRP superfamily. This strategy has been validated by the structural determination of another member of the TRP family, namely TRPA1 (Paulsen et al., 2015). Structural superposition of TRPV1 and TRPA1 induces a sequence alignment that is completely consistent with the one inferred from the profile–profile matching approach. We are thus confident that the alignment is accurate across the entire TRP family.
Relative entropy can serve as a proxy for evolutionary pressure
Guided by a structural alignment of TRPV1 to the Shaker-type paddle chimera Kv structure (Long et al., 2007), we were able to directly compare the distributions of amino acids observed at equivalent positions in the TM region of each superfamily (the colored regions in Fig. 1, B and C). In large, diverse protein families, the nature of the amino acid distribution observed at a certain position is indicative of a common evolutionary pressure acting throughout the family. Visual inspection of aligned sequence logos of the sensor domain show that despite substantial structural similarity, very few positions share similar amino acid distributions between TRP and Kv channels (Fig. 3 A). To make this observation more quantitative, we calculated the relative entropy for each equivalent position in the MSAs.
The relative entropy (or Kullback–Leibler divergence) is a nonsymmetric measure of the difference between two probability distributions (Thomas, 2006). For each structurally equivalent position between two MSAs, the relative entropy quantitatively describes how different the distributions of amino acids at those positions are (see Materials and methods). The relative entropy between TRP and Kv channels is lowest for lipid-facing residues and much higher within the core of the sensor domain, throughout the ion conduction pathway in the pore domain, and at the intersubunit interface (Fig. S2 A). The S3 and S5 helices are exceptions, showing low relative entropy between TRP and Kv channels (Fig. S2 B). In general, two limiting scenarios can give rise to the observed low mutual entropy: both the distributions are similar to the uniform distribution (high entropy), or they deviate from it in a similar fashion, i.e., they are both low entropy with coincident modes. Specifically, although we expect the distribution of the lipid-facing residues to be random-like, the low mutual entropy of S3 and S5 is potentially the signature of shared (although loose) patterns of conservation. To distinguish between these two alternative hypotheses, we calculated the mutual entropy with respect to a common reference distribution Qref, for the two MSAs separately, with the ultimate goal of detecting conserved patterns of residues. In fact, the low mutual entropy of S3 and S5 turns out to reflect conserved patterns. The functional implications of this conserved sequence feature are discussed in details in the following sections concerning the sensor and pore domains. Technical details of this analysis are described in Materials and methods section.
Few common sequence motifs are observed in the TRP- and Kv-sensor domains
It is unclear what role the sensor domain plays in the structure and function of TRP channels (Kalia and Swartz, 2013). Certain positions at the base of S3 show similar distributions of amino acids in both TRP and Kv channels (TRP positions 507, 513, and 516; Kv positions 307, 313, and 316). This motif is conserved in most sensor domains (Kumánovics et al., 2002; Palovcak et al., 2014) and has recently been implicated in the gating of TRPM8 and TRPM2 channels (Winking et al., 2012).
Other positions have high relative entropy in both superfamilies but different amino acid compositions, implying that the position is important in both families but for different reasons. For example, TRP S1 positions 441 and 444 are conserved aromatic residues, whereas the equivalent Kv positions 237 and 240 are conserved hydrophobic and polar residues, respectively. In the voltage-sensitive Kv channels, positions 237 and 240 occupy the gating pore through which the positively charged S4 gating charges are shuttled (Long et al., 2007). Their mutation strongly affects the voltage sensitivity and kinetics of gating charge transport (Lacroix et al., 2014). To our knowledge, positions 441 and 444 have not been mutated in TRPV1, but in TRPM8, mutations of the equivalent Y745 (erroneously described as an S2 residue in Bandell et al., 2006) specifically decrease menthol sensitivity (Winter et al., 2013).
To visualize the spatial distribution of evolutionary pressure through the sensor domains, we map relative entropy onto the TRP and Kv structures (Fig. 3, B and C). Both feature high relative entropy positions at the intracellular interior of the sensor, whereas in Kv channels, high relative entropy positions are also found at the extracellular interior (Kv 247, 248, 283, 286, 322, 324, 362, and 365) and correspond mostly to residues such as the acidic countercharges and the basic gating charges. As expected, no putatively lipid-facing positions in either family had relative entropy greater than our chosen threshold (1.25 nats).
Solvation of the TRP sensor is limited to the intracellular side
In the Kv channel–sensor domain, polar residues inside the sensor act to solvate the gating pore, with upper and lower crevices forming an hourglass-shaped vestibule of water, which focuses the membrane electric field onto the gating charges, allowing efficient transport of charge (Starace and Bezanilla, 2004; Krepkiy et al., 2009). In the TRPV1 structure, the lower vestibule is splayed open and occupied by the polar face of the C-terminal TRP helix, suggesting that it is solvated (Liao et al., 2013). The upper portion of the vestibule, on the other hand, is densely packed with bulky aromatic and hydrophobic residues and is therefore likely unsolvated.
It is not clear whether this asymmetric solvation is unique to TRPV1. To test this, we measured the average polarity in the TRP- and Kv-sensor domains (following a protocol described in Chowdhury et al., 2014; see Materials and methods) and mapped positions with conserved high polarity (on average, more polar than serine) onto the TRP- and Kv-sensor domain structures (Fig. 3, D and E). Although the Kv-sensor domain shows the expected hourglass shape, the TRP sensor recapitulates the pattern of polar residues observed in the TRPV1 structure, suggesting that an asymmetrically solvated sensor domain is a conserved feature of TRP channels.
Networks of evolutionarily coupled residues can be revealed through statistical inference
Although the relative entropy view of evolutionary pressure is intuitive, it misses a great deal of the information that can be extracted from such large MSAs. Evolution may act on more variable positions by strictly tuning their covariation with other residues (Süel et al., 2003; Halabi et al., 2009). This is especially true for residues that make physical contact in the native protein structure (Göbel et al., 1994; Morcos et al., 2011). In general, directly measuring this covariance with statistics such as the mutual information does not produce easily interpretable results, as transitive interactions between chains of residues can yield indirect correlations (Weigt et al., 2009). In Weigt et al. (2009), the authors showed that a probabilistic model assuming a direct interaction between coevolving residue pairs can explain the observed frequency distribution of amino acids in a large protein family MSA. Strong interactions, or “evolutionary couplings,” are highly predictive of native protein contacts (Ekeberg et al., 2013).
Only a subset of contacting residue pairs in a protein are detected as ECs. We recently examined the network of ECs present in VSDs and found that the structure of this network was highly informative about the general mechanism by which these molecular machines sense voltage (Palovcak et al., 2014).
We calculated EC networks for the TRP and the Kv channel superfamily using the aplmDCA (Ekeberg et al., 2013). We also performed an error analysis to assess the strength and robustness of the ECs we inferred. In our framework, the statistical significance of an evolutionarily coupled pair of amino acids is given by a Z-score (ECZ score). Only ECs exceeding a chosen significance threshold of ECZ > 4 are reported in this study; considering a single-tailed test and correcting for the multiple comparison, this corresponds to a Bonferroni-corrected p-value of ∼0.1. Technical aspects of the evolutionary network inference with aplmDCA and the error analysis are described in Materials and methods.
Before delving into the details of the TRP and Kv channel EC maps, we cautiously note the conceptual differences between ECs and energetic couplings in proteins. Energetic couplings, as measured by double-cycle mutant analysis, correspond to the nonadditive energetic effects of double mutants on the stability or functional activity of a macromolecule (Yifrach and MacKinnon, 2002). ECs are instead meant to describe how knowledge of the identity of a residue at a given position affects our expectations about a different one. Although it is reasonable to expect that an energetically favorable pairwise interaction is conserved along evolution through compensatory mutations, the “strength” of the associated statistical signal may vary: interactions that are of nonpolar nature can arguably rely on a larger repertoire of residue pairs compared with polar H-bonding ones that require a more precise match.
Chains of such pairwise interactions may underlie allosteric signal communication (Yifrach and MacKinnon, 2002; Süel et al., 2003). However, unlike protein structure, which is consistently conserved among homologues, the conservation of discrete allosteric pathways has not been demonstrated to the same extent (Fodor and Aldrich, 2004).
Here, we use the presence of ECs to draw attention to physically contacting residues that have coevolved unambiguously. Coevolution implies that those physical interactions have been important during evolution and are likely to be conserved among members of the superfamily; however, we cannot make a claim on the energetics of specific molecular processes. Finally, a comment is in order concerning coevolving pairs that are not in contact in the available experimental structures analyzed in this paper: in analogy with our previous investigation of the VSD (Palovcak et al., 2014), the coevolution observed for some of these pairs could be the signature of residue–residue contacts that exist only in channel conformations distinct from those captured by the experiments so far. Exploring this fascinating hypothesis amounts to distinguishing with accuracy bona fide coevolution signals from false positives, a task that is beyond the scope of the current investigation.
ECs involving TRP S4 are consistent with a static helix
Mapping the strong and robust ECs (ECZ > 4) onto the structures of the TRPV1- and paddle chimera–sensor domains, we observe that the upper S1–S2 and lower S2–S3 interfaces exhibit numerous ECs in both families, although the couplings in Kv channels are more extensive and run further down the S1–S2 interface (Fig. 4, B and D, right).
In the TRP sensor, several couplings are also observed between the aforementioned aromatic TRP S1 positions 441 and 444 and bulky residues on S3 (TRP 516 and 523; Fig. 4 A). Such interactions would obstruct the gating pore of a Kv channel VSD and are observed neither as residue contacts nor ECs in the Kv channels.
Most significantly, we find differences in the ECs of the S4 helices of TRP and Kv channels. In previous work, we demonstrated that VSDs in general do not feature significant coevolution between the hydrophobic contact interfaces of S4 and the neighboring S1 or S3 helices but do show coevolution between the gating charges and acidic or polar residues in the gating pore (Palovcak et al., 2014). These results are consistent with a recent investigation revealing that coevolution analysis can detect residue contacts unique to distinct configurations (Morcos et al., 2013). Here, the pattern of ECs is shown experimentally by the portability of S3b–S4 paddle motif among chimeric voltage-gated cation channels, where the gating charges are conserved among the transposed paddle motifs, but the hydrophobic interfaces are not (Alabi et al., 2007). In analogy with these results, we interpret an absence of extensive networks of coevolving pairs as the signature of the mobility of the structural elements involved. This association between coevolution and structural dynamics is, at this level, a working hypothesis. We argue that because the dynamic remodeling of the contact map associated with a conformational change entails the interaction of each residue with more than one partner, correlations among residues are potentially of higher order than the pairwise ones considered here. Accessing these higher order correlations and thus testing this hypothesis is challenging because of the limited number of sequences available, and is beyond the scope of this type of analysis.
These observations are largely recapitulated when considering just Kv channel VSDs, although we see more of the gating charges coevolving with gating pore residues, and we do see one hydrophobic interaction between S4 and S1 (Kv 244–366) (Fig. 4, B and D). In contrast, the S4 helix of the TRP-sensor domain is extensively coupled at the hydrophobic faces of the S1 and S3 helices, stabilized into a putatively stationary conformation by specifically tuned interactions (Fig. 4, A and C, left). Consistently, chimeragenesis studies show that the corresponding paddle motif is apparently not portable in TRP channels (Kalia and Swartz, 2013).
The TRP sensor–pore interface has shifted through evolution
In Lee et al. (2009), sequence analysis and electrophysiology illustrated the functional importance of a coevolved contacting interface between the VSD and pore of Kv channels. This interface involved residues at the top of the S1 helix and the pore helix of a neighboring subunit. Its disruption yields channels that either do not traffic to the membrane or are defective in the final concerted step of activation gating (Lee et al., 2009; Shem-Ad et al., 2013). Despite using a dataset of Kv channels with a different size (ours is an order of magnitude larger), phylogenetic scope (ours is more confined), and algorithm for coevolution detection (theirs is a perturbation-based approach, whereas ours infers structured probabilistic graphical models), we similarly observe strong ECs at this interface in Kv channels (Fig. 5, B and C). Additionally, Kv position 248, the residue most sensitive to mutation at this interface, shows up with high relative entropy, suggesting a specific evolutionary pressure that constrains this position to small polar residues (Fig. 3 A).
In TRPV1, this sensor–pore domain interface has shifted (Liao et al., 2013). With respect to the paddle-chimera structure, the TRPV1-sensor domain is rotated counterclockwise, breaking most of the S1–pore helix interactions (Fig. 5 B). Primary contact between the sensor and the pore is made instead through the S4 helix and the S5 helix of a neighboring subunit. This rolling of the sensor domain is similarly observed in the bacterial sodium channel NavAb (Payandeh et al., 2011). Accordingly, we see no positions with high relative entropy on the top of the TRP S1 analogous to Kv 248, and we see strong EC of S4 positions to contacting S5 positions (Figs. 2 A and 5, A and B). Interestingly, the structure of TRPA1 reveals an analogous shift of the interface, thereby confirming the generality of our coevolution analysis (Paulsen et al., 2015).
The pore architecture is more variable in TRP channels than in Kv channels
Aligning the sequence logos of the TRP and Kv pore domains, we see again that despite structural similarity, the underlying distributions of amino acids observed at the same structural position greatly differ in the two subfamilies (Fig. 6 A). This difference is perhaps most striking around the pore helix and selectivity filter (TRP 627–649 and Kv 427–448). Kv channels have the common constraint of selective potassium ion conduction, which is echoed by the high relative entropy of the pore architecture responsible for this function (Fig. 6 C, the region in green). TRP channels have more variable cation selectivity depending on the channel and may exploit the pore architecture for channel-specific modulation of gating and permeation (Owsianik et al., 2006). TRP channels have comparatively low relative entropy in this region (Fig. 6 B).
In both superfamilies, the pore itself is made up by the S6 helices of the four subunits. However, when comparing TRP and Kv channels, the relative entropy of the distribution of residues in the S6 helix is high. TRP channels feature a conserved asparagine (TRP 676) common to voltage-gated sodium and calcium channels, whereas Kv channels contain a characteristic “P[VT][PG]” motif (Guda et al., 2007). As observed previously in a smaller dataset, unlike almost all voltage-gated cation channels, TRP channels do not possess a gating hinge glycine at TRP position 670 (Kv 466) (Zhao et al., 2004).
The interface between the base of S5 (TRP 560–574 and Kv 380–394, the so-called “S4–S5 linker”) and the C-terminal region of S6 shows high relative entropy in both TRP and Kv channels, although in neither case are the residues strictly conserved. This interface has been shown to be crucial for VSD–pore coupling in Kv channels and allosteric modulation of pore opening in TRP channels (Lu et al., 2002; Labro et al., 2008; Boukalova et al., 2010). Curiously, TRP position 573 in this region represents an insert present only in TRPV channels (Fig. S1).
A coevolving network at the base of S5–S6 is observed in TRP but not Kv channels
A network of ECs bridges the top interface between the S5 and S6 helices of both families (Fig. 7, C and D, left in both). This network centers on a sensor-facing aromatic residue (TRP 591 and Kv 410) but also runs behind the pore helix and throughout the pore architecture (Fig. 7, C and D). In Kv channels, these couplings mount the pore helix to S5 and S6, the selectivity filter to the pore helix, and the pore loops to each other. In TRP channels, similar features are seen but with less density.
Similarly, the S4–S5 linkers are coupled to the C-terminal portion of S6 in both TRP and Kv channels. As mentioned above, this interface plays a substantial role in coupling the voltage-dependent motion of S4 to pore opening (Lu et al., 2002; Labro et al., 2008).
In contrast, the bent region in the middle of the pore domain has several very strong intrasubunit ECs in TRP channels and none in Kv channels. One mostly acidic position, TRP 576 (Kv 395), has three such ECs. In Kv channels, 395 coevolves with one nearby aromatic residue (Kv 481) on an adjacent subunit, but this intersubunit coupling is not in physical contact per se and does not correspond to an EC in TRP channels. A recent study also noted this coevolution between Kv 395–481 and investigated the interaction in detail (Pless et al., 2013). However, the equivalent positions (and their very different interacting partners) have not been characterized in TRP channels.
No intersubunit ECs are observed within the upper pore of TRP channels
Several ECs in the pore domain were found to be in contact across subunits. In both families, numerous intersubunit ECs are observed at the lower region of the pore, connecting N- to C-terminal regions of S4–S5 linkers and the S6 helices (Fig. 7, E–H). In TRP channels, several strong intersubunit couplings were also observed between the conserved hydrophobic residues of the lower S6 helix.
Intersubunit ECs are observed at the upper region of the pore in Kv channels. Most of these couplings link the pore helix to the S6 helix of a neighboring subunit. In TRP channels, there are no such intersubunit ECs, even though the intersubunit interface occupies a similar contact area.
Evolution has encoded a register-shifting helical distortion into the S6 helix of TRP channels
Close examination of the structure of TRPV1 reveals a helical distortion in the middle of S6, where Y671 appears to have popped out from the helix and shifted the helical register of the TRPV1 S6 (the helical repeat goes from i + 4 to i + 5 after Y671; see Fig. 8 A, top inset). Because the S6 helices line the ionic permeation pathway in TRP and Kv channels, and presumably contain the gate, subtle structural features are worth investigating in detail.
Structural alignments of the TRPV1 S6 helix to the paddle chimera S6 and the bacterial sodium channel NavAb S6 indicate that position Y671 is an insertion. However, considering the sequence alignment of TRPV1 to NavAb rules this out: if TRP 671 is an insertion, it would cause the conserved N676 to be misaligned from the obviously equivalent position in NavAb (N211) (Fig. 8 B). Instead, the observed bulging residue may be caused by an inherent geometrical distortion of the TRP S6 helix. Consistent with this hypothesis, we recall the aforementioned contacting couplings at the S5–S6 interface in TRP channels. Not only are the equivalent pairs not coevolving in Kv channels, most are not even in contact (Fig. 8 A, bottom inset). It is precisely the register shift of S6 that enables these evolutionarily coupled residues to make contact in TRPV1. Although the distortion of the S6 helix is unambiguous from the TRPV1 structure and its representative electron density map, this featured was not discussed previously (Liao et al., 2013). We hypothesize that these residues have been properly selected such that favorable, highly specific interactions pin the S6 helix to the S5 helix out of register, distorting the helix and breaking the backbone hydrogen-bonding network (Fig. 8 C). This helical defect has been recently observed in TRPA1, another member of the TRP family (Paulsen et al., 2015). Strikingly, a recently determined structure of the ryanodine receptor RyR1 (PDB accession no. 3J8H) shows a similar distortion of S6 with the bulging residue (Val4926) in the corresponding helix turn (Yan et al., 2015), thereby suggesting that this mechanism could be shared with other evolutionarily related members of the 6-TM superfamily.
Comparative sequence analysis of the TRP channel superfamily has illuminated subtle and unanticipated sequence and structural features that differentiate the TRPV1 structures from those of Kv channels, including the asymmetric solvation of the TRP sensor, the shift of the intersubunit sensor–pore domain interface away from the pore helix, and the register-shifting helical distortion in S6. Such structural differences are robustly represented in the ensemble statistics of the TRP and Kv sequences, implying their significance not only in TRPV1 but also throughout the TRP superfamily. Collectively, these results suggest that TRPV1 is a good representative of the TRP TM region. This insight, coupled with our nearly comprehensive MSA of TRP channel sequences, will facilitate the accurate comparative modeling of other TRP channels while more detailed structural information is collected.
TRP channels have diverged in sequence and structure for at least a billion years; common sequence motifs and patterns of amino acid covariation present in the ensemble of their extant sequences should therefore represent the action of common, constraining evolutionary pressures. We now attempt to rationalize these observations with previously conducted experiments in a model of TRP channel structure and function.
First, we note that despite structural similarity, the TRP sensor is not likely acting as a canonical voltage sensor in any TRP homologue. Unlike Kv channels, TRP channels feature asymmetric solvation of the sensor and a gating pore plugged by conserved aromatic residues. The primary dynamic gating element of Kv channels is the positively charged S4 helix (Aggarwal and MacKinnon, 1996). The requirement for voltage-dependent mobility in the Kv S4 has largely relaxed the pressure for highly tuned hydrophobic helical interaction with the contacting S1 and S3 helices, and we see no EC there (Fig. 3, D and E). On the other hand, in weakly voltage-sensitive TRP channels, only one positive position is found on S4 (TRP 551), and several positions of the S4 helix are evolutionary coupled to positions in the S1 and S3 helices, consistent with a stationary, structural role for the TRP S4 (Fig. 4, A and B).
Also, Kv channels possess an S1–pore helix interface whose disruption yields nontrafficking or substantially modulated channels (Lee et al., 2009; Shem-Ad et al., 2013). This interface is almost entirely absent in the TRPV1 structure, as the sensor has shifted with respect to the pore domain (Fig. 4 B). The equivalent S1 positions are significantly more variable in TRP channels, and mutation of S629A in TRPV1 produces functional channels with normal capsaicin sensitivity and weakened but not abolished proton sensitivity, consistent with a more relaxed pressure on this region in TRP channels (Ryu et al., 2007). Most significantly, strong ECs observed at this interface in Kv channels are absent in TRP channels, appearing instead at the more substantial S4–S5 intersubunit interface (Fig. 4 A).
Why shift this interface away from the pore helix? If the lack of ECs on the hydrophobic interfaces of the S4 helix in Kv channels is indeed an imprint of its dynamic role in channel gating, our data suggest that the pore helix has taken up its role in TRP channels. Although the extensive coupling of TRP S4 to the S1, S3, and S5 helices strongly suggests that it is not making large movements during gating, the pore helix of TRP channels is largely decoupled with respect to the Kv pore helix, losing intersubunit ECs to the aforementioned S1 helix and the contacting S6 helix of a neighboring subunit. Two intrasubunit couplings mount the pore helix to S5 and S6 with hydrogen-bonded interactions in TRPV1 (TRP 584–641 and 637–666), and disruption of these interactions in TRPV1 (T641S) or yeast TRPY1 (where TRP 666 corresponds to TRPY1/YVC1 Y442C) yields constitutively active or completely nonfunctional channels, respectively (Su et al., 2007; Myers et al., 2008). Mutation of a conserved aromatic residue making a ternary contact with both these pairs (TRP 640) also yields a constitutively active channel (Myers et al., 2008). Additionally, structures of TRPV1 in three conformations show a different location of the pore helix in each conformation, suggesting that it is a highly dynamic element (Cao et al., 2013). Several important gating modalities including proton potentiation, tarantula toxin binding, and possibly thermal gating localize to the TRPV1 pore architecture (Ryu et al., 2007; Bohlen et al., 2010; Grandl et al., 2010), and it has also been recently suggested that the weak voltage sensitivity of TRPV1 may be caused by the movement of the intrinsic dipole of the pore helix in response to membrane depolarization (Cao et al., 2013). Considering this body of evidence and the EC data, we speculate that the dynamic movement of the pore helix, mediated by the conserved and coevolving networks of residues described above, is a fundamental feature of gating in most, if not all, TRP channels.
TRP channels lack a glycine gating hinge (Kv 466, TRP 670) that is conserved in Nav, Cav, and Kv channels (Zhao et al., 2004). Even so, the resolved conformations of TRPV1 show the S6 helix pivoting in this region to open and close the intracellular gate (Cao et al., 2013). We traced the flexibility of the TRP S6 helix to a register-shifting helical distortion around TRP 671 and identified strong, TRP-specific ECs at the joint in S5–S6 that pin S6 into this distorted conformation. Like the half-twist of a Möbius strip, the broken backbone hydrogen bonds induced by this distortion can be shifted but not fully ironed out.
Free energy changes associated with shifts in frustrated backbone hydrogen bond networks have been estimated to be very small (<1 kcal/mol), consistent with the modest free energies measured for TRP channel activation (Brauchi et al., 2004; Yao et al., 2010; Cao and Bowie, 2012). We propose that the aforementioned motion of the pore helix is translated into pore gating by enhancing the relative stability of S6 backbone hydrogen bond networks that correspond to an open intracellular gate (Fig. 8 C). Importantly, the recently reported structure of another member of the TRP family, TRPA1, reveals an analogous distortion of S6 (Paulsen et al., 2015). This observation strongly supports our fundamental hypothesis concerning a common gating mechanism and lends credibility to the specific mechanism envisioned here. For both TRPV1 and TRPA1, binding of agonists or antagonists is able to induce subtle conformational transitions in the region surrounding the helical defect, thereby promoting a displacement of some of the pore-lining residues (Cao et al., 2013; Paulsen et al., 2015). Equivalently, this conformational equilibrium could be shifted by perturbations propagating from intracellular cytosolic domains. In this model, signals from distant domains detecting diverse stimuli could be integrated into one allosteric pore-opening step. Because this argument comes largely from the ensemble statistics of TRP channel sequences and not the specific features of a single channel, we expect this mechanism may underlie polymodal gating throughout the superfamily.
Although the features identified in this study are robust measurements of the statistical properties of sequences in the TRP superfamily, this model of TRP channel structure and gating is, of course, speculative. It should be seen as an early hypothesis-generating step in a full explanation of TRP channel gating rather than its culmination. Even so, the usefulness of the approximate answers obtainable through exploratory data analysis should not be trivialized (Tukey, 1962). We expect that this analysis will serve as a solid platform for the design and interpretation of experiments conducted across the diverse and daunting TRP channel superfamily.
We thank Sandipan Chowdhury and the members of the laboratory of Kenton Swartz for insightful critiques.
This work was supported by the National Institutes of Health via National Institute for General Medical Sciences (grant P01 GM55876 to M.L. Klein) and by the National Science Foundation (grant ACI 1440059 to M.L. Klein and V. Carnevale). L. Delemotte receives funding from European Union Seventh Framework Program “Voltsens” (grant PIOF-GA-2012-329534).
The authors declare no competing financial interests.
Kenton J. Swartz served as editor.