This study was carried out in accordance with the recommendations of Directive 2010/63/EU of the European Parliament and of the Council of September 22, 2010, on the protection of animals used for scientific purposes. Manipulations of animals and experimental procedures were performed in compliance with the international rules according to the “Guidelines for the care and use of mammals in neuroscience and behavioral research”¹. The research protocol was approved by the Interinstitutional Commission on Bioethics at the ICG SB RAS, 10 Lavrentyev Avenue, Novosibirsk, Russia.

Analysis of the inheritance of agonistic behavior indicators and social dominance levels was conducted on 230 adult male mice that are diallelic crosses of a set of five maternal inbred mouse strains (i.e., PT, DD, YT, A/He, and C57BL/6J) with two analytic inbred paternal strains (BALB/cLac and CBA/Lac) of the murine tendencies in dominance and subordination, respectively, as determined experimentally previously (Bragin et al., 2006).

All the mice were maintained under standard conditions of a conventional animal facility of the ICG SB RAS.

One can see all the 230 diallelic crosses in Table 1, where five rows and two columns present F1 males. In each row of this table, there are descendants of mothers of the same inbred strain. Thus, the maternal non-genetic (pre- and postnatal) and cytoplasmic effects are the same for males of the same row of this table. To exclude non-genetic paternal postnatal effects on offspring, pregnant female mice were isolated from male mice.

We made up groups of F1 hybrid male mice with the minimal society size, namely: two males each: one from each column of the same row of Table 1. In each pair, both male mice had identical age, weight, and body size, but visually differed from each other in color. This approach allowed us to estimate the influence of the paternal genotype on the social dominance level of the appropriate F1 crosses.

A total of 115 experimental pairs (230 F1 hybrids) were distributed into five groups, corresponding to the maternal inbred strains (see Table 1). For each mouse male pair tested, we performed 14 observations (20 min each) during 5 days. Each observation was recorded using a video camera in automatic mode with a fixed period. Next, we analyzed these video recordings using the protocols of software The Observer XT 7.0 (version: 7.0, Noldus Information Technology, license No. OB070-03670). This way, we identified the social rank for each male within the appropriate pair according to asymmetry in agonistic behavior, in particular, by means of attacks and submissive poses as described in the Supplementary Experiment (Supplementary File S6).

Both domesticated and laboratory animals are artificially selected using the known target traits (Belyaev, 1979; Kulikov et al., 2016), which can help in any computer-based genome-wide analysis of these animals (e.g., Zapata et al., 2016) in contrast to the human genome, which is the result of natural selection in favor of unknown unspecializing target traits (Dobzhansky, 1963). Hence, on the basis of our preliminary work (Chadaeva et al., 2017), we proposed the following heuristic decision-making rule: “IF {an excess of a protein relating to the nervous system is an experimentally known physiological marker of low pain sensitivity, fast post-injury recovery, or aggressive, fearless, impulsive, anxious, exploratory, risk/novelty-seeking, anesthetic-like, or similar agonistic-intolerant behavior} AND IF {a given SNP can cause overexpression of a gene encoding this protein} THEN {this SNP can be a SNP marker of predisposition to social dominance} WHILE {the underexpression corresponds to subordination} AND vice versa.” This whole study is devoted to evaluation of this decision-making rule.

Using the aforementioned basic decision-making rule (see subsection “The Basic Decision-Making Rule”), we analyzed all the 5052 SNPs retrieved from the dbSNP database (build 150, Sherry et al., 2001), which are found within the 70 bp promoter regions upstream of the protein-coding transcripts of all the 231 human genes of the neuropeptidergic, non-neuropeptidergic, and neurotrophinergic systems retrieved from database Ensembl (GRCh38/hg38 assembly, Zerbino et al., 2015), which are listed in the alphabetic order in the first columns of Supplementary Tables S1–S3, respectively (hereinafter: see Supplementary Files S1–S3, respectively). These genes encode proteins that are known as key factors altering human social behavior, namely, neurotrophic and growth factors, interleukins, neurotransmitters, receptors, transporters, and enzymes.

