The amino acid sequences that are related to the NR proteins were selected from the NCBI database (ncbi.nlm.nih.gov), using the representative keywords as described in Table I. Protein sequences that responded to the query, but did not include NR members were eliminated from the primary dataset, by using related keywords and regular expression techniques in the header information, and local alignments with reference protein sequences. Furthermore, a final dataset for each member class was produced by using internal protein alignments and protein identity score. Duplicated protein sequences in each phylum that were found to share <90% protein identity within the dataset were removed. In total, 109,613 NR protein sequences were identified from several members and species, and a dataset containing 333 representative, unique, non-duplicate protein sequences of phylum and kingdom was created.

A more specialized analysis was performed in order to identify the protein domains in each NR member. The analysis was performed using the InterPro Database and a representative protein sequence from each NR member class (8,9). In the majority of cases, the Homo sapiens species protein sequence was used as a representative, except in cases where no human sequences existed and thus, another species was selected.

MSA was performed using the MATLAB Bioinformatics Toolbox (https://www.mathworks.com/products/bioinfo.html), utilizing a guide tree and the progressive MSA method as previously described (9,10). Pairwise distances among sequences were estimated based on the pairwise alignment and followed by calculating the differences between each pair of sequences. The Neighbor-Joining method was used towards to estimating the guide tree by assuming equal variance and independence of evolutionary distance estimates (11). A more specific analysis was performed using the MATLAB Bioinformatics toolbox towards comparing the MSA results and the protein domains structural features from the previous step.

Consensus sequence was calculated and visualized through the Jalview platform (12) using the MSA result and parameters including amino acid conservation (9,10,12). The commentary section of Jalview, which presents the amino acid conservation using logos and histograms, was further observed to uncover innovative motifs within the protein domains.

The representative dataset of the steroid hormone members was used for the construction of the phylogenetic analysis. The primary dataset was filtered, and duplicated protein sequences in each species that were found to share <95% protein identity within the dataset were removed. MSA was performed using the MATLAB Bioinformatics Toolbox (13) alignment methods. Only unambiguous homologous regions were retained for phylogenetic analysis; manual masking, trimming and consensus multiple alignments were performed in MATLAB (13). The phylogenetic analysis was performed using the MATLAB bioinformatics toolbox. The distances between the sequences of the SR dataset were measured with the Jukes-Cantor pairwise distance method (14,15) and phylogenetic analysis was performed utilizing the average distance statistical method (16,17) with 100 bootstrap replicates and visualized using MEGA software radiation option (18).

The structural and functional characterization of the LBD domain conserved motifs was analyzed using representative 3D structures of NRs within the steroid hormone subfamily. In the present study, the 3D structures of the glucocorticoid receptor (GR) LBD [4P6X (19) and 5NFT (20)] and the androgen receptor (AR) LBD [6NWL (21)] from the Protein Data Bank (https://www.rcsb.org) were used. The structural and functional analysis at SR structures was gained using MOE (20,22-24), where all the interactions sites with the corresponding proteins and ligands were identified and studied. Specifically, each PDB entry was examined for ligand interaction using the ligand interaction function. MOE showcased the LBD amino acids that interacted with said ligands or co-activators.

The NR proteins are composed of a number of domains and sub-domains which are differentially conserved between the NR members. The NRs protein domains play diverse roles, including the variable N-terminal region, the conserved DBD, the variable hinge region, the conserved LBD and the variable C-terminal region. The representative protein sequences were selected for the analysis of NR domains belonging to the species Homo sapiens, except in a few members from which this species was absent, including FXRβ, Eip78C, Cnr14, EcR, HR96, USP, FTZ-F1β and kni. The FXRβ NR representative was selected from the species Mus musculus, and Eip78C, Cnr14, EcR, HR96, USP, FTZ-F1β and kni from the Drosophila melanogaster species. Sequence analysis of the NR protein domains revealed a clear separation between them (Fig. 2). The characteristic five-protein domains of the NR family were identified in the majority of the members. NR protein domains vary in length among members, and in a few cases, some of them are missing (Fig. 2). Based on these results, some members appear without a C-terminal region, while others do not even have a DBD (SHP, DSS) or LBD (kni) (Fig. 2).

MSA of proteins sequences from the NR superfamily, as well as the subfamily of SRs, was performed to identify highly conservative regions within all organisms of the animal kingdom. The visualization and analysis of the MSA results in both cases were performed using Jalview. As is known, the majority of NRs consist of five structural domains, including the N-terminal, DBD, LBD, C-terminal domains and the hinge region. Sequence analysis of the SR protein sequences dataset has revealed clear conservation in all species. The DBD is the most highly conserved region in the NR superfamily, as observed in the MSA of both NRs and SRs. Two zinc finger motifs can be identified, each of which consisted of a zinc ion linked to four cysteine residues. In the first zinc finger, a conserved pattern, termed P-box, can be identified, containing amino acid residues necessary for the specialized identification and binding of the receptor to the DNA of the target-gene (22). In the second zinc finger, another pattern of amino acids, can be identified, which is involved in the dimerization of the receptors and is termed D-box (23). The amino acids of the D-box appear to be quite diverse in the NR superfamily (Fig. 3), whilst in the SRs, they are highly conserved (Fig. 4) (23).

The LBD appears to be less conserved; nevertheless, some important conserved motifs can be observed, particularly within the SRs. Motif A can be found in positions 944-949 of the steroid receptor MSA (Fig. 5). This is an LLxxL motif (where L represents a leucine residue and x denotes any amino acid), which is an inverse NR box (LxxLL) (24). The NR box (NR interacting box) plays an important role in the interaction of NRs with co-activators and therefore, in the regulation of transcription. Motif B, occupying positions 960-970 of the SR alignment (Fig. 5) and motif C, occupying positions 972-982 (Fig. 5), as reported in the literature, also appear to be regions necessary for co-activator function (25). Motif D, which occupies positions 1051-1063 of the steroid receptor multiple alignment (Fig. 5), appears to be highly conserved and is an area of utmost significance for ligand binding, as mutations in this area have been proven to lead to the complete inability of the receptor to bind (22,24,25). Finally, in alignment positions 1115-1122, a LxxLL motif is observed (motif E), whereas in alignment positions 1144-1153, an inverse NR box (LLxxL) can be found, named motif F (Fig. 5) (24).

As for the NR superfamily (Fig. 6), only motifs B, C, D and E seem to be conserved, whilst in the alignment positions that motifs A and F would be located, a great amino acid variability is observed.

The phylogenetic tree of the steroid hormones comprises of five major monophyletic sub-clusters including the GR, mineralocorticoid receptor (MR), AR, progesterone receptor (PR) and estrogen receptor (ER) (Fig. 7). Based on the tree topology, steroid receptors have been clustered in distinct evolutionary clades. In the phylogenetic analyses, representative members from several kingdoms and phyla were identified, including mammals, Aves (birds), Pisces (fish), insectum (arthropods) and Reptilia.

In the present study, eight candidate conserved motifs that were identified in both DBD and LBD are summarized. Although no structural information is provided for the DBD of the NR, some ‘key’ amino acid residues have been identified and provided within the P-box and B-box motifs, which are mainly termed ‘inc-coordinating motif’ (Figs. 3 and 4) (24,26). The zinc-coordinating motif is a potential pharmacological target of the DBD of the NRs. It is characterized by two anti-parallel helices connected by loops at their amino-terminal ends, from which each helix-loop combine coordinates a single zinc ion using four ‘key’ cysteines (26). Through this specific mechanism, NR members are able to change conformation and activate the selection mechanism towards binding DNA molecules (Fig. 8).

On the other hand, conserved motifs that are formed the interaction site of the LBD of NR are of a major interest as potent pharmacological targets. The majority of these directly interact with ligands, and some others are contained critical patterns, including the ‘LLxxL’ pattern in back and forward direction (Fig. 9). The ‘LLxxL’ pattern is repeated >3 times in the LBD of the majority of the NR members and constitutes a critical conserved region that is directly related to their biological function (24,27). Mutations in the ‘LxxLL’ pattern, have been shown to be strongly associated with incorrect signaling in NR members, particularly in motif D (28). Mitsis et al (24) provided a structural and chemical analysis of the LBD of the NRs, which explains the possible role of this pattern and its significant as a potent pharmacological target.