In silico prediction of structure and function for a large family of transmembrane proteins that includes human Tmem41b

Pfam database screening

Searches using the sequences of VTT domain proteins Tmem41b, Yqja, Ydjx, Ydjz, Tvp38 and Mt2055 were made against the Pfam-A_v32.0 (RRID:SCR_004726) (El-Gebali et al., 2019) database using the HHPred (RRID:SCR_010276) v3.0 server (Zimmermann et al., 2018) with default parameters (-p 20 -Z 10000 -loc -z 1 -b 1 -B 10000 -ssm 2 -sc 1 -seq 1 -dbstrlen 10000 -norealign -maxres 32000 -contxt /cluster/toolkit/production/bioprogs/tools/hh-suite-build-new/data/context_data.crf) and eight iterations for MSA generation in the HHblits (Remmert et al., 2012) stage.

Contact map predictions

The DeepMetapsicov v1.0 server (Kandathil et al., 2019) was used to generate contact predictions with ConKit v0.12 (Simkovic et al., 2017b) utilised to visualise the contact maps. ConPlot (RRID:SCR_019216) was used to overlay additional prediction data (Sánchez Rodríguez et al., unpublished work).

Other prediction data

Transmembrane helical topology predictions were obtained from the Topcons server (Tsirigos et al., 2015). Secondary structure predictions were made employing a local installation of PSIPRED (RRID:SCR_010246) v4.0 (McGuffin et al., 2000). ConKit was also used to predict and visualise potential structural domain boundaries (Rigden, 2002; Simkovic et al., 2017a). Residue analysis of putative amphipathic regions were performed using HELIQUEST (Gautier et al., 2008) to determine the presence, direction and magnitude of any hydrophobic moment. Residue conservation was determined using the Consurf server (Ashkenazy et al., 2016).

Dataset for custom re-entrant database

A library of re-entrant loop sequences together with the putative re-entrant loop sequences from the query proteins were clustered to establish any visible relationships of the sequences. The library was built by obtaining a non-redundant set of 56 re-entrant helix sequences by first retrieving all 714 TM proteins that contain at least one re-entrant loop from the PDBTM (RRID:SCR_011962) (Kozma et al., 2013) and removing redundancy with a 40% identity threshold. The resulting 127 protein structures were split into their component chains, eliminating any chain lacking a re-entrant loop. The subsequent set of 188 unique re-entrant loop sequences were then filtered removing any sequences of less than 10 residues and more than 20, thereby ensuring the collection of sequences conformed to the length of typical (Yan & Luo, 2010) re-entrant loops. The remaining 56 sequences were clustered, supplemented by candidate re-entrant sequences from the proteins studied here. Clustering was performed using CLANS v1.0 (Frickey & Lupas, 2004) with the BLAST results (p-value cut-off threshold of 0.1) (Altschul et al., 1997) used to calculate strengths of similarity.

Model building

Ab initio models were built using the trRosetta (Yang et al., 2020) server with default settings. Conservation was mapped on to the models using the ConSurf server (Ashkenazy et al., 2016). Visualisation of models was achieved using PyMOL (RRID:SCR_000305) v2.3.0 (DeLano, 2002).

Structural alignments

Dali (RRID:SCR_013433) v4.0 (Holm & Laakso, 2016) was used to structurally align the output models and to query against the PDBTM (Kozma et al., 2013).