Using our public Web service SNP_TATA_Comparator (Ponomarenko et al., 2015), we compared the DNA sequences of the ancestral (wt) and minor (min) alleles of SNPs of the 70 bp promoter region of these genes. We applied it together with the public Web service UCSC Genome Browser (Haeussler et al., 2015) and two public databases dbSNP (Sherry et al., 2001) and ClinVar (Landrum et al., 2014), as described in the Supplementary Web-service (Supplementary File S5). As a result, we obtained two pairs of (–ln(K_D^(wt)) ± δ_(wt)) and (–ln(K_D^(min)) ± δ_(min)) values of TBP affinity for these alleles of the promoter being studied according to contextual, conformational, and physicochemical changes in its B-helical DNA under the influence of a given SNP, as described in the Supplementary Method (Supplementary File S4). Next, we calculated Fisher’s Z-score as follows: Z = abs[ln(K_D^(min)/K_D^(wt))]/[δ²_(min)+δ²_(wt)]^1/2, and in turn found the p-value of statistical significance of this score using package R (Waardenberg et al., 2015).

Finally, using this p-value, we discarded all the SNPs the effects of which were estimated as insignificant; otherwise, using decisions on the SNP-caused significant increase and decrease of the binding affinity of TBP for the analyzed promoters, we predicted the candidate SNP markers for over- or underexpression of the appropriate genes, respectively, as demonstrated experimentally (Mogno et al., 2010). Readers can find all our predictions within the columns “K_D, nM, prediction” of Supplementary Tables S1–S3. Their subcolumns “wt” and “min” contain K_D values of TBP’s binding affinity for the ancestral and minor alleles of the appropriate promoters, respectively. Furthermore, subcolumns “Δ” and “α” correspond to the human gene expression alterations and their statistical significance levels α, which are equal to (1 -p). In addition, subcolumn “ρ” presents a heuristic rank of our predictions varying in alphabetical order from the “best” (A) to the “worst” (E). Finally, Table 2 contains total numbers of our predictions (N_RES) as well as the numbers of the candidate SNP markers for either overexpression (N_>) or underexpression (N_<) of the human genes, as predicted by this work.

For each candidate SNP marker predicted, we manually performed a two-step keyword search in the PubMed database (Lu, 2011) as shown in Figure 1.

As presented in this figure, we handled each candidate SNP marker independently of the others, one by one. First of all, we checked whether the SNP in question was annotated by database ClinVar (Landrum et al., 2014) as depicted in Supplementary Figure S1C (hereinafter: see Supplementary File S5 “Supplementary Web service”) and boldfaced in both the first and third rightmost columns of Supplementary Tables S1–S3.

When this database associated the SNP under study with the human diseases, we manually carried out a primary keyword search for the literature data on the known physiological marker of these diseases, which corresponds to the gene expression alteration predicted for this SNP as described elsewhere (Lu, 2011). Figure 1 depicts this procedure as two boxes consisting of dashed lines. In the case of a successful finding of such a publication, the clinical data taken from database ClinVar (Landrum et al., 2014) indicated the adequacy of our predictions for the SNP under consideration. These confirmations of our predictions are italicized in both the first and third rightmost column of Supplementary Tables S1–S3.

Finally, two dotted boxes in Figure 1 depict our secondary keyword search for the known physiological markers for pain sensitivity, postinjury repair efficiency, or agonistic behavior, which correspond to underexpression of the human gene containing this SNP. This way, we tested the basic decision-making rule of this work (hereinafter: see subsection “The Basic Decision-Making Rule” “Basic decision-making rule”). As the main bioinformatic results, we predicted the candidate SNP markers for the human tendencies in dominance and subordination, which are in both the first and third rightmost column of Supplementary Tables S1–S3. Table 2 contains the total number of these candidate SNP markers (N_↑ and N_↓, respectively).

The section “References” lists the articles cited in Supplementary Tables S1–S3 and in section “Supplementary Method.”

We analyzed dichotomies via the equiprobable binomial distribution and χ² criteria taken from the standard statistical package Statistica (StatSoft^TM, Tulsa, United States).

In the genome-wide study in silico, using only Fisher’s Z-score test, we predicted the candidate SNP markers, the numbers of which for the human gene overexpression and underexpression were compared with one another using the binomial distribution as well as in the case of the human tendencies in dominance and subordination.

During in vivo validation in mice, by means of the χ² criterion, we compared the actual numbers of dominants and subordinates among male mice, which were the F1 hybrids of crossing females from inbred strains of an unknown tendency in social hierarchy with males from two inbred strains BALB∖cLac and CBA∖Lac of the previously experimentally identified tendencies in dominance and subordination, respectively (Bragin et al., 2006).

