Statistical Genetic Approaches to Investigate Genotype-by-Environment Interaction: Review and Novel Extension of Models

Statistical genetic models of genotype-by-environment (G×E) interaction can be divided into two general classes, one on G×E interaction in response to dichotomous environments (e.g., sex, disease-affection status, or presence/absence of an exposure) and the other in response to continuous environments (e.g., physical activity, nutritional measurements, or continuous socioeconomic measures). Here we develop a novel model to jointly account for dichotomous and continuous environments. We develop the model in terms of a joint genotype-by-sex (for the dichotomous environment) and genotype-by-social determinants of health (SDoH; for the continuous environment). Using this model, we show how a depression variable, as measured by the Beck Depression Inventory-II survey instrument, is not only underlain by genetic effects (as has been reported elsewhere) but is also significantly determined by joint G×Sex and G×SDoH interaction effects. This model has numerous applications leading to potentially transformative research on the genetic and environmental determinants underlying complex diseases.


Introduction
There is presently a robust body of approaches to modeling the genotype-by-environment (G×E) interaction in statistical genetics developed from within a linear mixed model framework, which can be divided into genotype-by-dichotomous environment [1][2][3][4][5][6][7][8] and genotype-bycontinuous environment [1,2,7,[9][10][11][12][13] interaction models.The former class has been applied to sex (male/female) [1,[4][5][6][7], disease affection status (affected/unaffected), or presence/absence of environmental exposures (e.g., smoking/non-smoking, basal/high cholesterol-high-fat diet) [2].The latter class has been applied to continuous environments such as the age continuum [7,9], physical activity levels and sedentary behavior [10], and socioeconomic status variables [8,13] (e.g., socioeconomic index, household income, or education levels) among others.It would be advantageous to biomedical investigations to combine the two approaches to address interesting questions like a dichotomous sex-specific response to a continuous index of social determinants of health (SDoH), as we do here, or disease affection-status-specific responses to continuous environmental exposures, for example.Deeper reflection shows a plethora of invigorating, potentially transformative research that could be carried out using such a joint interaction model.Here, we review the G × E interaction models for dichotomous and continuous environments and demonstrate for the first time that they can be integrated into a single unified model.The results of our review show that the joint model can uncover important patterns and phenomena that neither model can elucidate alone.
Since the work of Sir Michael Marmot and colleagues beginning in the seventies, SDoH has come to the forefront of biomedical studies aimed at elucidating the causal mechanisms underlying disease [14][15][16][17].Not surprisingly, the role of SDoH subsequently became an active area of research and application in mental health, particularly depression [18][19][20][21][22][23][24][25][26].In physiological research, stress in general and SDoH-associated stress have been linked to depression, mainly by the hypothalamus-pituitary-adrenal (HPA) axis [27][28][29][30], which is the primary physiological system that deals with stress [31,32].We theorize that SDoH may be conceptualized as a complex environment [33] in which the genetic factors underlying depression respond dynamically.That is, there is an observable genotype-by-environment (G × E) interaction between the genetic basis of depression on the one hand and the (composite multivariate) SDoH environment on the other.The SDoH environment may trigger our genetically based stress response system, which may lead to mental health disorders such as depression or anxiety.Moreover, there is an additional layer of complexity to consider because there is now ample evidence of sex-specific genetic effects underlying depression [34][35][36][37].Thus, it is theoretically possible that the sex and SDoH environments may jointly influence depression outcomes in non-trivial ways.Our novel joint interaction model seeks to investigate how the genetic effects underlying depression may be jointly influenced by the sex and SDoH environments.

Materials and Methods
The University of Texas Rio Grande Valley IRB approved the study protocol.All participants provided informed consent before participating in the study.We evaluated 522 Mexican American participants recruited from the community in an ongoing genetic study for the presence of obesity, diabetes, hypertension, hyperlipidemia, and depression.
The Beck Depression Inventory II (BDI-II) was used to assess the degree of depressive symptoms over two weeks [12,38].The BDI-II is a reliable screening tool for assessing the severity of depression, and it can be administered in both Spanish and English.The BDI-II assesses the severity of depression and is an acceptable screening instrument for depression when administered in both Spanish and English [39][40][41].
For our measure of SDoH, we used the survey of health-related social needs screen developed by Billioux et al. (2017)-known as the Accountable Health Communities Health-Related Social Needs (AHC HRSN) screen-comprising ten questions covering housing instability (2 questions), food insecurity (2 questions), transportation needs (1 question), utility needs (1 question), and interpersonal safety (4 questions) [42].The responses to the ten questions were summed and then divided by 10.
Our main variable of interest in this report is depression as measured by the BDI-II screening instrument.The dichotomous and continuous environments of interest are respectively gender and a derived SDoH environment described in greater detail below.

Polygenic Model
For a generic phenotype vector y, assumed to follow a multivariate normal distribution, we write the following: where X is a matrix of covariates augmented at the left by a column of 1s, β is a vector of the intercept parameter and corresponding regression coefficients, and g and e are unobserved random genetic and environmental effects, respectively [43][44][45].The phenotypic covariance matrix, denoted by Σ, is given as follows: where K and I, respectively, denote the genetic relationship and identity matrices, and σ 2 g and σ 2 e are correspondingly the additive genetic and environmental variance components.In statistical genetics, this base model is called the polygenic or additive genetic model [43][44][45][46].In our initial polygenic models, we accounted for age, sex, age-squared, age-by-sex, and age-squared-by-sex as covariates, and then subjected the residuals to an inverse normalization transformation to induce agreement with the normality assumption [44].The polygenic model is used to obtain estimates of trait heritability, defined as the ratio of the additive genetic variance to the total phenotypic variance, , and as a model reference point upon which more complex models can be elaborated.
Using the polygenic model, the initial SDoH variable (the AHC HRSN score) was found to have a moderate heritability (h 2 = 0.47; p = 4 × 10 −6 ).Given our objective of analyzing SDoH strictly as an environment, we, therefore, employed a best linear unbiased prediction (BLUP) approach to extract the genetic effects that then yielded a predicted value of the environmental effects, which, in turn, can be thought of as a representation of the SDoH variable no longer confounded by genetic effects [45,47].We term this latter BLUP-derived variable the SDoH Index (SDHI) and use it as our SDoH environmental variable.

Modeling the Genotype-by-Environment Interaction for Discrete and Continuous Environments
For a sample of related individuals, assuming fully uncorrelated genetic and environmental effects, the polygenic model posits that the phenotypic covariance is decomposable into additive genetic and residual environmental variance components and that interindividual covariances are determined strictly by the additive genetic variance weighted by the genetic relatedness coefficient (see the polygenic covariance equation).The latter feature of the polygenic model makes two implicit assumptions regarding the genetic covariance: that the pairwise genetic correlation is unity and that the additive genetic variance is homogeneous.Explicitly modeling these assumptions is key to our approach to modeling the G×E interaction.
For the simplest case of contrasting two different environments, the G × E variance is zero if the following two conditions are simultaneously true: homogeneity in the additive genetic variance, σ 2 g1 = σ 2 g2 = σ 2 g , where σ 2 g1 and σ 2 g2 are the additive genetic variances in environments 1 and 2 (for example, male and female status for a G × Sex model or unaffected and affected status for a G × dichotomous disease model), respectively; the genetic correlation (ρ g ) is one across environments, ρ g = 1 [1,6,10,12].Rejection of either or both is evidence that the phenotypic response to the environment has a genetic basis.
We formulate the G × E model for discrete environments in terms of the G × Sex model, with which our group has had great success in elucidating interesting G × Sex or G × dichotomous environment interaction effects (for more detail, see the Appendix A) [4,6].Under the G × Sex model, the total phenotypic covariance can be decomposed into (1) a within-female polygenic model; (2) a within-male polygenic model; and (3) the across-sex additive genetic covariance.The G × Sex model has parameters σ 2 g f , σ 2 gm , σ 2 e f , σ 2 em for the sex-specific additive genetic and environmental variances and ρ G( f ,m) for the across-sex genetic correlation.
We can extend this theory to an environmental spectrum to model G × E for arbitrary continuous environments as opposed to two levels of the environmental variable (for more detail, see the Appendix A).To this end, we employ variance and correlation functions [10,12,13], which we define as follows: , and ρ g = exp −λ g q i − q j , where the additive genetic variance is reparameterized as an exponential function of the value of the environmental variable q for the ith individual, q i , scaled against the sample mean, q, and where the genetic correlation is reparameterized as an exponential decay function of the difference of environmental variables for any pair of individuals i and j, and where α g , γ g , and λ g are parameters to be estimated.The statistical null hypotheses under the reparameterizations for variance homogeneity and genetic correlation stationarity at unity, respectively, are given as γ g = 0 and λ g = 0. To guard against model misspecification bias, we also model the residual environmental variance as a function of the environment in the same way as the additive genetic variance.Thus, the G × SDHI model has parameters α g , γ g , α e , γ e , and λ g .

Joint Genotype-by-Environment Interaction for Discrete and Continuous Environments
We now have the elements to construct a model to allow for a joint accounting of the G × E interaction under both a dichotomous environment and a continuous environment (sex and SDHI in the current case, respectively).The full model (detailed in the Appendix A) is decomposed into five equations specifying the within-female variance ( 1) and ( 2) covariance, the within-male variance (3) and ( 4) covariance, and across-sex covariance functions (all functions of the SDHI environment).The joint model has 11 total parameters, namely α g f , γ g f , λ g f , α e f , γ e f , α gm , γ gm , λ gm , α em , γ em , and ρ G( f ,m) .

Statistical Inferential Theory
For the basic polygenic model, denote the parameter vector by θ = β, σ 2 g , σ 2 e ′ and the residuals vector by ε = (y − Xβ).On assuming multivariate normality of the phenotype vector, the log-likelihood function is as follows: The statistical genetics package SOLAR was used to obtain the model likelihoods, maximum likelihood estimates (MLEs) of model parameters, and their standard errors (SEs) [48].Hypothesis tests were performed by way of the likelihood ratio test (LRT) statistic, which is given as follows: where for the simplest example, θ N denotes the parameter vector for which a single parameter is constrained to 0 and all other parameters are free to be estimated at their MLEs, and θ denotes the fully unconstrained parameter vector.In this case, the LRT is distributed as a chi-square random variable with degrees of freedom (df) given by the difference of the number of constrained and unconstrained parameters, which in this simplest case is 1 df: Λ ∼ χ 2  1 .It is necessary at this point to distinguish between so-called standard and non-standard conditions.Under standard conditions, the null hypothesis is not on a boundary of the acceptable parameter space, in which case the usual asymptotics for the limiting distribution of the LRT hold.For example, a regression coefficient, being essentially a slope term, takes values on the real line and, more to the point, the null hypothesis of β = 0 is not on a boundary.Under non-standard conditions, however, the null hypothesis is on a boundary of the parameter space, in which case the asymptotic limiting distribution for a 1 parameter difference can be shown to be given by a 1  2 : 1 2 mixture of a chi-square random variable with a point mass of 0, denoted by χ 2 0 and χ 2 1 [44,[49][50][51][52][53].Thus, for this non-standard case, we have Λ ∼ 1 2 χ For example, testing h 2 = 0 is such a 1-parameter testing scenario where the LRT follows this last mixture distribution.
Our first step is to establish that there is a genetic basis for the trait of interest, which amounts to testing and rejecting h 2 = 0. Moving forward from this point, we advocate a staged hypothesis testing approach that our team has successfully developed and implemented in analyses of the gene-by-environment interaction [10,12,13].In the first stage, we compared the polygenic model to the G × Sex and to the G × SDHI models.Both interaction models are 5-parameter models.Thus, there is an overall 3-parameter difference between the polygenic model (with parameters σ 2 g and σ 2 e ) and either G × E model.As made clear below in the next stage, 2 of the 3 parameters making up the difference may be considered to be under standard conditions, whereas the remaining parameter can be shown to be under a non-standard condition.Given that, in general, chi-squares are additive, these considerations give rise to Λ ∼ χ 2  1 In the second stage, we examine the specific sources of potential G × E. In particular, we examine the null hypothesis of additive genetic variance homogeneity (σ 2 g f = σ 2 gm and γ g = 0 under G × Sex and G × SDHI, respectively) and a genetic correlation equal to 1 (ρ G( f ,m) = 1 and λ g = 0 under G × Sex and G × SDHI, respectively).We note that these reparameterizations and their resulting hypothesis tests are predicated on the principle that the likelihood function is invariant under one-to-one reparameterization [46,[54][55][56][57][58].Under the null, the G × Sex model reduces to the polygenic model, whereas the G × SDHI model reduces to a reparameterized polygenic model: exp α g + exp(α e ).The additive genetic variance homogeneity null under the G × Sex model is a standard scenario because it is algebraically equivalent to testing the null hypothesis that their difference equals 0, σ 2 g f − σ 2 gm = 0, where this difference may take values on the real line and is thus not on a boundary.Regarding the G × SDHI model, γ g is essentially a slope term on the log-linear scale and, similar to the case for regression coefficients, the null hypothesis of γ g = 0 is thus not on a boundary.Therefore, for either model, Λ ∼ χ 2  1 .As for the null hypothesis of a genetic correlation equal to 1 under the G × Sex model, this is clearly on the right boundary of the acceptable parameter space for any correlation coefficient, which takes values in the closed interval [−1, +1].As for the G × SDHI model, it happens that λ g = 0 is on the left boundary of the permissible parameter space for the exponential decay function, which corresponds to the right boundary of the genetic correlation coefficient because for λ g = 0 we have ρ G = exp −λ g q i − q j = e 0 = 1.Thus, for both cases, we have Under the joint G × Sex and G × SDHI model, we now advocate a third stage, where we test the joint model against whichever of the two 5-parameter models (G × Sex or G × SDHI) has the higher likelihood (this difference need not be significant).The joint model has 11 total parameters, namely α g f , γ g f , λ g f , α e f , γ e f , α gm , γ gm , λ gm , α em , γ em , and ρ G( f ,m) .Therefore, there is a 6-parameter difference, whereupon focusing on the parameters relevant for G × E testing (or environmental variance heterogeneity), we have 2 parameters on a boundary (λ g f and λ gm ) and 4 not on a boundary (γ g f , γ gm , γ e f , and γ em ).The sum of 4 χ 2 1 variables gives χ 2 4 .Further, to jointly determine the mixture distribution for λ g f and λ gm we sum their individual mixture distributions as follows: 1  2 χ . The fourth stage would then consist of the individual tests with limiting distributions given by their corresponding parameters under either the G × Sex or G × SDHI models.Note that although ρ G( f ,m) was not relevant in deriving the LRT mixture distribution for the comparison of the joint model to the G × Sex model (because this is the only parameter present under both models), it is necessary in this last stage to test if it is significant by using the same procedure given earlier for the G × Sex model.
Before concluding this section, we note that two models, namely the G × SDHI and joint interaction models, had a parameter with a standard error (SE) greater than its corresponding maximum likelihood estimate (MLE).Likelihood theory shows that the MLE of a parameter should be greater than twice its SE if it is significant [57].Given this principle, we elected to constrain the parameter in question (where the MLE was less than the SE) to its null value following our previously published G × E investigations where this issue was addressed [10,12,13].To be sure, we also formally tested the parameter by the appropriate test discussed above and confirmed it to be non-significant.In the ensuing sections of this study, such models are referred to as a reduced (abbreviated Red.) version of the model.More importantly, the resulting LRT distributions were accordingly modified by one less parameter following the principles detailed above.

Comparison of Sex-Specific Additive Genetic Variance Functions
In order to enable statistical comparison of the sex-specific additive genetic variance functions at fixed values of the SDHI continuum (minimum, mean, and maximum), we computed their sampling variances using a second-order Taylor expansion approximation for a multivariable non-linear function (see Appendix A) [59].
These sampling variances were then used to compute the adjusted confidence intervals for a two-sided hypothesis test [60].We note that we have adjusted our confidence interval to correspond to a two-sided test of the hypothesis of no difference [61][62][63].It can be shown that the proper confidence level corresponding to a significance level of α ≤ 0.05 is conservatively given as 84%, yielding a confidence interval of 8% and 92% for the lower and upper bounds, respectively.We also performed a modified Wald test [46,54,56,57,64] denoted by W, as follows: where σ 2 g f |q and σ 2 gm|q are the sex-specific additive genetic variances for females and males, respectively, expressed as a function of a fixed SDHI value, denoted by q, and where Var σ 2 g f |q and Var σ 2 gm|q are their respective variance approximations.Given that we are testing if the difference is different from 0, W ∼ χ 2 1 .

Results
In Table 1, we report the demographic characteristics of the study sample by sex.The heritabilities of both BDI-II and AHC HRSN, each estimated under a polygenic model, were both found to be significant (Table 2).In all analyses, the focal variable is BDI-II.Given that we were interested in using the AHC HRSN score as a measure of the SDoH environment, we used BLUP to statistically extract its genetic effects, at which point we termed it the SDoH index (SDHI).We proceeded to separately test the performance of the G × Sex and G × SDoH models of BDI-II against the polygenic model for this variable.From the formal comparison in Table 3, we infer that the G × Sex model provides a significantly better fit to the data than the polygenic model.Thus, in Table 4, we proceeded to test the individual hypotheses of homogeneity in the additive genetic and residual environmental variances, and of the genetic correlation equal to 1.We found that there is significant G × Sex interaction due to the cross-sex genetic correlation being significantly different from 1.There was also evidence of significant residual variance heterogeneity.The LRT statistic for this comparison is distributed as a 50:50 mixture of chi-squares with 2 and 3 degrees of freedom.
Similarly, we found that the G × SDoH model performed significantly better than the polygenic model in explaining the data (Table 5).We thus proceeded to test the slopes of the additive genetic and residual environmental variance functions, and whether the genetic correlation function was different from 1 (Table 6).We found that there was evidence of G × SDoH interaction due to both additive genetic variance heterogeneity and a genetic correlation function that decays away from 1 with increasing differences in the environmental index, SDHI (Figures 1 and 2).The distributions of the LRT statistics are a chi-square with 1 df and a 50:50 mixture of a chi-square with a point mass at 0 and a chi-square with 1 df."---" indicates the parameter constrained to the null hypothesis under the reduced G×E interaction model.The distributions of the LRT statistics are a chi-square with 1 df and a 50:50 mixture of a chisquare with a point mass at 0 and a chi-square with 1 df."---" indicates the parameter constrained to the null hypothesis under the reduced G×E interaction model.Up to this point, we have demonstrated evidence of G × Sex and G × SDoH interactions.However, each interaction analysis was conducted separately from the other, and supposing that our investigation ended here, we could not be certain if both significant interactions hold when considered together.Therefore, we performed a joint interaction model analysis to address whether both interactions are still important when considered together.We selected the G × Sex interaction model to be the baseline model to assess the  Up to this point, we have demonstrated evidence of G × Sex and G × SDoH interactions.However, each interaction analysis was conducted separately from the other, and supposing that our investigation ended here, we could not be certain if both significant interactions hold when considered together.Therefore, we performed a joint interaction model analysis to address whether both interactions are still important when considered together.We selected the G × Sex interaction model to be the baseline model to assess the significance of the joint interaction model because it had a higher likelihood than the G × SDoH interaction model.On comparing the joint interaction model to the G × Sex interaction model, we found that it had a significantly better fit to the data (Table 7).We further investigated the potential significance of the sex-specific variance and correlation functions and the cross-sex genetic correlation (Table 8).We found that there was still evidence of the G × Sex interaction in this case due to the genetic correlation being significantly less than 1.There was evidence of the G × SDoH interaction in males due to additive genetic variance heterogeneity.Red.G × E interaction model −207.858 "---" indicates the parameter constrained to the null hypothesis under the reduced G×E interaction model.
Our final analysis compared the sex-specific additive genetic variance functions at the minimum, mean, and maximum SDHI values (Table 9; Figure 3).The male and female additive genetic variance functions significantly differed at the minimum SDHI value but not at the mean or maximum.Consistent with the finding of male-specific additive genetic variance heterogeneity (Table 8), the adjusted confidence intervals for the male-specific additive genetic variances do not show overlap when comparing the variances at the minimum to the mean SDHI values and when comparing the variances at the mean to the maximum SDHI values, indicating a sustained and significant increase in the male-specific additive genetic variance with increasing SDHI (Table 9; Figure 3).

Conclusions
In this study, we developed a novel model to jointly account for genotype-by-environment interactions for dichotomous and continuous environments.In particular, we applied this joint interaction model to account for the G × Sex and G × SDoH interaction influencing depression.A motivating factor for developing this model is that it allows us to establish if both types of interaction are important independent of one another, similar to the rationale underlying multivariate logistic regression models.We were able to show that there is G × SDoH interaction but only in males and that there is G × Sex interaction due to the cross-sex genetic correlation being significantly different from 1, which indicates that depression is underlain or influenced by different sets of genes in males and females.This model has potentially critical applications in medical research.For example, the dichotomous environment component of the model could be used to model the affected and unaffected states of a given disease of interest, while the continuous environment component could be used to investigate how the genetic response of individuals in both states may be different across a continuous environment of interest, such as physical activity or sedentary behavior.Whatever the case, we are confident that this novel model can lead to fruitful investigations on the genetic basis of response to the environment.

Conclusions
In this study, we developed a novel model to jointly account for genotype-by-environment interactions for dichotomous and continuous environments.In particular, we applied this joint interaction model to account for the G × Sex and G × SDoH interaction influencing depression.A motivating factor for developing this model is that it allows us to establish if both types of interaction are important independent of one another, similar to the rationale underlying multivariate logistic regression models.We were able to show that there is G × SDoH interaction but only in males and that there is G × Sex interaction due to the cross-sex genetic correlation being significantly different from 1, which indicates that depression is underlain or influenced by different sets of genes in males and females.This model has potentially critical applications in medical research.For example, the dichotomous environment component of the model could be used to model the affected and unaffected states of a given disease of interest, while the continuous environment component could be used to investigate how the genetic response of individuals in both states may be different across a continuous environment of interest, such as physical activity or sedentary behavior.Whatever the case, we are confident that this novel model can lead to fruitful investigations on the genetic basis of response to the environment.
where K and I, respectively, give genetic relationship and identity matrices, and σ 2 g and σ 2 e are correspondingly the additive genetic and environmental variance components.In statistical genetics, this base model is called the polygenic or additive genetic model [43][44][45][46].
It is instructive at this point to provide the following scalar phenotypic covariance equation, which specifies the elements of the covariance matrix for all possible ij-comparisons: where σ p y i , y j denotes the phenotypic covariance, k ij is the ij-th element in K, and For a sample of related individuals, assuming fully uncorrelated genetic and environmental effects, the polygenic model posits that the phenotypic covariance is decomposable into additive genetic and residual environmental variance components and that interindividual covariances are given strictly by the additive genetic variance weighted by the genetic relatedness coefficient (see the polygenic covariance equation).The latter feature of the polygenic model makes two implicit assumptions regarding the genetic covariance: that the pairwise genetic correlation is unity and that the additive genetic variance is homogeneous.Explicitly modeling these assumptions is key to our approach to modeling the G × E interaction.
For the simplest case of contrasting two different environments, the G × E variance is zero if the following two conditions are simultaneously true: homogeneity in the additive genetic variance, σ 2 g1 = σ 2 g2 = σ 2 g , where σ 2 g1 and σ 2 g2 are the additive genetic variances in environments 1 and 2 (for example, male and female status for a G × Sex model or unaffected and affected status for a G × dichotomous disease model), respectively; the genetic correlation (ρ g ) is one across environments: ρ g = 1.Denoting the G × E variance as σ 2 g∆ , we have the following expression: There is G × E evidence if either null hypothesis is rejected [1,6,10,12].Rejection of either or both is evidence that the phenotypic response to the environment has a genetic basis.
We formulate the G × E model for discrete environments in terms of the G × Sex model, with which our group has had great success in elucidating interesting G × Sex or G × dichotomous environment interaction effects [4,6].Let there be an indicator variable specifying the sex of individuals in the study, denoted by s i , as follows: where Z is the set of males or, for generic clinical dichotomous variables, individuals in the "1-class" or affected class.The phenotypic covariance can be decomposed under the G × Sex model as follows: and s j = 0 or s i = 0 and s j = 1 Note that the decomposition covers three general classes, which in descending order are the (1) within-female polygenic model; (2) within-male polygenic model; and (3) the across-sex additive genetic covariance.To write the matrix model, let there be a n × 1 indicator vector s with elements s i , as defined above.To construct incidence matrices for pairwise comparisons consisting of both individuals in the male class and both individuals in the female class, we write, respectively, M n×n = s nx1 s ′ 1×n and F n×n = r nx1 r ′ 1×n , where r = 1 − s and 1 is a vector of 1s [4].For the across-class pairwise comparisons, we write C n×n = s nx1 r ′ 1×n + r nx1 s ′ 1×n [4].The covariance matrix for this model may be given as follows: where denotes the Hadamard matrix multiplication operator.We can extend this theory to an environmental spectrum to model G × E for arbitrary continuous environments as opposed to two levels of the environmental variable.To this end, we employ variance and correlation functions [10,12,13], which we define as follows: σ 2 g = exp α g + γ g (q i − q) , and ρ g = exp −λ g q i − q j , where the additive genetic variance is reparameterized as an exponential function of the value of the environmental variable q for the ith individual, q i , scaled against the sample mean, q, and where the genetic correlation is reparameterized as an exponential decay function of the difference of environmental variables for any pair of individuals i and j, and where α g , γ g , and λ g are parameters to be estimated.These functions can be interpreted as the variance and correlation functions of a Gaussian stationary stochastic process [8,13,[67][68][69][70], where the index variable of the stochastic process is the SDHI environment.The statistical null hypotheses under the reparameterizations for variance homogeneity and genetic correlation stationarity at unity, respectively, are given as γ g = 0 and λ g = 0. To guard against model misspecification bias, we also model the residual environmental variance as a function of the environment in the same way as the additive genetic variance.The phenotypic covariance function, with components given by the variance and covariance functions of a Gaussian stationary stochastic process, can be partitioned as follows: e ∀ i = j; δ ij = 1; q i = q j σ i σ j ρ g ∀ i ̸ = j; δ ij = 0; exp α g + γ g (q i − q) + exp[α e + γ e (q i − q)] ∀ i = j; δ ij = 1; q i = q j exp α g + γ g (q i − q) 1 2 exp α g + γ g q j − q 1 2 • • • × exp −λ g q i − q j ∀ i ̸ = j; δ ij = 0; q i ̸ = q j =      exp α g + γ g (q i − q) + exp[α e + γ e (q i − q)] ∀ i = j; δ ij = 1; q i = q j exp α g + 1 2 γ g q i + q j − 2q − λ g q i − q j ∀ i ̸ = j; δ ij = 0; q i ̸ = q j We define a genetic covariance matrix Ψ = ψ ij , with elements given as follows: exp α g + γ g (q i − q) ∀ i = j; δ ij =; q i = q j exp α g + 1 2 γ g q i + q j − 2q − λ g q i − q j ∀ i ̸ = j; δ ij = 0; q i ̸ = q j We posit a diagonal matrix Ω = diag{ω ii } with diagonal elements containing the residual environmental variance function, ω ii = exp[α e + γ e (q i − q)].The covariance matrix for this G × E model for continuous environments is then given as follows:

. Joint Genotype-by-Environment Interaction for Discrete and Continuous Environments
We now have the elements to construct a model to allow for a joint accounting of the G × E interaction under both a dichotomous environment and a continuous environment (sex and SDHI in the current case, respectively).A necessary first step is to posit sex-specific variance and correlation functions subscripted by f for females and m for males.Further,

Figure 1 .
Figure 1.(A) Additive genetic variance function.(B) Genetic correlation function.The do ed line represents the expected value of the genetic correlation under the null.The red en dash line represents the SDHI value at which the genetic correlation has decayed to half the maximum value.

Figure 1 .
Figure 1.(A) Additive genetic variance function.(B) Genetic correlation function.The dotted line represents the expected value of the genetic correlation under the null.The red en dash line represents the SDHI value at which the genetic correlation has decayed to half the maximum value.Genes 2024, 15, x FOR PEER REVIEW 9 of 18

Figure 2 .
Figure 2. The genetic covariance function showing how the additive genetic variance as a function of SDHI and the genetic correlation as a function of SDHI differences jointly influence the overall covariance.

Figure 2 .
Figure 2. The genetic covariance function showing how the additive genetic variance as a function of SDHI and the genetic correlation as a function of SDHI differences jointly influence the overall covariance.

Figure 3 .
Figure3.Sex-specific additive genetic variances as functions of SDH Index (SDHI) at the minimum, mean, and maximum SDHI values.The horizontal dashed black lines demonstrate no overlap between the confidence intervals for the male-specific additive genetic variances at the minimum and the mean and at the mean and the maximum.The different color en dash lines indicate the adjusted confidence intervals corresponding to the sex-specific colors (pink for females and blue for males in this case).

Figure 3 .
Figure 3. Sex-specific additive genetic variances as functions of SDH Index (SDHI) at the minimum, mean, and maximum SDHI values.The horizontal dashed black lines demonstrate no overlap between the confidence intervals for the male-specific additive genetic variances at the minimum and the mean and at the mean and the maximum.The different color en dash lines indicate the adjusted confidence intervals corresponding to the sex-specific colors (pink for females and blue for males in this case).
denotes the Kronecker delta[65,66].The top case on the right-hand side gives the diagonal elements of the covariance matrix under the polygenic model, and the bottom case gives the off-diagonal elements.Appendix A.2. Modeling the Genotype-by-Environment Interaction for Discrete and Continuous Environments

Table 1 .
Demographic characteristics of the sample.

Table 2 .
Heritability analysis of BDI-II and AHC HRSN screening on residualized normalized data.

Table 3 .
Testing the G × Sex interaction model against the polygenic model for BDI-II.

Table 4 .
Testing the critical parameters of the G × Sex interaction model for BDI-II.
G × Sex interaction model −229.230

Table 5 .
Testing the G × SDHI interaction model against the polygenic model for BDI-II.

Table 6 .
Testing the critical parameters of the G × SDHI interaction model for BDI-II.

Table 7 .
Testing the joint G × Sex and G × SDHI interaction model against the G × Sex for BDI-II.

Table 8 .
Testing the critical parameters of the joint G × Sex and G × SDHI interaction model for BDI-II.

Table 9 .
Comparison of sex-specific additive genetic variances at the minimum, mean, and maximum SDHI values (95% confidence interval at the lower and upper bounds).
* See text on the adjusted lower and upper bounds.Genes 2024, 15, x FOR PEER REVIEW 11 of 18