The amino acid sequences of H. sapiens cytoglobin (Accession #Q8WWM9), neuroglobin (Accession #Q9NPG2), myoglobin (Accession #P02144), and hemoglobin subunits alpha (Accession #P69905), and beta (Accession #P68871) were downloaded from the ExPASy Proteonomic Server of the Swiss Institute of Bioinfomatics (http://www.expasy.org; http://www.uniprot.org). About 98% of the protein sequences provided by UniProtKB are derived from the translation of the coding sequences (CDS) which have been submitted to the public nucleic acid databases, the EMBL-Bank/Genbank/DDBJ databases (INSDC). Amino acid sequences were compared using the Pairwise Sequence Alignment software (LALIGN) at http://www.ebi.ac.uk/Tools/services/weblalign/ to find internal duplications by calculating non-intersecting local alignments. The Emboss Water protocol (version 36.3.5e Nov, 2012; preload8) used here employs the Smith-Waterman algorithm (with modified enhancements) to calculate the local alignment of two sequences.

PSIPRED is a simple and accurate secondary structure prediction method, incorporating two feed-forward neural networks which perform an output obtained from Position Specific Iterated—Blast. PSIPRED 3.2 achieves an average Q3 score of 81.6% (Nugent et al., 2011) and is available at http://bioinf.cs.ucl.ac.uk/psipred/. Alpha-helical transmembrane (TM) channel proteins play key roles in a variety of cellular processes essential in movement of ions and molecules across membrane lipid bilayers. As noted by Nugent and Jones (2012), technical difficulties in obtaining high quality crystals have led to an under-representation of channel proteins in structural databases. These workers have trained a support vector machine classifier to predict the likelihood of a TM helix being involved in pore formation. This approach has a prediction accuracy of 72%, while a support vector regression model is able to predict the number of subunits participating in the pore with 62% accuracy. Pore-lining regions in transmembrane protein sequences in this study were predicted using the method of Nugent and Jones (2012).

The Consensus Constrained TOPology prediction (CCTOP; http://cctop.enzim.ttk.mta.hu) server is a web-based application providing transmembrane topology prediction. In addition to utilizing 10 different state-of-the-art topology prediction methods, the CCTOP server incorporates topology information from existing experimental and computational sources available in the PDBTM, TOPDB, and TODOM data bases using the probabilistic framework of hidden Markov model (Dobson et al., 2015).

Based on X-ray crystallography data (PDB files), Pellegrini-Calace et al. (2009) have developed a computational approach (PoreWalker 1.0) for the identification and characterization of channels and cavities in proteins based on their three-dimensional structure. The software is available as a web-based resource at http://www.ebi.ac.uk/thornton-srv/software/PoreWalker/.

Meruelo et al. (2011) have identified distinct residue preferences in kinked vs. non-kinked helices and have exploited these differences and residue conservation to predict kinked helices using a neural network algorithm. The kink predictor, TMKink, is available at http://tmkinkpredictor.mbi.ucla.edu/.

Using a GST-fusion protein containing the caveolin scaffolding domain as a receptor to select peptide ligands from a bacteriophage display library, Covet et al. (1997) identified at least three related but distinct caveolin binding motifs, ΦxxxxΦxxΦ, ΦxΦxxxxΦ, and ΦxΦxxxxΦxxΦ (where Φ represents an aromatic amino acid, W, Y, or F), that have been shown to interact with caveolin in most proteins. The caveolin scaffolding domain is ascribed a role in the recognition and binding of cholesterol, which is highly concentrated in caveolae and essential for their formation. More precisely, formation of caveolae strictly requires tight binding of caveolin to cholesterol (reviewed in Covet et al., 1997). This functionality is thought to localize to a cholesterol recognition/interaction amino acid consensus motif in residues VTKYWFYR (reviewed in Fantini and Barrantes, 2013) defined as the CRAC motif. Despite its short length, this 20-residue segment appears to incorporate an array of critical functionalities in overlapping sequence elements.

CRAC is a short linear amino acid motif that mediates binding to cholesterol and stands for Choles-terol Recognition/Interaction A mino acid Consensus sequence (Li and Papadopoulos, 1998; reviewed in Fantini and Barrantes, 2013). In a C-terminus to N-terminus direction the motif consists of a polar Leu (L) or Val (V) residue, followed by a segment containing 1–5 of any residues, followed by a mandatory aromatic Tyr (Y) residue, a segment containing 1–5 of any residues, and finally a basic Lys (K) or Arg (R). In the one letter amino acid codes the algorithm is (L/V) – X1–5 – (Y) – X1–5 – (K/R). A second cholesterol recognition domain similar to the CRAC domain (CARC) has been identified (Baier et al., 2011) but exhibits the opposite orientation along the polypeptide chain (“in- verted CRAC”), i.e., (K/R) – X1–5 – (Y/F) – X1–5 – (L/V). CARC is distinct from CRAC in that the central amino acid can be either Y or F.

Transmembrane (TM) helices and pore-lining regions span a large variety of membrane systems including cell-surface receptors, redox proteins, ion channels, and transporters (reviewed in Hildebrand et al., 2006). TM helices are found to be amino acid sequences with hydrogen-bonded helical configurations, including α-, 3₁₀− and π-helices. The α-helix is very common; while the 3₁₀ helix is found at the ends of the α-helix and π-helices are rare. Based on data for 160 transmembrane helices of 15 non-homologous high resolution X-ray protein structures, Hildebrand et al. (2004) found the average TM helix to be a 17.3 ± 3.1 (SD N = 160) residues in length. However, it is difficult to obtain a representative sampling of the membrane proteins.

Many membrane proteins have proven to be difficult to crystallize owing to their partially hydrophobic surfaces, flexibility, and lack of stability as well as physiologically bound lipids (reviewed in Carpenter et al., 2008).

Nugent et al. (2011) have developed a method (MEMSAT-SVM) to predict pore-lining regions in transmembrane (TM) proteins and to distinguish pore-lining regions from TM helices (Nugent and Jones, 2012). Pore-lining regions are usually enriched in negatively (e.g., E, D) or positively charged residues (e.g., K, R, H). Four features were used to train a support vector regression model: sequence length, the number of pore lining residues, with the target value set to the number of subunits contacting the pore within the membrane region. Nugent et al. (2011) were able to identify pore-lining regions and to predict both the likelihood of transmembrane helices involved in pore (channel) formation and to determine the number of subunits required to form a complete pore or channel. As shown in Figure 1, topological analyses of H. sapiens Cygb (top), Ngb (upper middle), Mb (middle), Hbb (subunit β, lower middle), and Hba (subunit α, bottom) predict that, with the exception of Hba, each globin contains a 16 residue pore-lining (channel-forming) structure in the C-terminal region of the protein (blue squares). In contrast, Hba contains a 15 residue TM helix (black squares). Typically, channel proteins contain a cavity (or pore) which spans the entire membrane protein with an opening on each side of the membrane. The pore often runs parallel to TM helices, forming a path along which ions or molecules travel, with adjacent structural features determining pore specificity. As shown in Figure 1, the N-terminal region of both Cygb and Mb was predicted to be extracellular (orange boxes) whereas the N-terminal region of Ngb and the Hba and Hbb subunits is cytoplasmic (white boxes).

Figure 2 compares the topology for Hba (upper panel) and Hbb (lower panel) using the Consensus Constrained TOPology server (Dobson et al., 2015). Colors are based on the localization of peptide sequences: gray (transit sequence), black (signal peptide), blue (extracytosolic), red (cytosolic), yellow (transmembrane), and orange (re-entrant loop regions). CCTOP analysis of Hba predicts a TM helix in the 100–120 residues region, consistent with that shown by the MEMSAT-SVM algorithm in Figure 1. Similar TM helices were predicted by the Octopus, Prodiv, and Scampi algorithms. In contrast, CCTOP analysis of Hbb predicted a helix-turn-helix in the 20–56 residue region but failed to detect pore-lining regions in the C-terminal region. The findings in Figures 1, 2 suggest that Cygb, Ngb, Hba, and Hbb, as well as Hba (Shaklai et al., 1977), may also be integral to cell membranes.

Cybg, Ngb, Mb, Hba, and Mbb amino acid sequences were compared using LALIGHN (see Section Methods). As shown in Table 1, a comparison of H. sapiens Cybg and Ngb indicates a 29.4% identity (63.4% similar) in a 154 amino acid overlap (19–171: 3–152) with a Waterman:Eggert score of 251. A similar comparison of Cybg and Mb indicates a 28.9% identity (63.1% similar) in a 149 amino acid overlap (19–167: 3–142) with a Waterman:Eggert score of 227. A comparison of Mb and Hba indicates a 39.3% identity (46.4% similar) with an amino acid overlap of only 28 amino acids (17–39: 57–84) with a Waterman-Eggert score of only 41. Finally, a comparison of Hba and Hbb indicated a 25.5% identity (58.6% similar) in a 145 amino acid overlap (4–146:3–147) with a Waterman-Eggert score of 161. Despite the limited amino acid sequence identities and similarities shown here, all four clearly belong to the globin superfamily (Pesce et al., 2002), all contain the typical globin fold (Laberge and Yonetani, 2007) and reversibly bind O2, CO, and NO (reviewed in Pesce et al., 2002). Hb and Mb have been considered almost exclusively as transporters of metabolic gases and nutrients whereas erythrocytes also take up and inactivate endothelium-derived NO to form methemoglobin and nitrate. As noted by Brunori et al. (2005), Ngb is unlikely to be involved in O2 transport (like Mb) but has been proposed to act as a sensor of the O2/NO ratio in the cell, possibly regulating the GDP/GTP exchange rate, forming a specific complex with the Gαβγ–protein when oxidized, but not when bound to a gaseous ligand (reviewed in Chu et al., 2008). Tomita et al. (2009), using pulsed- laser at cryogenic temperatures to study CO migration in Mb found that the migration of the CO molecule into each cavity induces structural changes of the amino acid residues, resulting in the sequential motion of the ligand, and the cavity. They propose that a breathing motion of internal cavities is a general mechanism of ligand migration.

Pellegrini-Calace et al. (2009) have developed an improved computational approach (PoreWalker 1.0) for the identification and characterization of channels and cavities in proteins, based on their three-dimensional structure. Given a set of 3D crystallography coordinates, this method can detect and identify the pore centers and axis using geometric criteria and the biggest and longest cavity/channel through the protein is identified. Pore features, including diameter profiles, pore-lining residues, the size, shape, and regularity of the pore are used to provide a quantitative and visual characterization of the channel. Since pore-forming proteins interact with lipid bilayers to generate “proteolipidic” pores (Gilbert et al., 2014), the “crystal structure” of a lipid-free protein crystal may not correspond precisely to the physiological cavity/channel.

The left column of Figure 3 illustrates the PoreWalker structure profiles for Cygb (top), Ngb (upper middle), Mb (lower middle), and Oxy T State Hb (bottom), whereas the right column presents the corresponding pore diameter profiles. The structure profiles are visualizations of a pore section showing pore-lining residues and pore centers, using 3À steps: the section of each structure was obtained by cutting the protein structure along the xy-plane, where the x-axis corresponds to the pore axis, and the y-coordinate >0 only are displayed. The PoreWalker graphics identifies and characterizes the biggest and longest cavity in a set of 3D crystallography coordinates in a cholesterol and lipid-free preparation. The pore-lining atoms and residues are colored in orange and blue, respectively. The remainder of the protein is shown in green. Red spheres represent pore centers at given pore heights and their diameters correspond to 1/10 of the pore diameter calculated at that point. As can be seen in the right hand column, a cavity of varying diameters extends across all four globins. Cygb (1UMO) and Ngb (1OJ6) each contain a sequence of cavities with varying diameters, whereas Mb (3RGK) exhibits a large single cavity localized left (Mb) whereas Oxy Hb (1GZX) contains a large centrally located cavity(s).

Maximum cavity diameters varied from about 6 to 11 angstroms in the 4 globin crystals, or about 3–5 times the diameter of a water molecule. The cavities and channels may be larger in fully hydrated proteins embedded in lipid bilayers.

The shortest sequence of amino acids in proteins that contains functional and structural information is termed a “motif,” whereas a conserved part of a given protein that can evolve, function, and exist independently is termed a “domain.” Domains and motifs contribute to the structural basis of the physiological functions of proteins and each domain can be considered as a semi-independent structural unit of a protein capable of folding independently (e.g., Batey et al., 2008; Dill et al., 2008; Pang et al., 2008). A variety of different methodologies have been employed to predict domains but many are fraught with problems since domain assignment is difficult even when the structures are known. Bryson et al. (2007) have developed a useful method for computer-assisted protein domain boundary prediction, using the DomPred server (see Section Methods). This server uses the results from two completely different categories (DPS and DomSSEA). Each is individually compared against one of the latest domain prediction benchmarks to determine their respective reliabilities.

Figure 4 illustrates a DomPred analysis of the five different globins from H. sapiens. As can be seen, the protein domain patterns differ, with the domains of Hba and Hbb being the most complex. As noted by Marsden et al. (2002), the modular nature of many domains means that they can often be found in proteins with a similar content, but with different orders, or in different proteins in combination with entirely different domain structures. In addition to the TM helices, pore-lining regions, and domain structures described above, functional domains such as CRAC or CARC cholesterol binding motifs (reviewed in Fantini and Barrantes, 2013) can be identified within the Cygb, Ngb, Mb, and Hb subunits α (Hba) and β (Hbb). Figure 5 compares the distribution of cholesterol-binding (CRAC/CARC) domains (highlighted in red) with caveolin-binding motifs (underlined in blue) and both TM helices (double underlined) and pore-lining regions (single underlined) identified. A cholesterol binding domain similar to CRAC has also been characterized and is termed a “CARC” motif (reviewed in Fantini and Barrantes, 2013). The so- called “CARC motif” is comparable to the CRAC domain but exhibits the opposite orientation along the polypeptide chain (“inverted CRAC”; Baier et al., 2011) and is distinct from CRAC in that the central amino acid can be either Y or F.

Molecular modeling studies have shown that CRAC/CARC motifs have a good fit for cholesterol (cit. Fantini and Barrantes, 2013). Although cholesterol is concentrated in sphingolipid-enriched membrane microdomains such as “lipid rafts” (Silvius, 2003), it is also present in the lipid disordered phase of the plasma membrane that contains high amounts of glycerolphospholipids such as phosphatidylcholine (Heberle and Feigenson, 2011). Artificial lipid membranes also contain rather large voids, the volume of which can be an order of magnitude larger than the largest spherical cavities present in biological membranes (Alinchenko et al., 2005). Studies have shown that addition of cholesterol to artificial lipid membranes produces elongated voids lying perfectly parallel with the membrane in normal axis (Alinchenko et al., 2005), reminiscent of the cavities shown for globins (see Figure 3).

Since the PoreWalker algorithm identifies both the amino acid residues and/or domains lining the largest cavity within each globin, it is possible to identify potential cholesterol binding (CRAC/CARC) sites lining the globin cavities. Caveolin (blue underline in Figure 5) is an unusual protein that can exist both as an integral membrane protein and as a soluble protein in multiple cell compartments (reviewed in Liu et al., 2002). As shown, Ngb and Hbb each contain a caveolin binding motif (⁴²FQYNCRQF⁴⁹ and ³⁶YPWTQRFFESF⁴⁶) overlapping CRAC/CARC motifs. Topological analysis of all three isoforms of H. sapiens caveolin indicates that each caveolin molecule also contains a pore-lining region, a TM helix-turn-TM helix, and a CRAC/CARC motif (Morrill et al., 2014). Thus, the presence of CRAC/CARC motifs as well as pore-lining regions within the caveolins in proximity to either the pore-lining regions or TM helices in Ngb and Hba could contribute to channel formation.

As shown in Figure 5, two to six CRAC/CARC domains are associated with the Cybg, Ngb, Mb, Hba, and Hbb monomers. Only Cybg (top) and Hbb (bottom) exhibit pore-lining regions with minimal (one residue) overlapping CRAC/CARC motifs. The Cybg dimer would contain up to 6 bound cholesterol molecules whereas the Hb tetramer would contain up to 14 bound cholesterols. Judging from the topology for Cybg and Hbb shown in Figure 1, the CRAC/CARC motifs are at the protein surface and may involve surface docking sites (e.g., Listowski et al., 2015). However, since removal of cholesterol and/or other lipids is often essential in obtaining protein crystals (reviewed in Carpenter et al., 2008), the protein data base for globin preparations used in x-ray crystallographic analysis list only globin, heme and sulfate or phosphate ions. Cholesterol, phospholipids and water in kinks have been removed.

Similarly, it should be noted that Hb used in biophysical studies (e.g., Roche et al., 2013) has been prepared by anion exchange chromatography methods (Andrade et al., 2004) that largely removes endogenous lipids and may thus alter the physiological properties of the Hb molecule.

It should also be noted that methionine is the N-terminal amino acid in each of the 5 globins shown in Figure 5. In a biological system, protein synthesis is initiated universally with the amino acid methionine (see Drabkin and RajBhandary, 1998) which is subsequently removed by enzymes in the cell. Since x-ray crystallography methods require large quantities of a pure protein (Watanabe et al., 2010), crystallography pdb files of cell-free preparations such as those analyzed by PoreWalker (Figure 3) contain N-terminal methionine. The codon initiator methionine is therefore listed in Figure 5, although it may not be present physiologically. A contribution from N-terminal methionine to the crystallographic structure has not been ruled out.

Most plasma membrane and secretory proteins contain one or more carbohydrate chains and the addition and subsequent processing of carbohydrates is the principal chemical modification of such proteins (Cit. Moremen et al., 2012). Glucose reacts non-enzymatically with α- and ε-amino groups on proteins in vitro and was shown to occur in vivo with the discovery of a naturally existing human Hb component, Hb A1c (Shapiro et al., 1980). Hb glycosylation is elevated in patients with diabetes mellitus, reflective of the increased blood glucose levels. The major Hb glycosylation sites in vitro and in vivo include β-Lys(K)-66 and β-Lys(K)-120. As shown in Figure 5, after correcting for the initiator methionine, Lys⁶⁶ and Lys¹²⁰ occur in two of the four Hbb CRAC/CARC motifs. This suggests that insulin directly affects cholesterol binding to the Hb β subunit. Mb is similar to the β-chain of hemoglobin (Hbb) and has a sufficient number of ε-amino groups to be glycosylated like Hbb. Banerjee and Chakraborti (2014) have shown that glyoxal, a highly reactive oxoaldehyde, is elevated in diabetes and modifies Lys¹⁴⁵ (a residue within one of two Mb CRAC/CARC motifs), suggesting that insulin may also effect Mb function.

A link between insulin and cholesterol binding is also demonstrated by the finding that, under hypoglycemic conditions, hyperinsulinism causes an increase in cholesterol binding to Hb (Tomasevic et al., 2003). Using a fluorescent lipid analog, these workers have shown that lipid internalization occurs in erythrocytes in response to elevated insulin levels. A further link between glycosylation and diabetes is the report that glycosylation per-se leads to blood pressure reduction in type 2 diabetics patients untreated for hypertension (Cabrales et al., 2008).

There is a long history concerned with Hb binding sites on the red cell membrane (see Shaklai et al., 1977). Current evidence indicates that anion exchange protein “band 3,” a multipass membrane protein, is the major membrane binding site for Hb in the human erythrocyte (e.g., Chu et al., 2008). The bicarbonate/chloride exchanger “band 3” (Anion Exchanger, AE1) is an important structural protein in the erythrocyte cell membrane (van den Akker et al., 2010). Each human erythrocyte is estimated to contain 1 million copies of “band 3.” “Band 3” is also thought to play an important role in transmembrane gas exchange and function as a point of attachment for the cytoskeletons that maintain the osmotic properties of the cell (cit. van den Akker et al., 2010). The oxygenation state of Hb also plays an important role, since deoxyhemoglobin, but not oxyhemoglobin, binds “band 3” reversibly with high affinity (van den Akker et al., 2010).

Figure 6 demonstrates that H. sapiens “band 3” (Accession #P48751) contains a large central cavity (Figures 6A,B) as well as 12 transmembrane helices located in the C-terminal region. The 705 N-terminal residues are extracellular with the subsequent 491 residues containing 12 closely packed helices and/or pore-lining regions. The five helices indicated by blue squares (2, 4, 6, 8, 10) are predicted to be pore-lining regions (Figure 6C). It should also be noted that about half of all transmembrane helices contain bends and other deviations often referred to as “kinks” (reviewed in Meruelo et al., 2011). As proposed by Meruelo et al. (2011), distortions in helix geometry may facilitate conformational changes required for protein function by providing sites of flexibility and can be important for positioning key residues precisely in the protein structure. Kinks that open the polar backbone to alternative hydrogen bonds often attract water molecules, thus providing a polar region within the hydrophobic core (Hall et al., 2009). TMkink analysis indicates that TM-1 of “band 3” contains a large kink (⁷¹⁰A – A⁷²³), thus providing a potential polar-water filled region. Anion exchange protein “band 3” thus contains numerous pore-lining regions as well as a water-containing kink suggesting that the anion exchange protein (“band 3”) may merge with pore-lining regions in hemoglobin molecules to form bicarbonate/chloride exchange channels (see Hildebrand et al., 2006; Nugent and Jones, 2012 for discussion of channel/pore formation).

The regularity of the pore cavity is deduced from the positions of the pore centers calculated by PoreWalker, aimed at the optimization of the pore axis. If all the pore centers within a channel are co-linear, then the pore is linear and the cavity is symmetrical. If a pore centers in only certain areas of the channel, but not the entire channel, and are co-linear, the cavity will be non-symmetric, and the channel irregular. Figure 7 compares the protein orientation of Cybg (top), Ngb (upper middle), Mb (lower middle), and Hb (bottom) for the XZ plane section when Y > 0 coordinates only. Cybg exists as a dimer, Ngb and Mb are monomers and Hb is a 4-mer (ABAB). As noted above, cholesterol and other lipids have been removed prior to crystal formation. The lowest coordinate is along the pore axis (=x-axis) at the bottom. The red spheres represent pore centers at 1 À steps along the pore axis. Cavities traverse each globin, indicating that all four globins without bound cholesterol contain a single major continuous pore or channel.

Table 2 identifies distal and proximal Fe sites for Cybg, Ngb, Mb, and Mbb at the top of each column and underlines the cavity-lining Fe binding sites within Cybg (¹¹³H, heme proximal ligand), Ngb (⁹⁶H, heme proximal ligand), Mb (⁶⁵H, heme distal ligand), and Hbb (⁶⁴H, heme distal ligand). The cavity-lining Fe binding site for Hba was ⁸⁸H, heme proximal ligand (data not shown). Cavity-lining residues are defined as amino acids contributing at least one atom to the inside surface of the pore. The pore surface can be considered as a continuum made up of horizontal and vertical layers of atoms. As also shown in Table 2, the cavity-lining regions involve 57 of 190 Cybg residues, 38 of 151 of Ngb residues, 55 of the 154 Mb residues, or 35 of 144 residues in Hbb (subunit β). Residues associated with cholesterol binding (CRAC, CARC) motifs are highlighted in red. As indicated, each globin contains one or more bound cholesterol molecule within the predicted cavity and/or channel. For comparison, Table 3 compares the specific residues within the pore-lining regions that correspond to the cavity-lining residues identified by PoreWalker. Ten pore-lining residues within Cybg correspond to cavity-lining residues identified by PoreWalker. Similarly, 9 pore-lining residues are within the cavity of Ngb, and 10 and 8 residues were within the cavities of Mb and Hbb, respectively. Many of the residues identified by PoreWalker as pore-lining residues were confirmed as accessible to water (hydrophilic sites) using a water accessibility predictor (Lingwood et al., 2009) for α-helical transmembrane proteins (MPRAP).

Compared to the dimeric (Cybg) or monomeric forms (Ngb and Mb), Hb is a heterotetramer consisting of two α chains and two β chains. Table 4 compares the predicted cavity-lining residues of the T state Hb heterotetromer (1GZX). Columns A and C represent the two α subunits whereas B and D correspond to the two β subunits. The CRAC/CARC motifs are highlighted in red and the pore-lining or TM helix residues are highlighted in blue. Table 4 indicates that all four tetramers contribute pore-lining and/or TM helix residues as well as cholesterol-binding (CRAC/CARC) domains to the major cavity within Hb. As shown, each heterotetramer contributes a total 29 blue-highlighted (pore-lining and TM helix) residues and 30 red-highlighted (cholesterol-binding) residues to the largest Hb cavity.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.