We applied our experimentally verified public Web service (Ponomarenko et al., 2015) to analyze 395 SNPs in 70 bp proximal promoter regions of 27 human genes encoding neuropeptidergic-system–related proteins, namely: arginine vasopressin receptors (AVPRs), C-X-C motif chemokine receptors (CXCRs); neuropeptide Y and its receptors (NPYs), opioid growth factor receptor (OGFR), opioid receptors (OPRs), oxytocin and its receptor (OXTs), prodynorphin (PDYN), proenkephalin (PENK), prepronociceptin (PNOC), proopiomelanocortin (POMC), and tachykinins together with their precursors and receptors (TACs). The results obtained can be found in Supplementary Table S1.

The human PDYN gene, i.e., the opioid polypeptide hormone prodynorphin, which is a basic building block of endogenous opioid neuropeptides, so-called endorphins, that can inhibit the pain signals peripherally and cause a feeling of euphoria (when acting in the brain) as neurotransmitters of happiness and joy. SNP rs886056538 of this gene’s promoter was annotated within database ClinVar (Landrum et al., 2014), where it is associated with spinocerebellar ataxia as shown in Supplementary Figure S1C. Supplementary Figure S1D illustrates our prediction for this SNP, which is the line “Decision: excess significant” accompanied by the line “Z-score = 2.51, p > 0.95” within the textbox “Result.” This outcome means that this SNP can statistically significantly cause overexpression of this gene. Our primary keyword search (hereinafter: two dashed boxes in Figure 1) produced an original experiment (Smeets et al., 2015) involving a mouse model of the human diseases, which has identified the prodynorphin excess as a physiological marker for spinocerebellar ataxia. As one can see, these in vivo experimental data independently support our prediction for SNP rs886056538 (Supplementary Figure S1). This observation indicates the suitability of our Web service (Ponomarenko et al., 2015) for computer-based analysis of the human genes encoding neuropeptidergic-system–related proteins as italicized in Supplementary Table S1.

After this validation, we manually conducted our secondary keyword search (hereinafter: two dotted boxes in Figure 1) and found the original experiment (Szklarczyk et al., 2012) in a mouse model of human behavior, which associated the prodynorphin excess with reduced conditioned fear. Using our basic decision-making rule within the limitations of the above experimental model of human behavior (Szklarczyk et al., 2012), we predicted that the analyzed SNP rs886056538 can be a candidate SNP marker for the human tendency in dominance (Supplementary Table S1).

Near this clinically characterized SNP marker, we found two unannotated SNPs (rs371345545 and rs557431815), which can also cause overexpression of the human PDYN gene (hereinafter: according to our predictions shown in Supplementary Tables S1–S3). That is why we suggest them as two candidate SNP markers of the same genetic tendencies, namely: spinocerebellar ataxia with limitations (Smeets et al., 2015) and social dominance within the framework of the model (Szklarczyk et al., 2012) as presented in Supplementary Table S1.

This way, we predicted 66 and 31 candidate SNP markers for excess and deficiency of the proteins of the human neuropeptidergic system, respectively, which are also 51 and 46 candidate SNP markers predicted by this work for the human tendencies in dominance and subordination (Table 2 and Supplementary Table S1). First of all, readers can see that the numbers of the candidate SNP markers predicted for the human tendencies in dominance and subordination markers are not statistically significantly different from one another according to equiprobable binomial distribution criterion (P(N_↑ ≡ N_↓ ≡ N_RES/2) > 0.6). This finding is in agreement with our preliminary estimate (Chadaeva et al., 2017), namely: P(N_↑ ≡ N_↓ ≡ N_RES/2) > 0.9.

On the contrary, the numbers of the candidate SNP markers predicted for excess and deficiency of the proteins of the human neuropeptidergic system are significantly different from one another according to the equiprobable binomial distribution criterion (P(N _> ≡ N _< ≡ N_RES/2) < 0.0005) in line with our preliminary observations (Chadaeva et al., 2017), as presented in Table 2: N > = 66, N < = 26 (P(N _> ≡ N _< ≡ N_RES/2) < 0.0005). According to a number of studies, various molecular phenomena can shift frequencies of mutations – e.g., influence of the nucleotide context on the occurrence and repair of pre-mutational damage to genomic DNA, gene conversion, pleiotropic and epistatic effects – Kasowski et al. (2010) first noticed that SNPs decreasing the protein–DNA affinity are much more frequent than SNPs increasing this affinity within the human genome. Next, the authors of ref. (1000 Genomes Project Consortium et al., 2012) quantitatively characterized this mutational shift, namely: there are ∼800 SNPs damaging the transcription factor binding sites and ∼200 SNPs improving these sites per random individual human genome as shown in Table 2. According to Haldane’s dilemma (Haldane, 1957) and neutral evolution theory (Kimura, 1968), this genome-wide estimate can correspond to the neutral mutational drift as a norm. Indeed, we observed 37 clinically proven SNP markers of the human hereditary diseases, which decrease the TBP–promoter affinity, and 14 such SNP markers increasing this affinity (Ponomarenko et al., 2015) in agreement with the above-mentioned genome-wide estimate (Table 2). This pattern matches the commonly accepted opinion on these diseases as a genetic load of the neutral mutational drift in the norm.