An earlier version of this article can be found on bioRxiv (doi: https://doi.org/10.1101/2020.06.27.174763)

Sequence comparisons suggest Pfam families PF09335 and PF06695 are related

HHpred (Zimmermann et al., 2018) was used to screen a selection of VTT proteins against the Pfam database (El-Gebali et al., 2019). Hits were observed in the same region against both PF09335 and the Pfam domain PF06695 (‘Sm_multidrug_ex’) which is strongly indicative of homology: a probability of 99.4% with an E-value of 9E-17 for the PF09335 hit and 98.3% and 2E-10 respectively for PF06695. A HHpred search against the Pfam database using a member of PF06695 - the short archaeal sequence Mt2055 (UniProt code W9DY28) (Apweiler et al., 2004) - returned similar results (Table 1). The Mt2055 sequence originates from the unpublished draft genome of the archaebacterium Methanolobus tindarius DSM 2278. For many of the subsequent analyses, the shorter archaeal sequence was used initially but the clear homology among this set of proteins means that inferences can be drawn across the group.

There are no known experimental protein structures representing PF09335 or PF06695, but both Gremlin and DMPfold have constructed ab initio models for these Pfam domains (Greener et al., 2019; Ovchinnikov et al., 2017).

The predicted Pfam domains are inconsistent with a structural domain

Analysis of the HHpred results obtained for the archaeal protein Mt2055 revealed the presence of additional hits for both PF06695 and PF09335 Pfam domains, in which the C-terminal half of the domains aligned with the N-terminal half of the Archaea protein. For example, residues 1-69 of the archaeal protein aligned with residues 52-117 of the Pfam PF09335 profile with a probability of 74.15%. Interestingly, contact density analysis (Rigden, 2002; Sadowski, 2013) supported the existence of a domain boundary around residue 60, in broad agreement with the HHpred results (Figure 1). Both the HHpred and contact density results therefore pointed to a specific domain structure being present.

Figure 1. Mt2055 domain analysis.

(a) Contact density profile constructed by ConKit (Simkovic et al., 2017b) utilising DeepMetaPSICOV contact prediction. Solid black line represents contact density and dotted red lines mark density minima corresponding to possible domain boundaries. (b) HHalign alignments for the N-terminal and C-terminal Mt2055 halves, formatted using Jalview (Waterhouse et al., 2009) and coloured according to the ClustalX scheme. Red bars represent helical secondary structure. (c) Maps of predicted contacts generated by DeepMetaPSICOV and plotted using ConKit; left is N-terminal half (residues 1-84) and right is C-terminal half (residues 85-168). Black points represent predicted intramolecular contacts.

Sequence & contact prediction map analysis indicate that PF06695 is made up of a tandem repeat

When the Mt2055 sequence was split at residue 60-61, the resulting N-terminal region of 60 residues and the C-terminal section of 79 residues could be aligned using HHalign (Soding, 2005) with a 78% probability and an E-value of 1.9E-3. Examination of the map of predicted contacts for Mt2055 reveals features that are present in both the N- and C-terminal halves of the protein (Figure 1c). Taken together, these data strongly support the existence of a tandem repeat within the Mt2055 protein and hence across the PF06695 and PF09335 protein families.

Interestingly, however, an equivalent sequence analysis with HHpred of other PF09335 homologues including Tmem41b itself does not reveal a repeat. However, inspection of their corresponding predicted contact maps does reveal features repeated when N- and C-halves of the protein are compared (Figure 2). Apparently, evolutionary divergence has removed all trace of the repeat sequence signal in bacterial and eukaryotic proteins, although the feature remains visible by evolutionary covariance analysis.

Figure 2. Tmem41b Contact map constructed using DeepMetaPSICOV and plotted using Conkit.

The highlighted areas represent repeat units that have been revealed through evolutionary covariance analysis.

Ab initio modelling of Mt2055 reveals an unusual topology

Several authors have deposited structures of uncharacterised Pfam families in databases (El-Gebali et al., 2019); however, Pfam domain boundaries for PF09335/PF06695, which define the limits of these previous modelling exercises, do not reflect the conserved structural domain that we predict. Given the fact that the available ab initio models were inconsistent with the transmembrane helix, secondary structure and contact predictions (data not shown), we constructed our own models of Mt2055 as well as Tmem41b and Yqja with trRosetta.

The Mt2055, Tmem41b and Yqja models had estimated TM scores from the trRosetta server of 0.633, 0.624 and 0.635 respectively, suggesting that they were likely to have captured the native fold of the family. All-against-all pairwise structural superposition of the models with DALI gave a mean Z-score of 11.9 confirming their strong similarity. We also used satisfaction of predicted contacts to validate the models (Figure 3) (Simkovic et al., 2017a). This showed that 80% of the top L predicted contacts (where L is the length of the protein) are satisfied by the model contacts for both Mt2055 and Yqja and a value of 60% was achieved for Tmem41b suggestive of good quality models (de Oliveira et al., 2017).

Figure 3.

(a) trRosetta model of MT2055 - amphipathic helix (green) and a re-entrant loop (orange) packed with a TM helix (red) (b) Superposition of DMP predicted contact map for Mt2055 and contacts from the Mt2055 model. Black points are matching contacts, red are mismatches and grey are contacts predicted but not present in the model. Diagonal is a visual representation of transmembrane helix and secondary structure prediction – central diagonal is the visualisation of the TopCons transmembrane prediction (orange being a TM helix) and the outer diagonals are the visual representation of the PsiPred secondary structure prediction (pink – alpha helix and yellow – coil). Red boxes highlight the re-entrant loop and TM helix packing. c) trRosetta model of Tmem41b only showing the conserved structural domain (residues 39-217) d) trRosetta model of Yqja only showing the conserved structural domain (residues 14-176). e) Proposed topology for (extended) VTT domain.

The models (Figure 3) contained interesting features: two inversely symmetrical repeated units each possessing an amphipathic helix (green) and a re-entrant loop (orange) packed with a TM helix (red).

The presence of a re-entrant loop packed against each TM helix can also be seen on predicted contact maps for these proteins (Figure 3b). Interestingly, each of the re-entrant helices is predicted as a single transmembrane region in the TopCons predictions (see the diagonal of Figure 3b) with a two-residue region of coil in the centre. Such a prediction would more obviously be treated as indicative of some kind of kink in the helix (Law et al., 2016) but the explanation here is these regions form re-entrant helices. Similar contact map features, indicative of re-entrant loops packing against TM helices, can be seen clearly on the contact maps of other VTT proteins (data not shown).

The analysis was performed by HELIQUEST (Gautier et al., 2008) which constructed helical wheel diagrams and provided a quantitative measure of the hydrophobic moment for the region being analysed (Figure 4).

Figure 4. Helical wheel diagrams generated using the HELIQUEST server.

