All human canonical sequences have been obtained from UniProt. The full list of UniProt IDs and the corresponding gene names are listed in the Supplementary Table S1.

Data for the genes was downloaded from GnomAD versions 3.1.1 and 2.1.1. We used the former for the analysis and the latter portion of the data for cases of variants for which we could not find an amino acid match in the reference sequences. The data was fetched manually and then analyzed using Python 3 (Van Rossum and Drake Jr, 2009) and Pandas (McKinney, 2010; team, T.p.d., 2020). As the variants have a wide spread of screened alleles we filtered the outliers (Supplementary Figure S2). Next, we filtered out the variants, which have the frequency lower than 0.001%. Furthermore, we retained only missense variants. Those we transferred onto a reference numbering based on an alignment of all reference protein sequences.

All experimental structures analyzed in this study have been taken from the PDB (Supplementary Item S1) (Unwin and Fujiyoshi, 2012; Miller et al., 2017; Basak et al., 2018a; Basak et al., 2018b; Chen et al., 2018; Gharpure et al., 2019; Noviello et al., 2021; Zhang et al., 2021; Zhu and Gouaux, 2021).

The molecular structures of the majority of the 168 investigated drugs (Supplementary Item S2) were downloaded from DrugBank (https://go.drugbank.com/) (Wishart et al., 2018) (154 in total). 12 of the 168 drugs were not available on DrugBank and have been retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) (Kim et al., 2016) instead. The remaining two drugs—ivermectin and picrotoxin represent receptor-bound ligand structures and were extracted from the PDB complexes 5VDH and 6X40 (Berman et al., 2000), respectively. Structures of drugs taken from PubChem and DrugBank were all downloaded and stored as SD-files containing 2D atom coordinates. This format was chosen because several of the investigated drugs get administered as mixtures of stereoisomers and information about undefined stereocenters is not present in SD-files containing 3D atom coordinates. The individual drug SD-files were then concatenated to a single SD-file for further processing by the software Flipper (https://docs.eyesopen.com/applications/omega/flipper.html, version 3.1.1.2, with default settings for all non-mandatory options) in order to enumerate all possible stereoisomers of drugs having one or more undefined stereocenters. The resulting SD-file contained a total of 246 structures which were then subjected to molecule conformer ensemble generation using the program OMEGA (http://www.eyesopen.com/, version 3.1.1.2, with default settings for all non-mandatory options). For the generation of the pharmacophore screening database, a simple KNIME workflow (Berthold et al., 2008) using nodes from the LigandScout KNIME extension (Wolber and Langer, 2005) (version 1.8.0, http://www.inteligand.com/download/LigandScout_Knime_Nodes.pdf, http://www.inteligand.com/ligandscout3/) was created which consisted of a SDF-reader node connected to a LDB-writer node. In the SDF-reader node configuration dialog the multi-conformer SD-file generated by OMEGA was chosen as input file, molecule conformer detection got enabled and “Compare by Name” has been selected as associated conformer matching strategy. After execution of the workflow, a pharmacophore screening database with 167 records has been obtained. Ethanol was not included in the output database since its pharmacophores do not comprise at least three features with distinct spatial positions (a prerequisite for LigandScouts pharmacophore alignment algorithm). In the generated database, every entry corresponds to a particular drug of the dataset and comprises all of its stereoisomers and associated conformers. This means, a drug will be considered as a hit in the screening process if at least one of its stereoisomers matches the query pharmacophore. This mimics the situation in physiological systems when a drug gets administered as a mixture of stereoisomers and not all isomers of the drug are active towards the target receptor.

Pharmacophores of ligands bound to the assessed target binding sites were generated using the software LigandScout (http://www.inteligand.com/ligandscout3/, version 4.4.8) (Wolber and Langer, 2005) by employing the following general procedure: First, the PDB file providing the structure(s) of the receptor bound ligand was loaded into the “structure-based” perspective of LigandScout. Afterwards, for every binding site of interest (Supplementary Table S2), a structure-based pharmacophore of the ligand was generated (using default settings) by zooming into the corresponding site and initiating the pharmacophore generation process. The obtained structure-based pharmacophores were continuously collected in LigandScouts “alignment” perspective and, after all PDB files have been processed, stored in a single file for further processing.

Pharmacophores representing combinations of multiple structure-based pharmacophores were generated in the “alignment” perspective of LigandScout using the “merged” pharmacophore generation functionality. The pharmacophores “[3/Ach (α7+/α7−)_18 Hits]” (merged pharmacophores of PNU-120596 in PDB entry 7EKT), “[3/Gly (α3+/α3−)_119 Hits]” (merged pharmacophores of ivermectin in PDB entry 5VDH), “[3/GABA (α1+/β2−); (γ2+/β2−)_36 Hits]” (merged pharmacophores of phenobarbital in PDB entry 6X3W) and “[4/(GABA (α1+/α1−)_74 Hits” (merged pharmacophores of pregnanolone in PDB entry 5O8F, tetrahydrodeoxycorticosterone in PDB entry 5OSB and alphaxolone in PDB entry 6CDU) (Supplementary Table S2) were all obtained via a feature-based alignment strategy since the ligands in each of these complexes displayed identical binding modes. However, slight differences were observed for the structure-based pharmacophores of ivermectin in 5VDH. Two of the ivermectin features (14 in total) that were not present in all sites of 5VDH have therefore been marked as optional for matching in the performed pharmacophore screening runs.

For the generation of the pharmacophore models “[3/GABA (γ2+/β2−)_103 Hits]” (merged pharmacophores of diazepam in 6X3X [E] and phenobarbital in 6X3W [E]), “[3/GABA (β(2/3)+/α1−)_84 Hits_]” (merged pharmacophores of propofol in 6X3T [D], etomidate in 6X3V [A] and diazepam in 6HUP [E]) a reference point (Cɑs of proximal binding site residues) alignment strategy has been chosen. This was appropriate since the ligands do not share a common binding mode but nevertheless bind to overlapping regions of homologous receptor sites.

The parallel screening of the generated structure-based and merged pharmacophores against the prepared drug database was performed by means of KNIME workflows using the activity profiling node provided by the LigandScout KNIME extension (version 1.8.0). The corresponding inputs of the activity profiling node were connected to a PMZ-reader node providing the screened pharmacophore dataset and a LDB-reader node specifying the drug database LDB file. For screening with the structure-based pharmacophore dataset, the retrieval mode in the activity profiling node configuration has been set to “Get best matching conformation,” exclusion volume checks were enabled and the omitted feature count was set to zero. For screening with the merged pharmacophore dataset, the retrieval mode was set set to “Stop after first matching conformation,” exclusion volume checks were enabled and the omitted feature count has been set to 9 in order to require only 3 matching drug features for any of the query pharmacophores (even in the case of the largest pharmacophore “[3/Gly (α3+/α3−)_119 Hits]” with 12 mandatory features). After execution of the screening workflows, the obtained results were saved in CSV format via the heatmap view of the activity profiling node.

Based on the data from GnomAD versions 3.1.1 and 2.1.1 and the canonical sequences, an alignment was generated: The alignment was produced using MOE (Inc, C.C.G, 2016) by first aligning AlphaFold-predicted structures (Tunyasuvunakool et al., 2021) for all the main isoforms. For the section of the sequences before strand 1 the published structures were used for reference. The intracellular part of the proteins between TM3 and TM4 was separately aligned for minimization of gaps. After finishing the alignment of the major isoforms, other isoforms were added to the alignment with the exception of the ones too short and diverging for a confident mapping to be possible. Alignment was subsequently used in the Python environment through Biopython (Cock et al., 2009). The alignment is available as Supplementary Item S3.

Sequence alignments have been produced with Clustal X (Larkin et al., 2007) and Promals3D (Pei et al., 2008). Sequence-to-structure alignments were produced with Promals3D, PDBeFold (Krissinel and Henrick, 2004) and MOE. Specific alignments and supporting information are provided as individual FASTA files (Supplementary Item S3) and Supplementary Figures S3, S10. Pairwise sequence similarity was calculated by using the sequence similarity monitor script from MOE (ULC, 2022) with PAM180 substitution matrix.

Binding sites have been identified through the 3D structures of the pentameric ligand-gated ion channels that are found in the PDB. These structures with ligands in the binding sites of interest were further analyzed in MOE. The ligplot function and visual inspection of the binding sites were used to identify the pocket-forming amino acids for subsequent local sequence similarity calculations and the descriptor analysis. For the TMD, the chosen amino acids (colored in the sequence alignment in the Supplementary Item S4) were used for all subsequent procedures (similarity and descriptor calculation, mapping of amino acids and variants to a representative ribbon structure, computation of Cɑ distances for the conformational analysis).

For the proteochemometric (PCM) modeling to model the ligand-target interaction space the z-Scales(3) (Sandberg et al., 1998) and the MSWIM (Zaliani and Gancia, 1999) descriptors have been calculated using the R package “peptides” (Osorio et al., 2015) for each binding site in RStudio 1.4.1717 (http://www.rstudio.com).

The z-Scale(3) and MSWHIM have been plotted in a 3D scatter plot and all pairwise euclidean distances have been calculated by using the Plotly (Plotly Technologies Inc., 2015) and Numpy (Harris et al., 2020) libraries in Python 3.8. Through these distances, a dendrogram and a heatmap have been created with hierarchical clustering by using Seaborn (Waskom et al., 2021). The clusters have been visualized by using Plotly.

All the following visualizations were produced using Matplotlib (Hunter, 2007).

For each position in the reference sequence the sum of minor allele count per 100,000 was calculated. The values were transferred to 7EKT [A] (Zhao et al., 2021). The conserved regions of the proteins were also mapped onto 7EKT [A] and for those positions the values were transformed into RGB values. The values for the rest of the sequence were set to yellow (0xb3b300). The colors were then transferred to MOE and applied to the 7EKT [A] structure.

Variable regions were mapped on the reference alignment. For each gene and each region, we calculated the sum of minor allele count per 100,000, divided by the length of the region in the gene’s major isoform.

The amino acids involved in the formation of binding sites were first assigned to the GABRA1 sequence used in AlphaFold predictions. Then they were transferred to the reference numbering. For each of the binding sites of interest (and the pore facing amino acids on TM2) the per-amino acid number of variants and the sum of minor allele count per 100,000 were calculated and separately visualized.

The lists of pocket forming amino acids for each binding site were used to extract the Cɑ coordinates for each individual amino acid in the original PDB file using a Python script. To identify the different local conformations of the protein binding sites, distances and angles formed by the selected Cɑ were calculated using the Python library Biopython (https://biopython.org/). For the site present on a single subunit (5OSC), single subunits were used. For the sites at subunit interfaces, all dimers were used. All datasets of distances and angles (per subunit, or per dimer) have been processed using the principal component analysis (PCA) to reduce the dimensionality to three dimensions. The PCA was performed by using the Python library scikit-learn (Pedregosa et al., 2012).

All structures that are listed in Supplementary Item S1 were used in the analyses, except those that have missing segments or amino acids (4X5T, 6QFA, 4BOT, 4BOR, 4BOI, 4BON, 4BOO, 4AQ5) or are partly unfolded or misfolded (6D6T, 6D6U).

Homology modeling of GABA_AR α6 was performed by the use of MODELLER 9.21 (Webb and Sali, 2016). The alignment (which is free of INDELS) was done with Promals3D with 6HUO as the template. All the dockings were performed using GOLD (Jones et al., 1997) from Hermes 2020.2.0. The ligand was set flexible by allowing to flip ring corners and detecting internal H bonds. 300 poses were generated per docking run with the option to generate diverse solutions turned on. The CHEMPLP was set as the main scoring function and chemscore for the rescoring. The best 20 docking poses from both scoring functions were used for further analysis.

site 2 (GlyR α3+/α3−): 5TIO: F13, R27, I28, R29, F32, Y78, L83, D84, L85, D86

The ECD features a modified immunoglobulin fold (Brejc et al., 2001), consisting of several highly conserved segments with rather large variable regions at the N-termini, and interspersed into the domain in multiple places (Supplementary Figures S3, S4). The canonical agonist site is localized at ECD interfaces of specific subunits in all families. Allosteric sites are formed by other specific subunit interfaces, such as the high affinity Bz-binding site of GABA_ARs (Figure 1). For the agonist sites of all families, as well as for the Bz-site, high affinity radioligands exist for efficient in vitro screenings. The binding site forming segments comprise conserved and variable parts (Supplementary Figures S3, S4), and the variable parts contain fairly large INDELS. Existing structures provide a glimpse into the structural and conformational variability of the canonical interface pockets, which renders them very challenging targets for rapid in silico methods due to the ambiguities of the INDEL positioning and the appropriate choice of ligand bound conformations (Puthenkalam et al., 2016). While there is considerable interest in allosteric sites for which radioligands are lacking (Ramerstorfer et al., 2011; Puthenkalam et al., 2016), the need for more elaborate in silico methods for suitable INDEL placement go beyond the scope of this study.

For a subset of GABA_AR receptor subunits which align well, an explorative computation of pocket descriptors was performed. For this, we used existing Bz bound structures and applied the pocket descriptor calculation in order to test and compare [z-Scales(3) and MSWHIM]. As the arrangements of α-β-γ containing receptors are generally assumed to be identical with the known situation for α1-β-γ2 as depicted in Figure 1, and α-β receptors are assumed to have a β subunit in the place of the γ subunit, we constructed all hypothetical α+/β−, β+/α−, α+/γ−, γ+/β− and β+/β− interfaces. Based on the available bound states we extracted the binding site forming sidechains and examined the resulting z-Scales(3) and MSWHIM descriptors (Supplementary Figures S5–S7). Both descriptors separate the GABA binding β+/α− interfaces from the Bz-binding α+/γ− interfaces, thus recapitulating the known pharmacology (Supplementary Figure S5). In line with experimental findings, the modulatory α+/β− interfaces share more properties with the Bz-sites than with the GABA sites (Supplementary Figure S6).

In addition to the canonical interface pocket, a novel pocket has been observed in a crystal structure of a GlyR α3 homopentamer (Figures 1, 3) (Huang et al., 2017). Structural superposition with other superfamily members and sequence to structure alignments indicate that this site comprises conserved and variable segments and considerable structural diversity (Figures 3B,C). Sequence similarity analysis for the conserved parts of this site suggest that homologous sites likely exist at interfaces of other subunits across the entire superfamily. The principal component of GABA_AR β, δ and θ subunits is most similar, the much more variable complementary component of GABA_AR β, δ, θ and γ subunits are the strongest candidates for sites that might be used by ligands common to multiple family members. However, the overall size and shape is strongly affected by the variable region 3. A pharmacophore screen into this site provides first insights into the chemical space of candidate ligands. Of the small library we used here (Supplementary Item S2; Supplementary Figure S11), 46 drugs with a wide range of chemical scaffolds generated hits in this site (Supplementary Figures S12, S13). This result is suggestive of a possible contribution from this novel binding site to the polypharmacology within the pLGICs.

Generally, the large number of variable segments present at and near ECD binding site forming regions (Supplementary Figures S3, S4) indicates that these binding sites are of concern only for limited polypharmacology (e.g., the interactions of Bz-site binders with the homologous GABA_A receptor α+/β− interfaces (Iorio et al., 2020), or for the interactions of e.g., taurine with glycine- and GABA_A receptors. The mapping of variants to a representative structure is also not possible for large parts of the ECD and needs to be considered per- subunit. An overview of the variable segments is provided in Supplementary Figure S8.

The most conserved domain is the TMD, where all anion channels share >75% sequence similarity. Among nAChR subunits, conservation is as high, but 5-HT₃Rs and the ZAC- protein feature remarkable variability within the 5-HT₃ family and also relative to nAChR subunits (Supplementary Figure S9). The similarity index for the whole domain is not informative about the individual binding sites, it is only suggesting that the TMD is likely to contain binding sites conserved across many superfamily members. This is also evidenced by the observation that the channel blocker picrotoxin binds to most superfamily members with comparable affinity (Erkkila et al., 2004; Das and Dillon, 2003; Das and Dillon, 2005; Yang et al., 2007; Thompson et al., 2010). The channel blocker site is the only site in the TMD for which radioligands exist, thus, there is a big need to better characterize the other allosteric sites present in this domain and to develop reliable in silico screening methods for these.

The best characterized site in the TMD is localized at the upper portion of the subunit interfaces, where 9 ligand bound structures exist which cover GABA_ARs, nAChRs and GlyRs (Figure 4; Supplementary Item S1). The TMD is nearly free of INDELS, and structurally very conserved. Thus, it is ideally suited to compare pockets on the basis of local sequence similarity, and with the descriptor methods that have been demonstrated as computationally very inexpensive tools to estimate overlapping chemical spaces across pockets. The two descriptors used here (z- Scale(3) and MSWHIM) broadly agree on the similarities within and across family members for both the principal and the complementary components of this site (Figure 4B). The principal component contributing subunits fall into two clusters, one containing all GABA_AR subunits and individual members from all other families. The smaller cluster comprises the β subunits and the majority of the α subunits from the nAChR family, alongside with the GlyR β subunit. The complementary component segments cluster differently, with almost all subunits sharing high similarity for the MSWHIM descriptor with only the ZAC protein and the GLRB protein forming a small separate cluster (Figure 4B; Supplementary Figure S14B). This suggests that the upper TMD interface is likely to have very broadly overlapping ligands for many superfamily members. This analysis was limited to the amino acids that contribute to the binding sites for small molecules as determined from 7EKT, 6HUO, 6X3V, 6X3X, 6X3T and 6X3W. In order to gauge the performance for more detailed predictions, we re-analyzed the amino acid lists for the upper TMD pockets of specific GABA_ARs, as was done for the ECD interface pockets. The resulting clusters follow closely the known distinctions between barbiturate- preferring and etomidate- selective pockets, see Supplementary Figure S14A: Both descriptors generate a homogenous α+/β− cluster, reflecting one of the barbiturate preferring groups. The other barbiturate preferring interface, localized at γ+/β− interfaces, forms clusters containing several β+/α− interfaces for the z-Scale(3) and MSWHIM, suggesting that barbiturate- and etomidate- preferring pockets might have common ligands. This is in fact observed for the low affinity diazepam sites (Kim et al., 2020). The etomidate insensitive β1+/α− form a separate cluster, again recapitulating the known pharmacology.

The ivermectin binding site, also localized at the upper TMD interface, is larger and overlaps partly with the sites for the small molecule ligands discussed above, see Figure 4. In our pharmacophore screen, we observed that a stringent screen into ivermectin bound state, as expected, cannot retrieve small molecule binders, but a less stringent screen will produce many weak hits and is prone to false negatives. Applying the descriptor analysis to the amino acids that form ivermectin contacts results in partly controversial findings between the two used methods, see Supplementary Figure S14C: For example, the principal component of ZACN clusters with GABA_AR subunits in the z-Scale(3) and MSWHIM dataset. It would be interesting to have experimental data with which the results could be validated, requiring more subtypes of pLGICs to be tested. At the present time, several high and low affinity interactions are known suggesting a smaller group of high affinity sites, and a very broadly promiscuoius profile in the micromolar range (Lynagh and Lynch, 2010).

The high similarity for the pockets not only within, but also across families suggests common ligands. Most of the experimental structures that can be employed for structure-based pharmacophore screens are ligand-bound GABA_ARs. Since homologous upper TMD binding pockets exist within receptor families (low affinity diazepam/etomidate site and both barbiturate sites) as well as in different pLGIC family members (7EKT for nAChR and 5VDH for glycine receptors), multiple pharmacophore screens to this binding site were performed (Supplementary Figures S11, S12). Most pharmacophores derived from ligand bound GABA_ARs predict a wide range of drugs to be potentially accommodated by upper TMD interface sites (Supplementary Figure S12), in agreement with the promiscuity of this pocket according to literature (Iorio et al., 2020). All of the 18 substances that produced drug-target hits from pharmacophore screens derived from the homologous site in the nAChR α7 receptor (PDB ID: 7EKT) also generated a hit in at least one homologous site of GABA_ARs or GlyRs. As an example, the plant compound dronabinol (Δ⁹-THC), was predicted to potentially bind at the upper TMD sites of nAChRs, GABA_ARs and GlyRs (Supplementary Figures S11, S12).

Conformational analysis indicates, in line with previous work (Puthenkalam et al., 2016), substantial flexibility of this region. Structures of all family members share overlapping conformations, irrespective of the occupancy of the upper TMD binding site (Figure 5A), but for the cation channel branch of the superfamily, additional conformations have been observed. In these, the RMSD of the pocket backbone differs by up to 2.8 Å from the more frequently observed conformations (Supplementary Figure S17A). The majority of structures populate a conformation space shared by receptors with ligands present at one or several upper TMD interfaces and receptors without ligands in this pocket (apo- site 3, Figures 5B,C). This indicates that ligands of site 3 do not lead to induced fit conformations. Intriguingly, the picrotoxin bound 6UD3 and an ivermectin bound 6VM2 GlyRα1 structure display nearly identical conformation not only of site 3 in the apo- state versus ivermectin- bound state, but of the entire TMD. In contrast, a very different conformation of the upper TMD interface is featured in the GlyR α3 5CFB in complex with the orthosteric antagonist strychnine, see Figures 5A,C. Following on this observation we selected further structures from the region around 5CFB and noted many instances of orthosteric antagonist bound structures in this region of the conformation space, including the bicuculline bound 6HUK β3+/α1− interface, Figures 5A,C and Supplementary Figure S17A. For nAChRs, we observe the same dominant effect of orthosteric antagonists on the TMD pocket conformation, see Supplementary Figure S17A. The strongest effect appears to be a contraction of the upper TMD pocket, mainly on the “long” diagonal between complementary TM2 and principal TM3, see Figure 5C. These observations suggest that orthosteric ligands are dominant drivers of the conformation in the upper TMD, which is further weakly influenced by other factors. Of further interest is the high degree of asymmetry relative to a pseudo symmetrical pentagon that occurs in many structures, such as the 6X3X GABA_AR with diazepam present in the upper TMD, Figure 5C. Interestingly, we find the conformations of chimeric constructs with ECDs from bacterial homologues to occupy distinct regions of the conformation space (Figure 5A; Supplementary Figure S17A).

The lower TMD interface also contains a binding site, which has been characterized as a modulatory site for neuroactive steroids (Figure 6) (Laverty et al., 2017; Miller et al., 2017; Chen et al., 2018). There is big interest in the development of neurosteroid-mimetics as putative drug candidates (Blanco et al., 2018), which also raises the question how prone to cross-reactivity such ligands might be. The analysis we performed based on local sequence similarity and pocket- descriptors indicates that the 45 subunits fall into two clusters, which are similar for the principal and the complementary component (Figure 5C; Supplementary Figure S15). In GABA_ARs receptors, α subunits contribute a glutamine sidechain from TM1, which drives the high efficacy modulation and is absent from all other subunits. The remainder of the site is quite conserved for GABA_ARs and GlyR α subunits, and a few members from the cation permeable family members fall into this cluster. Thus, the unique high efficacy determinant present exclusively on GABA_AR α subunits might render this site sufficiently unique to be of less concern for toxicological alerts.

An additional steroid binding site has been described, which is localized in a position between the lower end of TM3 and TM4 of a GABA_AR α1− subunit and the lipid collar, the TM3/TM4 lipid associated site (Miller et al., 2017, see Figure 7). Interestingly, in another structure, a PIP2 molecule was found in an overlapping localization (Laverty et al., 2019). This region has been subjected to mutational studies aimed at the identification of binding sites for endo- and phytocannabinoids (Sigel et al., 2011; Xiong et al., 2011; Bakas et al., 2017), and is thus considered to be an interaction site for both endogenous and exogenous lipophilic ligands.

The descriptor analysis does not give a consensus picture, as the z-Scales(3) predicts two clusters (one comprising all nAChR subunits, the A and C subunits of the 5-HT₃R and the ZAC protein), while MSWHIM predicts high similarity across the entire superfamily except for ACHD, ACHE and ACHG (Figure 7C; Supplementary Figure S16).

Conformation analysis of the lower TMD reveals that the majority of structures of all superfamily members cluster very tightly for both sites, indicative of very similar conformations irrespective of the presence or absence of ligands in this site, see Figure 8. Additionally, a few structures of nAChRs and 5-HT₃ARs adopt a very different conformation with a local RMSD up to 2.6 Å (Supplementary Figure S17C). In the case of the lipid associated site, this reflects very different conformations of the TM4 helices in the different experimental structures, resulting in local RMSD differences > 3 Å (Supplementary Figure S17D). For this site, the 5-HT₃R structures are uniquely different from all others, see Figure 8B and Supplementary Figure S17D. Surprisingly, the lower TMD ligands seem to not induce any induced fit at all, as the structure with the inhibitory steroid preganolone sulfate in site 5 (5OSC) and all structures with modulatory steroids in sites 4 (5OSB, 6CDU, and 5O8F) superpose very tightly, see Figure 8C.

An unexpected finding in our pharmacophore screens is the extensive overlap of hits for the novel ECD site 2 and the upper TMD interface sites. This raises the question whether it reflects spurious hits, or hints at a genuine overlap of ligand space shared by non-homologous pockets. A similar situation is well established for GABA_A receptors: Diazepam binds with high affinity at the canonical ECD α+/γ− interfaces of certain receptor subtypes, and additionally with lower affinity at two of the upper TMD sites (Masiulis et al., 2019; Kim et al., 2020).

In this study, we did not include the canonical ECD interface sites, as the large number of distinct interfaces and the multiple segments with INDELS would have exceeded the scope of this paper. Our results, however, suggest that the multi-site usage of ligands as displayed by diazepam is not limited to the sites 1 and 3, but may also affect other combinations of binding sites. For the sites localized in the TMD, structural data is yet not very suggestive of ligand promiscuity between the different sites (3, 4, and 5). In contrast, photoaffinity studies strongly suggest that steroid derivatives have multiple binding sites not limited to the lower TMD, but also in the upper TMD (Chen et al., 2019; Sugasawa et al., 2019). It would be very interesting to find out whether endogenous ligands also interact with multiple sites in the ECD and TMD.

While our limited pharmacophore screen does not suggest steroids as hits for the upper TMD sites, we note that 17 compounds, representing considerable chemical diversity, are predicted with moderate to strong scores to match with the pharmacophore models from all three types of sites. This is the case for dronabinol (Δ⁹-THC), another example is the 5-HT₃R interacting antiemetic metoclopramide (Supplementary Figures S11, S12). The latter thus is potentially another example for a compound which interacts with the canonical ECD interface site and additional sites in the TMD, akin to the Bzs. In the screens performed here, we also note that a fraction of the screened Bzs, as well as the Bz-site ligand alpidem, score with 0.7 or higher for all three TMD sites in at least one family member (Supplementary Figures S11, S12).

For phytocannabinoids, the literature demonstrates interactions with multiple family members (Barann et al., 2002; Oz et al., 2004; Yang et al., 2008; Yang et al., 2010; Xiong et al., 2011; Bakas et al., 2017; Schmiedhofer, 2017). Effects differ depending on the protein subtype, as well as the cannabinoid molecule. For example, Δ⁹-THC is relatively ineffective at inhibiting α7 nAChRs, while inhibiting 5-HT₃Rs with similar potency to CB₁ receptors (Barann et al., 2002; Oz et al., 2004). On the other hand, CBD and Δ⁹-THC are potentiating GlyRs and most GABA_AR subtypes (Xiong et al., 2011; Schmiedhofer, 2017). One candidate binding site has been delineated by mutational studies: In GlyR α1 and α3 subunits, a site which was proposed to correspond with the later described pregnanolone sulfate site on GABA_A R α1 subunits in the 5OSC structure, was shown to mediate a large fraction of a strong PAM effect elicited by CBD (Xiong et al., 2011). For 5-HT₃Rs an allosteric modulatory binding site has been suggested (Barann et al., 2002).

Our PH4 screen indicates matches for Δ⁹-THC with the TM3/TM4 lipid associated site 5, as well as with most upper TMD interface sites (3), and with the steroid bound lower TMD-interface (site 4) (Supplementary Figure S12). In addition, another match occurred for the novel upper ECD interface site. We thus performed an exploratory computational docking screen into all matched sites. In turn, we ruled out the upper ECD interface as it displays multiple clashes, and the docking results cannot recapitulate the overlap of the pharmacophore features at all. This leaves the binding sites in the TMD as candidates. We again turned to the literature to search for evidence on binding sites, and to prioritize docking for family members with known effects. For GABA_ARs, both CBD and Δ⁹-THC have been investigated in multiple subunit combinations (Bakas et al., 2017; Schmiedhofer, 2017). Both phytocannabinoids elicit similar effects in most GABA_AR subunit combinations that were tested, enhancing GABA currents in most of them and reducing currents in β1-containing assemblies (Schmiedhofer, 2017). A mutation in TM3 of the GlyR α3 subunit (Xiong et al., 2011) is localized in a position shared by sites 4 and 5, and in close proximity to site 3 (Figure 9).

In order to select subunits for further computational docking, we examined the predicted similarity between GlyR α3 and other candidate subunits. The descriptor analysis of the candidate sites is ambiguous for the upper TMD interface, where the two descriptors place GlyR α3 into different clusters. The placement by MSWHIM with nAChR subunits (Supplementary Figure S14b) would be consistent with the strong hit in the 7EKT structure’s site 3, and was thus followed up in more detail. For the lower TMD interface site, both descriptors cluster the GlyR α3 subunit with GABA_AR subunits (Supplementary Figure S15), in line with the hit for this steroid binding site. Interestingly, for the TM3/TM4 lipid associated site, where MSWHIM suggests broad conservation, the z-Scales(3) places GlyR α3 into a cluster with the GABA_AR α6 subunit (Supplementary Figure S16). This is interesting because the current enhancement by Δ⁹-THC that was observed in a pilot experimental study is stronger for α6-containing assemblies compared to other α- isoforms (Schmiedhofer, 2017).

In order to derive testable structural hypotheses for pockets with combined experimental and computational evidence, we generated docking poses for site 3 in the nAChR α7 homomer, in the steroid bound site 4 of the GABA_AR lower TMD α1+/α1− (original PDBID: 5OSB), and in site 5 in GlyR α3 [to compare with the CBD docking from (Xiong et al., 2011)] as well as in a homology model of GABA_AR α6. The docking results for these four cases were then subjected to two scoring functions to obtain poses with highly ranked consensus scores (Figure 9).

For all four dockings, multiple binding modes with high scores were obtained with both scoring functions. With additional restraints coming from indirect experimental evidence, poses were identified for site 5 which feature GlyR α3 Ser307 in close proximity to the benzo[C]-chromen-1-ol hydroxy group on the ligand, which is in agreement with the pharmacophore screen as well. In addition, this pose has a high overlap with a pose in the GABA_AR α6 model, where both are among the top ten of both scoring functions, see Figure 9. For the docking into the site 3 of the nAChR α7 homomer, the top pose of both scoring functions features high overlap with the PAM ligand of the experimental structure, and thus is in agreement with all available evidence. Lastly, for the results obtained for site 4, the best ranked binding mode overlaps very closely with the ligand from the experimental structure (Figure 9E). An alternative binding mode brings the ligand closer to the TM3 of the principal subunit and the site of mutagenesis (Supplementary Figure S20) (Xiong et al., 2011). These structural hypotheses can be used to guide further mutational studies, see Supplementary Figure S21.