Nevertheless, in the case of human reproductive potential, which is considered the target of natural selection, we observed a diametrically opposite pattern, namely: five candidate SNP markers were decreasing the TBP–promoter affinity and 19 candidate SNP markers were increasing this affinity (Chadaeva et al., 2018). Besides, we found (Ponomarenko P. et al., 2017) only a minority (12 of 28) of candidate SNP markers of familial Alzheimer’s disease that can decrease the TBP–promoter affinity; this finding is consistent with natural selection for its very slow pathogenesis, whose clinical manifestation is observed only at the age of over 65 (Table 2). In addition, in the case of core genes of the circadian clock (Ponomarenko et al., 2016), which are naturally selected for continuous coordination between the functioning of systems of the human body and daily fluctuations of the environment, we found 13 candidate SNP markers that can decrease the TBP–promoter affinity and 39 candidate SNP markers increasing this affinity (Table 2).

Looking through Table 2, we noticed that our predictions for the neuropeptidergic gene system are more similar to those for natural selection cases than to those for neutral drift within the normal range. That is why here we predict that the human genes encoding neuropeptidergic-system-related proteins are under natural selection pressure, which equally supports the human tendencies in subordination and domination, as was preliminarily estimated elsewhere (Chadaeva et al., 2017). This way, we followed the semicentennial bioinformatic tradition to compare the actual frequencies of natural mutations within their various dichotomies [e.g., transitions versus transversions (Kimura, 1980) as well as synonymous versus non-synonymous changes (Li et al., 1985)].

Using our public Web service (Ponomarenko et al., 2015), we analyzed 2226 SNPs located within the TBP-binding regions of 109 human genes encoding proteins that are related to the non-neuropeptidergic system, e.g., adenosine receptors (ADORs), adrenoceptors (ADRs), muscarinic cholinergic receptors (CHRMs), nicotinic cholinergic receptors (CHRNs), central cannabinoid receptor 1 (CNR1), catechol-O-methyltransferase (COMT), dopamine D receptors (DRDs), GABA type A receptor-associated proteins (GABARAPs), γ-aminobutyric acid type B receptor subunits (GABBRs), γ-aminobutyric acid - type A receptor subunits (GABRs), G protein–coupled receptors (GRPs), glutamate ionotropic receptor AMPA–type subunits (GRIAs), glutamate ionotropic receptor NMDA-type subunits (GRINs), glutamate metabotropic receptors (GRMs), 5-hydroxytryptamine (serotonin) receptors (HTRs), dopamine transporter DAT (SLC6A3), Na⁺/Cl^--dependent serotonin transporter SERT (SLC6A4), tyrosine hydroxylase (TH), and tryptophan hydroxylase 2 (TPH2). Table 2 and Supplementary Table S2 list the results.

The human COMT gene for catechol-O-methyltransferase has, in its promoter, a clinically annotated SNP, rs777650793, whose association with human cardiovascular disease was documented by database ClinVar (Landrum et al., 2014). Figure 2A presents our prediction for this SNP, which is an excess of this protein. As a non-statistical validation of this prediction, we manually performed our primary keyword search, which resulted in an experimental study (He et al., 2011) on a rat model of human pathologies, which has identified COMT overexpression as a physiological marker of cerebral vasospasm. This correspondence between our prediction (Figure 2A) and these experimental data (He et al., 2011) can support the suitability of the results of our Web service (Ponomarenko et al., 2015) in the case of a study of the human non-neuropeptidergic system as italicized in Supplementary Table S2.

As for our secondary keyword search, it resulted in an in vivo experiment in a rat model of human behavior (Wilhelm et al., 2013), where a catechol-O-methyltransferase excess was a physiological marker of depression. Within the framework of the behavioral animal model (Wilhelm et al., 2013), we predicted the candidate SNP marker of the human tendency in subordination (Supplementary Table S2).