Hydrophobic residues are shown in yellow, serine and threonine in purple, basic residues in dark blue, acidic residues in red, asparagine and glutamine in pink, alanine and glycine in grey, histidine in light blue and proline in green circles. Arrows represent direction and magnitude of the hydrophobic moment and residue marked with ‘N’ is the N-terminal end of the putative amphipathic helix with the residue marked ‘C’ being the C-terminal end. (a) Mt2055 putative amphipathic helix 1 (hydrophobic moment of 0.298). (b) Mt2055 putative amphipathic helix 2 (hydrophobic moment of 0.546). (c) Tmem41b putative amphipathic helix 1 (hydrophobic moment of 0.471). (d) Tmem41b putative amphipathic helix 2 (hydrophobic moment of 0.420). (e). Yqja putative amphipathic helix 1 (hydrophobic moment of 0.295). (f) Yqja putative amphipathic helix 2 (hydrophobic moment of 0.396).

Mapping conservation onto the models, using the Consurf server, indicates that the re-entrant loops are highly conserved and therefore likely to be functionally and/or structurally important (Figure 5).

Figure 5.

trRosetta models with Consurf conservation mapping for (a) Mt2055 (b) Tmem41b (c) Yqja. Conservation is shown as a spectrum from blue (highly conserved) to red (not conserved).

Re-entrant loops are also present in Cl^-/H⁺ Antiporters

The presence of re-entrant loops and the high density of conserved residues within them caused us to examine experimentally characterised re-entrant loops in the PDBTM database. A total of 56 non-redundant re-entrant helices were identified (see Methods). All 56 were clustered with the putative re-entrant loops from Mt2055 and four PF09335 homologues (Tmem41b, Tvp38, Ydjx and Ydjz) using relative E-values derived from an all-against-all BLAST run in CLANS (Frickey & Lupas, 2004) with a 0.1 p-value cut-off. The largest cluster contained 14 sequences, of which four were putative re-entrant sequences from the query proteins (Mt2055 C-terminal re-entrant, Ydjx C-terminal re-entrant, Ydjz N-terminal re-entrant & Ydjz C-terminal re-entrant), seven (3org, 5tqq, 3nd0, 3det and 6coy) were re-entrant loop sequences from Cl^-/H⁺ antiporters, one was from a boron exchanger (5l25), one from an electron transporter (2n4x) [albeit classified as a member of the lysine exporter superfamily (Saier et al., 2016)] and one from a mechanogated channel (5z10).

Analysis of the Cl^-/H⁺ antiporter structures show that they contain a similar inverted repeat as we infer for the VTT homologues, resulting in pseudo-2-fold axis of symmetry running along the membrane (Duran & Meiler, 2013). Again similarly, the Cl^-/H⁺ antiporter 3orgA also contains the amphipathic helices on the N-terminal side of the re-entrant loops. The fact that the presence of the amphipathic helices is restricted only to 3orgA and not found in all homologues suggest that these features are not essential for function (Figure 5).

A possible antiporter role for VTT proteins

The presence of re-entrant loops in a transmembrane protein strongly indicates a transporter or pore functionality since this structural feature has, hitherto, only been found in proteins of this kind (Yan & Luo, 2010). The structural similarities between the VTT proteins and the Cl^-/H⁺ antiporters raise the possibility that the families studied here are, in fact, unsuspected distant homologues having this putative pore feature in common. In that regard it is relevant to recall a hypothesis that DedA proteins are H⁺ antiporters resulting from SDM experiments (Justice et al., 2016).

A recent study has identified key residues (Figure 6) in the E. coli DedA protein Yqja that, when replaced in site directed mutagenesis experiments, resulted in decreased proton motive force across the E. coli inner membrane (Panta et al., 2019). Highlighting the essential residues (E39, D51, R130 and R136) on the Yqja model show that they come together in three-dimensional space with the N-terminal side of the first re-entrant possessing E39 and the C-terminal side possessing D51. R130 and R136 are similarly positioned on the second re-entrant loop (Figure 7). Re-entrant loops are known to form pores and here we have two proton-titratable residues (E39, D51) in close proximity to essential basic residues (R130 and R136) within a putative pore. This three-dimensional arrangement of key residues could serve a role in the coupling of the protonation status with the binding of a yet to be characterised substrate as is postulated for the multi-drug H⁺ antiporter MdfA (Heng et al., 2015) where these same residues are located inside a central cavity.

Figure 6. Essential residues determined by SDM experiments highlighted in pink on the YqjA model.

Figure 7.

(a) Left - Predicted Contact map with repeating units highlighted in yellow boxes, contact map signature of re-entrant loop packed with TM helix in red boxes.; Right - The Experimental Contact map obtained from the PDB structure with repeating units highlighted in yellow boxes, contact map signature of re-entrant loop packed with TM helix in red boxes. (b) Actual 3orgA topology; grey: TM Helices that are additional to the core; red: TM helices contributing to the formation of the core; orange; re-entrant loops contributing to the formation of the core; green: amphipathic helices contributing to the formation of the core. (c) The 2-fold pseudo symmetry of the amphipathic/re-entrant loop/TM helix core inverted repeat structure of 3orgA with membrane positions shown as grey planes obtained from PDBTM.