The human DRD3 gene (dopamine receptor D3) carries SNP rs36211802 annotated by database ClinVar (Landrum et al., 2014), which associates it with hereditary essential tremor. This SNP can cause an excess of this receptor, according to our prediction given in Figure 2B. We validated this prediction by our primary keyword search, which found the original experimental data (Kosmowska et al., 2016) on resistance to the high-dose DRD3-agonist treatment of tremor in a laboratory rat model of this human pathology as italicized in Supplementary Table S2.

In addition, our secondary search revealed (Supplementary Table S2) that a DRD3 excess reduced both motor activity and behavioral motivation in a mouse model of human motor activity (Ikeda et al., 2013). This finding allows us to predict rs36211802 as a candidate SNP marker of the human tendency in subordination (Supplementary Table S2).

The human HTR2C gene encodes 5-hydroxytryptamine (serotonin) receptor 2C and carries SNP rs3813929, manifestation of which is an abnormal response to olanzapine (antipsychotic) according to database ClinVar (Landrum et al., 2014). For this SNP, we predict an excess of this serotonin receptor as shown in Figure 2C. Our primary keyword search pointed to the clinical data (Ellingrod et al., 2005) on an HTR2C excess caused by this SNP, whose manifestation is a resistance to olanzapine-caused increase in body mass. It is noteworthy that Tecott et al. (1995) reported that knockout mice (5HT2C^(-)/(-)) are obese, whereas Stahl (1998) observed eating behavior downregulation with a 5HT2C level increase. With this in mind, our prediction of the rs3813929-related 5HT2C excess (Figure 2C) fits the clinical observation of the rs3813929-related resistance to olanzapine-caused increase in body mass (Ellingrod et al., 2005). This agreement between our prediction shown in Figure 2C and the clinical observations (Tecott et al., 1995; Stahl, 1998; Ellingrod et al., 2005) is consistent with our verification of our predictions of this type by electrophoretic mobility shift assays (EMSAs) under equilibrium (Savinkova et al., 2013) and non-equilibrium conditions (Drachkova et al., 2014) in vitro. Besides, this result is in agreement with our verification of our predictions on this subject using biosensor ProteON^TM (Bio-Rad Lab, United States) (Drachkova et al., 2012) and stopped-flow spectrometer SX.20 (Applied Photophysics, United Kingdom) (Arkova et al., 2014, 2017) in real-time mode. In addition, it fits our verification of our analogous predictions using human cell lines transfected with the pGL 4.10 vector (Promega, United States) (for a review, Ponomarenko M. et al., 2017). Finally, it is in line with our verification of our predictions on this subject using independent data from 60 experiments (for a review, see Ponomarenko et al., 2010) and by means of 43 known clinical SNP markers of human diseases (Ponomarenko et al., 2009) and 38 known genetic SNP markers of the breeding traits of animals and plants (Suslov et al., 2010). All these verification data can be a reason for the applicability of our Web-service (Ponomarenko et al., 2015) when the human genes relating to the non-neuropeptidergic system are studied, as italicized in Supplementary Table S2.

Our secondary keyword search yielded empirical data on two laboratory rat strains, which were bred for 60 generations for the presence and absence of high levels of stress-evoked aggression toward humans (Popova et al., 2010). According to these data, increases in both mRNA and protein levels were seen in the brains of non-aggressive rats in comparison with the aggressive ones (Supplementary Table S2). On this basis, we propose the candidate SNP marker for human tendency in subordination (Supplementary Table S2).

In total, we predicted 342 and 163 candidate SNP markers that can increase and decrease, respectively, the expression of the human proteins related to the non-neuropeptidergic system. Besides, these 505 predictions can be clustered as 240 and 265 candidate SNP markers for the human tendencies in dominance and subordination (Table 2 and Supplementary Table S2). As readers can see in Table 2, these results are again consistent with our preliminary estimates (Chadaeva et al., 2017) that natural selection equally supports the human tendencies in dominance and subordination.

We applied our public Web service (Ponomarenko et al., 2015) to study 2431 SNPs in 70 bp regions in front of the TSSs of 95 human genes encoding neurotrophinergic-system–related proteins, namely, adenylate cyclase-activating polypeptide 1 and its receptor (ADCYAP1s), artemin (ARTN), brain-derived neurotrophic factor (BDNF), cerebral dopamine neurotrophic factor (CDNF), ciliary neurotrophic factor (CNTF), fibroblast growth factors and their receptors (FGFs), Fms-related tyrosine kinases and their ligand (FLTs), glial-cell-derived neurotrophic factor (GDNF), GDNF family receptors (GFRs), glia maturation factors (GMFs), insulin like growth factors and their receptors (IGFs), interleukins as well as their receptors and signal transducers (ILs), leukemia-inhibitory factor (IL6-family cytokine) and its receptor (LIFs), nerve growth factor and its receptor (NGFs), neuregulins (NRGs), neuropilins (NRPs), neurturin (NRTN), neurotrophins (NTFs), neurotrophic receptor tyrosine kinases (NTRKs), oncostatin M and its receptor (OSMs), platelet-derived growth factor subunits and receptors (PDGFs), placental growth factor (PGF), persephin (PSPN), Ret receptor tyrosine kinase (RET), transforming growth factors β, its receptors and associated protein 1 (TGFBs), and vascular endothelial growth factors (VEGFs). We show our results in Table 2 and Supplementary Table S3.

The human FGFR2 gene (fibroblast growth factor receptor 2) contains two SNPs rs387906677 and rs886046768, which were clinically detected in patients with bent bone dysplasia syndrome and craniosynostosis, respectively, as documented by database ClinVar (Landrum et al., 2014). Readers can see in Figure 2B, 3A how we predicted the FGFR2 deficiency in the case of rs387906677, whereas rs886046768 corresponds to an FGFR2 excess.

At first, our primary keyword search revealed an experimental report (Merrill et al., 2012) on a mouse model of human embryonic development, which linked bent bone dysplasia with reduced levels of FGFR2. Next, in the same way, we found the original experiment (Mansukhani et al., 2000) on mouse osteoblast cell culture ex vivo that points to FGFR2 as an inducer of apoptosis in these cells and an inhibitor of their differentiation, hyperactivity of which causes craniosynostosis-linked alterations in cell culture. As depicted in the figures, these independent findings confirm the validity of our predictions (Figure 3A,B) in the case of the neurotrophinergic system analysis, as italicized in Supplementary Table S3.

After this validation, our secondary keyword search yielded an article (Meyer et al., 2012) on FGFR2 deficiency as a physiological marker of delayed post-injury skin wound healing. Analogously, we found a biomedical paper (Baatar et al., 2002) on the injections of recombinant human FGFR2 around ulcers, which have accelerated ulcer healing in rats as an animal model of the human pathologies. On this basis, we predicted rs387906677 and rs886046768 as candidate SNP markers of the human tendencies in subordination and dominance, respectively (Supplementary Table S3).

The human PDGFRA gene encodes platelet-derived growth factor receptor α and contains SNP rs183431225 annotated by database ClinVar (Landrum et al., 2014) in connection with both idiopathic hypereosinophilic syndrome and gastrointestinal stromal tumor. Figure 3C presents our prediction for this SNP: overexpression of this receptor. Our primary keyword search revealed two biomedical papers, one of which (Score et al., 2006) reports the PDGFRA excess as a marker of patients with hypereosinophilia, and another one (Hayashi et al., 2015) reveals reduced proliferation of gastrointestinal stromal tumor cells under the influence of a selective inhibitor of PDGFRA. Thus, these independent literature data support applicability of our predictions to the study of human genes encoding neurotrophinergic-system–related proteins as italicized in Supplementary Table S3.

Then, we did our secondary keyword search and found a mouse model of human behavior indicating that the PDGFRA overexpression causes oligodendrocyte-associated nociceptive hypersensitivity to neuropathic pain (Shi et al., 2016). That is why we assumed that rs183431225 is a candidate SNP marker of the human tendency in subordination (Supplementary Table S3).

The human RET gene codes for the Ret proto-oncogene, where two SNPs (rs10900297 and rs10900296) have been associated with three human diseases (renal adysplasia, Hirschsprung disease, and pheochromocytoma) as documented in database ClinVar (Landrum et al., 2014). As readers can see in Figure 3D,E, our predictions for these SNPs surprisingly correspond to over- and underexpression of this gene. Nevertheless, using our primary keyword search, we learned that both an excess (Sarin et al., 2014) and deficit (Bridgewater et al., 2008) of RET are known as physiological markers of renal adysplasia. In addition, both overexpression (Ishii et al., 2013) and underexpression (Zhan et al., 1999) of the RET gene can contribute to the pathogenesis of Hirschsprung disease. Finally, both increased (Huang et al., 2003) and decreased (Moore and Zaahl, 2012) levels of this proto-oncogene are often seen in pheochromocytoma. Thus, the above publications additionally validate our results (Figure 3D,E) as italicized in Supplementary Table S3.

Accordingly, we conducted a secondary keyword search and thus selected two animal models of human behavior. The rat model (Wang et al., 2017) associated the RET excess with hypersensitivity to neuropathic pain. In the mouse model (Golden et al., 2010), the RET deficit reduced epidermal innervation. Within the limitations of these models, we predicted two candidate SNP markers (rs10900297 and rs10900296) of the human tendency in subordination (Supplementary Table S3).

The human TGFBR2 gene (transforming growth factor β receptor 2) contains SNP rs138010137, which occurs in patients with thoracic aortic aneurysm as documented in database ClinVar (Landrum et al., 2014). According to our prediction illustrated in Figure 3F, this SNP can reduce levels of receptor TGFBR2 in humans. Using a primary keyword search, we found an original work about the TGFBR2-deficient aortic aneurysm and aortic dissection as the specific forms of these pathologies (Angelov et al., 2017). As one can see, this is one more argument in favor of the applicability of our Web service (Ponomarenko et al., 2015) to research on the human genes related to the neurotrophinergic system as italicized in Supplementary Table S3.

Next, our secondary keyword search yielded a transgenic mouse model of human health (Martinez-Ferrer et al., 2010), in which the TGFBR2 deficit accelerates healing, closure, and resurfacing of skin wounds. For this reason, we suggest rs138010137 as a candidate SNP marker of the human tendency in dominance (Supplementary Table S3).

Summarizing all the above, we can see 506 candidate SNP markers predicted by this work in the case of human genes encoding the neurotrophinergic-system-related proteins (Table 2 and Supplementary Table S3). These predictions can be grouped into 346 and 160 candidate SNP markers of the excess and deficiency of these proteins, respectively, as well as into 265 and 241 candidate SNP markers of the human tendencies in dominance and subordination (Table 2). Notably, the first of these dichotomies of SNPs in the human genome is statistically significantly uneven, whereas the second one is uniform. This is one more actual piece of evidence for the pressure of natural selection on the human neuron-specific genes, which equally supports the human tendencies in dominance and subordination, in agreement with our preliminary estimates (Chadaeva et al., 2017) as well as with all the other predictions of this work.

Altogether, we analyzed 5052 SNPs within all the TBP-binding regions of all the promoters in front of all the protein-coding transcripts of all the 231 known human neuron-specific genes and selected 1108 candidate SNP markers that can significantly affect the affinity of TBP for these promoters (22%) as shown in the bottom row of Table 2. This result of our exhaustive whole-genome analysis of three systems of the human body (neuropeptidergic, non-neuropeptidergic, and neurotrophinergic) is consistent with both Haldane’s dilemma (Haldane, 1957) and Kimura’s neutral evolution theory (Kimura, 1968). Our in silico fivefold reduction in the number of unannotated SNPs for their subsequent in vivo studies is in line with the current need for reducing the cost of both experimental and clinical searches for valuable SNP markers in the human genome by trial and error through preliminary computer analysis of the known SNPs (Deplancke et al., 2016).

With this in mind, we selected all the 10 among 1,108 candidate SNP markers predicted in this work (Figure 2, 3 and Supplementary Figure S1), which are currently linked to the human diseases by public database ClinVar (Landrum et al., 2014). As described above, we non-statistically validated this set of our selected predictions by our primary keyword search in the public PubMed database (Lu, 2011). Essentially, this match between our 10 selected predictions and the found literature data is statistically significant at the level of α < 0.001 according to the criterion of the equiprobable binomial distribution.

It is important to note that most of the candidate SNP markers that were marked in database ClinVar (Landrum et al., 2014) had a “Clinically insignificant” label because the number of patients with these candidate SNP markers varied from one to six, whereas for clinical significance it is necessary to use cohorts of several hundred patients. This observation supports subsequent verification (using clinical protocols) of the candidate SNP markers predicted by this work. In this way, genotyping for the elite combat athletes in addition to the widely used textual psychological questionnaires for them (Tiric-Campara et al., 2012) could enrich personalized sports medicine.

In addition, we used the semicentennial bioinformatic tradition of comparing the actual frequencies of mutations for their various dichotomies [transitions versus transversions (Kimura, 1980), synonymous versus non-synonymous changes (Li et al., 1985), etc.]. To this end, we grouped all the 1108 predictions into 754 and 354 candidate SNP markers for the increase and decrease in the TBP binding affinity for promoters of the human neuron-related proteins, respectively (Table 2: N_RES, N_> and N_<). This dichotomy contradicts the binomial distribution of the whole-genome ratio 4:1 of the SNPs reducing versus SNPs increasing affinity of the transcription factors for the human gene promoters (1000 Genomes Project Consortium et al., 2012) as neutral drift according to Haldane’s dilemma (Haldane, 1957) and neutral evolution theory (Kimura, 1968), Table 2: p(N _< = 4N _> = 4N_RES/5) < 0.000001. This significant contradiction means the adaptive pressure of natural selection on the human neuron-specific genes is in line with the commonly accepted opinion about the adaptive role of both the nervous system and social behavior in the course of human origin and evolution. That is one more evolutionary argument for the reliability of our predictions made in this work.

Finally, by the same reasoning, we grouped all the 1,108 predictions into 556 and 552 candidate SNP markers for the human tendencies in dominance and subordination, respectively (Table 2: N_RES, N_↑, and N_↓). In contrast to the above dichotomy, this one corresponds to the highly probable H₀ hypothesis about the equiprobable binomial distribution of these candidate SNP markers for human social hierarchy [Table 2: p(H₀: N_↑ = N_↓ = N_RES/2) > 0.9]. This correspondence means that the pressure of natural selection proven above equally supports the human tendencies in dominance and subordination.

Notably, so that natural selection can control the human tendencies in dominance and subordination, it is necessary that this human tendencies can be inherited from generation to generation from parents to offspring. That is why, we in vivo validated our in silico predictions of this work in a mouse model of human inheritance as described below.

Each public Web service addresses a specific sort of regulatory SNP analysis (e.g., Bendl et al., 2016), and each has its specific advantages and disadvantages. Therefore, a comparison between the particular predictions and experimental data as an independent commonly accepted uniform platform (rather than between predictions of various Web services) needs to be a necessary step for prediction of candidate SNP markers in silico (Yoo et al., 2015; Ponomarenko M. et al., 2017). Keeping this in mind, we in vivo validated our in silico predictions on the equal natural-selection support of the human tendencies in dominance and subordination using a mouse model of human inheritance as described in the section “Materials and Methods.” The obtained results are given in Figure 4 and Table 3.

Figure 4 indicates that we completely reproduced the temporal pattern of both formation and maintenance of the social hierarchy in mouse pairs by means of both the number and duration of attacks and submissive poses.

As one can see in the first row “PT” of this table, 21 of 31 mouse males of the F1 hybrids carrying the PT × BALB∖cLac genotype dominated over the male F1 hybrids of the PT × CBA/Lac genotype, and 10 mouse males of the PT × CBA/Lac genotype were dominant in the remaining pairs of the same combination. This actual difference between the F1 male hybrids PT × BALB/cLac and PT × CBA∖Lac is characterized by the χ²-score equal to 3.9, which is statistically significant at the level of α < 0.05. In addition, we observed the same significant dominance of the BALB/cLac-related F1 hybrids over the CBA/Lac-related ones, when the maternal inbred strains were DD and YP (Table 3). In addition, in the cases of maternal inbred strains C57BL/6J and A/He, we found only a tendency for the same dominance, which was insignificant, possibly because of the insufficient number of the appropriate mouse male pairs studied regarding these maternal genotypes. Finally, the last row of Table 3 represents the final result: the statistically significant majority of 79 among 115 BALB/cLac-related male hybrids achieved their dominant social status within their pairs with the CBA∖Lac-related males of the same maternal inbred strains. This finding means that this mouse model of human inheritance reveals an ability of the tendencies in dominance and subordination to be inherited from generation to generation from parents to offspring and, therefore, to be an object of natural selection. This is the main genetic in vivo argument in favor of the reliability of our in silico predictions in this work.

Finally, looking through Figure 4, one can see that, in contrast to the first day of microsocial observation of a pair of adult male mice, which was characterized by numerous and lasting attacks of one mouse on the other, by the end of the second day a social hierarchy is established, with rare short-term ritualized attacks of dominant and/or ritualized submissive poses of a subordinate without any injuries and dangers for their lives and health (Lorenz, 2002). This is the main ecological benefit of establishing and maintaining social hierarchy, as a result of which natural selection equally supports the human tendencies for both dominance and subordination.