A SNP, Gene, and Polygenic Risk Score Approach of Oxytocin-Vasopressin Genes in Adolescents' Loneliness.

Not much is known regarding underlying biological pathways to adolescents' loneliness. Insight in underlying molecular mechanisms could inform intervention efforts aimed at reducing loneliness. Using latent growth curve modeling, baseline levels and development of loneliness were studied in two longitudinal adolescent samples. Genes (OXTR, OXT, AVPR1A, AVPR1B) were examined using SNP-based, gene-based, and polygenic risk score (PRS) approaches. In both samples, SNP- and gene-based tests showed involvement of the OXTR gene in development of loneliness, though, significance levels did not survive correction for multiple testing. The PRS approach provided no evidence for relations with loneliness. We recommend alternative phenotyping methods, including environmental factors, to consider epigenetic studies, and to examine possible endophenotypes in relation to adolescents' loneliness.

Evolutionary theories suggest that inclusion in social groups is an essential prerequisite for survival, because this brings mutual opportunities for (social) care, protection and assistance (Cacioppo, Cacioppo, & Boomsma, 2014; MacDonald & Leary, 2005). One of the consequences of the subjective experience of lacking those significant social relations is loneliness. Framed otherwise, loneliness is the experienced discrepancy between actual and desired social relations (Perlman & Peplau, 1981). Loneliness is an aversive psychological condition and a known risk factor for several physical and mental health outcomes (Hawkley & Cacioppo, 2010; Holt‐Lunstad, Smith, Baker, Harris, & Stephenson, 2015). The absence of important personal and close relationships could negatively affect well‐being. Despite the overwhelming amount of possibilities to connect with others, both offline and online, loneliness has high prevalence numbers. Numbers suggest a peak of loneliness during adolescence, a decline in middle age, and a slight increase among the elderly (Perlman & Landolt, 1999; Victor & Yang, 2012). Generally speaking, the highest prevalence rates have been found in periods of developmental shifts that are accompanied by changes in one's social experiences (Heinrich & Gullone, 2006; Qualter et al., 2015). Adolescence is eminently a period in which youngsters strive to accomplish more autonomy and individuality. Next to bodily changes and identity formation, relations with peers become increasingly important (Brown & Klute, 2003; Steinberg & Morris, 2001). In this maturation phase with challenging developmental tasks, changes within social environments, and transitions to high school, more feelings of separateness can occur. If their actual needs for social affiliation are not fulfilled in this sensitive period, adolescents are especially vulnerable to experiencing feelings of loneliness (Heinrich & Gullone, 2006). This manifests itself in a peak of loneliness during adolescence, with up to 20% of adolescents reporting feeling lonely “sometimes” or “often,” which even increases to 70% in late adolescence (Qualter et al., 2015). However, studies have shown different developmental trajectories in adolescence (e.g., high increasing levels, or persistent high loneliness; Dulmen & Goossens, 2013).

The aversive feelings associated with loneliness alert people to restore belongingness levels in order to re‐establish social relationships. This is comparable with the aversive feeling of being hungry, telling one to eat (Cacioppo et al., 2014). From an evolutionary perspective, this reconnection tendency should ideally contribute to decreasing levels of loneliness. Being in groups and having significant relationships come with the biological adaptive benefits of increased chances of survival and passing on genes to the next generation (Cacioppo et al., 2006, 2014). Indeed, it has been argued that increased attention to various social stimuli promotes social learning in lonely people and that experiencing isolation leads to behavioral adaptions (e.g., other‐oriented motivations) aimed at increasing social inclusion (Cacioppo et al., 2014; Qualter et al., 2015). This evolutionary advantage of loneliness would predict that loneliness has a substantial degree of heritability (Cacioppo et al., 2014). Indeed, twin studies have shown heritability estimates ranging between 40% and 48% in adolescent samples (Boomsma, Cacioppo, Slagboom, & Posthuma, 2006; Boomsma, Willemsen, Dolan, Hawkley, & Cacioppo, 2005; Goossens et al., 2015; Waaktaar & Torgersen, 2012). Given the high prevalence and the heritability of loneliness among adolescents, this study focuses explicitly on the genetics of loneliness in adolescents.

Genes within the oxytocin‐vasopressin pathway could play a role in loneliness. Genes in this pathway influence levels of oxytocin and arginine vasopressin in the brain. These proteins have been shown to modulate the limbic system both in mammals and humans (e.g., amygdala activation and coupling of amygdala to brainstem regions; Bale, Davis, Auger, Dorsa, & McCarthy, 2001; Domes et al., 2007; Huber, Veinante, & Stoop, 2005; Kirsch et al., 2005), explaining their influence on various social and emotional behaviors related to social inclusion (Insel, 2010), bonding and attachment (Gordon, Zagoory‐Sharon, Leckman, & Feldman, 2010), and behaviors that comprise social skill deficits (Dolen, 2015). The oxytocin‐vasopressin (OT‐AVP) neural pathway is a promising avenue for investigation with regard to loneliness, as one should be able to bond with others to establish or maintain close relationships, and having social skill deficits would definitely impede such bonding (Eronen & Nurmi, 1999; Inderbitzen‐Pisaruk, Clark, & Solano, 1992; Jones, Hobbs, & Hockenbury, 1982; Jones, Sansone, & Helm, 1983; Segrin, 1999; Segrin & Flora, 2000; Spitzberg & Hurt, 1987). This line of reasoning suggests that variations within genes involved in this pathway (i.e., genes encoding the OT‐AVP neuropeptides as well as their receptors) could have an impact on loneliness.

Functional genetic variants that influence regulation of the OT‐AVP brain system include variations in the oxytocin gene (OXT), the oxytocin receptor gene (OXTR), and the arginine vasopressin receptor genes (AVPR1A and AVPR1B). Studies on adolescents’ loneliness so far only examined genetic variations within the OXTR gene. In a longitudinal study, a significant association was found between an OXTR genetic variant (rs53576) and development of loneliness in girls (but not boys) over time (van Roekel, Verhagen, Engels, Goossens, & Scholte, 2013). In an experience sampling method (ESM) study, the same genetic variant was found to be significantly associated with state levels of loneliness (i.e., fluctuations over the day) in girls, but not in boys (van Roekel, Verhagen, Scholte et al., 2013). Significant associations for other single nucleotide polymorphisms (SNPs) in the OXTR gene were also found in relation to social and emotional loneliness in adolescents, but not in adults (Lucht et al., 2009). In a female adult sample, an association between the rs53576 variant of OXTR and emotional loneliness was described, although this finding did not survive correction for multiple testing (Connelly et al., 2014). Despite the fact that these OXTR associations were of small magnitude, the effect sizes were comparable to other candidate gene studies.

Associations between AVPR1A, AVPR1B, and OXT genes and loneliness have not been specifically assessed. However, previous work showed evidence of involvement of these genes in human social functions related to loneliness, such as social reciprocity and empathy (for overview of studies see Ebstein, Knafo, Mankuta, Chew, & Lai, 2012; Feldman, Monakhov, Pratt, & Ebstein, 2016). In adult samples, studies have shown involvement of AVPR1A in pair bonding (Walum et al., 2008), cognitive empathy (Uzefovsky et al., 2015), and amygdala activation while participating in a facial affect task (Meyer‐Lindenberg et al., 2009). Human AVPR1B studies are scarce, but there is some evidence that this gene affects aggression in children (Luppino, Moul, Hawes, Brennan, & Dadds, 2014; Zai et al., 2012), whereas animal studies showed AVPR1B involvement in social motivation and social memory (Wersinger et al., 2004). The OXT gene has been associated with social withdrawal in children (Francis et al., 2016) and epigenetic modification of OXT revealed associations with adult attachment style and recognition of emotional faces (Haas et al., 2016). These alternate phenotypes have all been found to be related to loneliness to a greater or lesser extent (Boivin & Hymel, 1997; Cassidy & Asher, 1992; Davis, 1983; Gardner, Pickett, Jefferis, & Knowles, 2005; Renshaw & Brown, 1993; Rubin, Chen, McDougall, Bowker, & McKinnon, 1995; Vanhalst, Gibb, & Prinstein, 2017; Zysberg, 2012). It is evident that attachment styles and social withdrawal tendencies are closely related to loneliness, but the other phenotypes are also related to loneliness. For example, recognition of facial expressions is a prerequisite for adequate approach or withdrawal behavior when trying to connect to others (Vanhalst et al., 2017). Thus, these studies provide subsidiary evidence for possible involvement of these genes in loneliness. All in all, these findings indicate that genetic variations in the OT‐AVP pathway could partly explain individual differences in loneliness.

However, a recent genome‐wide association study (GWAS) of loneliness did not detect genome‐wide significant associations for variants in the OXTR, AVPR1A, AVPR1B, and OXT genes, nor any other gene (Gao et al., 2017). Although the sample size of over 10,000 participants is quite large, the power was still insufficient to detect variants genome‐wide that are significantly associated with loneliness. It should also be noted that the GWAS included individuals aged 50 years and older, making it hard to draw conclusions regarding the generalizability to adolescent samples. Whereas GWASs are very useful to identify specific SNPs underlying a trait using a hypothesis‐free method (and may therefore inform about novel pathways explaining the biological etiology of a trait), a general disadvantage of this single SNP approach is that complex traits, such as loneliness, are polygenic (i.e., influenced by many genes of small effects). For such polygenic traits, the investigation of candidate gene pathways with known biological function could lead to greater insight in the underlying molecular genetics. The idea underlying this approach is that genes within particular networks interact with each other and that their combined influence would explain a greater portion of the variance than single SNPs do (Purcell et al., 2009).

This study is the first to adopt this approach in relation to adolescent loneliness, with a focus on several genes within the OT‐AVP pathway (OXT, OXTR, AVPR1A, and AVPR1B). More specifically, genetic associations within the OT‐AVP pathway with longitudinal measures of loneliness will be examined in two adolescent samples (N Study1 = 1,103, N Study2 = 404). Single SNP association analyses, gene‐based association tests, and a polygenic risk score (PRS) approach will be used. In single SNP analyses, associations between a genetic variant and the outcome of interest are examined. The gene‐based approach allows to jointly interrogate all SNPs within a gene. Taking the gene as the unit of analysis is expected to increase the statistical power relative to single SNP analysis, because it accounts for multiple independent functional variants while decreasing the number of statistical tests. The strength of a PRS approach is that SNP effect sizes from previous large association results are used to calculate a polygenic risk score per individual that captures the combined effects of multiple genetic variants for a predefined functional pathway. We hypothesize that genetic variations within genes involved in the OT‐AVP pathway could be associated with individual differences in loneliness, both at baseline levels and in the development of loneliness over time.

Because both oxytocin and vasopressin are influenced by gonadal hormones in a sex‐specific manner (Gabor, Phan, Clipperton‐Allen, Kavaliers, & Choleris, 2012), and previous research has documented sex‐specific OXTR gene effects (Kogan et al., 2011; Lucht et al., 2009; van Roekel, Verhagen, Engels, et al., 2013; van Roekel, Verhagen, Scholte, et al., 2013), we included gender as a covariate in the analyses. Additionally, depressive symptoms and social anxiety symptoms were included as covariates, given that these measures are highly interrelated with loneliness in adolescents (Lim, Rodebaugh, Zyphur, & Gleeson, 2016; Vanhalst et al., 2012). However, for both prevalence and genetic architecture, no systematic gender differences have been observed in relation to loneliness (Bartels, Cacioppo, Hudziak, & Boomsma, 2008; Boomsma et al., 2005; Heinrich & Gullone, 2006).

This is the first study investigating the association between genes of the OT‐AVP pathway and loneliness. Using these different genetic approaches could aid in unraveling the genetic background of this highly prevalent phenomenon amongst adolescents and inform future studies into the biological underpinnings of loneliness.

Method

Procedure

Sample 1

Adolescents were recruited through nine high schools which consented to participate in the study. Active informed consents were obtained from both parents and adolescents. The consents inquired separately about participation in the psychological and the biological parts of the study (i.e., DNA collection). The study protocol was approved by a biomedical internal review board (IRB).

Data collection took place between 2012 and 2014. Three annual assessments were conducted during regular class hours, using questionnaires that included a loneliness measure. Saliva samples were collected (using DNA Genotek Oragene kits), under supervision of a researcher at the first measurement (Wave 1). Adolescents were paid approximately 6 U.S. dollars, for their participation at each wave.

Sample 2

High schools were approached through written and personal communication to participate in the present longitudinal study. Six schools agreed to participate and contacted the parents of their first‐year students (comparable to U.S. Grade 7) or provided the researchers with addresses and phone numbers to contact the parents. Parents received information letters and were asked to contact the researchers if they did not want their child to participate. Ethical approval for data collection was obtained from the university's IRB and the collection of saliva for genetic material (at W2) was approved by an independent review board.

Data were collected during regular school hours, in three consecutive years, starting in the first year of high school. In class, adolescents were asked to give signed informed assent, after which they filled out questionnaires about loneliness via computers (W1 and W3). During the second wave of data collection (W2), saliva samples of the adolescents were collected for DNA analysis purposes. Both the parents and adolescents were asked to provide active informed consent for saliva collection. The adolescents were given instructions to provide a small sample of saliva (using DNA Genotek Oragene kits) and were asked to fill out a questionnaire. In exchange for their participation, adolescents could choose a small gift (e.g., a pen).

Sample Characteristics

A total of 2,254 adolescents were invited to take part in this study. The response rate was 49.5%, resulting in 1,116 adolescents that had given active consent and also received active parental consent. After careful inspection of the questionnaires, five subjects were excluded from further analyses due to unreliably filled‐out questionnaires (Janssens et al., 2015). From the 1,111 remaining participants of Wave 1 (mean age at W1 = 13.79 years; SD = .94; 51.0% boys), 986 participants agreed to participate in Wave 2 and 880 participated in Wave 3. A total of 1,103 adolescents agreed to take part in the genetic part of the study (99.3%). In case genetic data (i.e., identity by descent) showed individuals to be related as siblings, one sibling was randomly included in the study (see Figure 1 for flowchart). Ancestry was assessed by the origin of grandparents. A total of 95.3% of the participants reported to be of Caucasian descent (88.8% of grandparents born in Europe, 6.5% Mediterranean Non‐European: predominantly Turkey and Morocco). Other countries were reported by 2.2% of adolescents, and for the remaining 2.5% this information was missing. A small majority of parents completed higher education (58% of the mothers; 52% of the fathers). The remaining parents finished (some years of) high school (30% of mothers; 34% of fathers) or completed primary education (2% of mothers; 2% of fathers). For this sample, we obtained complete data for the loneliness measure, that is, we had no missing values in the data.

Figure 1

Flowchart Sample 1.

Numbers of participants varied across waves, due to new students enrolling in schools and others dropping out of schools, but a subsample of N = 972 adolescents participated in all three waves of data collection. At the first wave of data collection, the participating classes consisted of 1,366 adolescents. Of this eligible group, 89 adolescents were absent at the day of data collection (i.e., they could not participate). From the remaining sample (N = 1,277), only 4.2% did not have parental consent for participation (N = 47) or did not give active consent themselves (N = 7). This resulted in a sample of 1,223 actually participating in Wave 1 (mean age = 12.81 years; SD = .43; 50.0% boys). In the second wave 1,469 adolescents gave consent to participate. From this sample, 28.5% had both parental active consent for DNA collection and provided active consent themselves. All participants from whom we obtained saliva were of European ancestry. In this study, participants that provided saliva samples for DNA analyses were included (N = 418; see Figure 2 for flowchart). The mean age of this subsample was 12.77 years at W1 (SD = .40; 47.0% boys); 20.3% attended preparatory secondary school for technical and vocational training, 33.4% attended preparatory secondary school for college, and 46.3% attended preparatory secondary school for university. For this sample, we obtained complete data for the loneliness measure, that is, we had no missing values in the data.

Figure 2

Flowchart Sample 2. *W2 was time of DNA collection.

Measures

Loneliness

Loneliness was measured in both samples using the Peer subscale of the Loneliness and Aloneness Scale for Children and Adolescents (LACA; Marcoen, Goossens, & Caes, 1987). This 12‐item scale (e.g., “I feel abandoned by my friends”) has a 4‐point Likert scale ranging from 1 (never) to 4 (always). Higher mean scores indicate higher levels of loneliness. Cronbach's alphas across the three waves ranged from α = .90 to α = .91 in Sample 1 and from α = .89 to α = .91 in Sample 2.

Depression

In Sample 1, depressive symptoms were measured with the Center for Epidemiologic Studies Depression Scale (CES‐D; Radloff, 1977). This scale consists of 20 items measuring the prevalence of depressive symptoms during the past week. An example item is “During the past week I thought my life had been a failure.” Respondents provided answers on a 4‐point scale ranging from rarely or none of the days (0 = less than 1 day) to most or all of the days (3 = 5 to 7 days). After reverse coding of three items, higher mean scores indicate a higher prevalence of depressive symptoms in the past week. Cronbach's alphas were good (ranging from α = .78 to α = 93).

In Sample 2, depressive symptoms were measured with the Iowa short form of the Center for Epidemiological Studies Depression scale (CES‐D; Kohout, Berkman, Evans, & Cornoni‐Huntley, 1993). This scale consists of 11 items measuring the prevalence of depressive symptoms during the past week. The items were rated on a 4‐point scale ranging from rarely or never (0 = less than 1 day) to usually or always (3 = 5 to 7 days). After reverse coding of two items, higher mean scores indicate a higher prevalence of depressive symptoms in the past week. Cronbach's alphas were good, ranging from α = .81 to α = .85

Social anxiety

In Sample 1, social anxiety was assessed using the short version of the Social Anxiety Scale for Adolescents (SAS‐A; La Greca & Lopez, 1998). This questionnaire assesses social anxiety with 12 questions that are scored using a 5‐point Likert scale ranging from 1 (not at all) to 5 (all the time), with higher mean scores reflecting higher levels of anxiety. Cronbach's alphas showed high reliability across all three waves (α = .92).

In Sample 2, social anxiety was measured using the Social Phobia subscale of the Screen for Child Anxiety Related Emotional Disorders (SCARED; Bodden, Bogels, & Muris, 2009). This questionnaire consists of 9 items scored on a 3‐point scale ranging from almost 1 (never) to 3 (often) (e.g., “I don't like to be with people I don't know well”). Higher mean scores indicate a higher prevalence of social anxiety symptoms (Cronbach's alphas ranged from .80 to .85).

Gene Selection and Genotyping

The genetic variants genotyped in this sample were selected based on an extensive literature search. Candidate genes that were suspected to be involved in neurotransmitter systems were identified. This selection was extended using information available on known neurotransmitter‐related protein interaction networks and pathways (Franceschini et al., 2013). A total of 344 genes represented by 6,325 unique SNPs were selected for genotype analysis, representing pathways related to oxytocin, serotonin, dopamine, HPA‐axis, GABA, glutamate, choline, and noradrenergic neurotransmission, and the clock‐gene related network. The selection included known candidate SNPs, as well as tagging SNPs identified using Haploview (Barrett, Fry, Maller, & Daly, 2005) allowing for analysis of the underlying linkage disequilibrium (LD) structure of the genes. SNPs were genotyped using an Illumina iSelect custom beadchip and the Illumina iScan. Only SNPs within the OXTR, OXT, AVPR1A, and AVPR1B genes were included in this study (see Table S1 in the online Supporting Information).

For quality control purposes, all SNPs were manually reviewed using GenomeStudio (Genotyping Module), according to the developer's guidelines (Illumina, 2008). Further quality control was performed according to standard GWAS protocols as described by Anderson et al. (2010) and Purcell et al. (2009). Details can be found in Appendix S1 in the online Supporting Information. After quality control on the SNPs related to the genes of interest, the sample comprised 1,030 individuals and 76 SNPs within the OXTR, OXT, AVPR1A, and AVPR1B genes (see Table S1) (see also Van Assche et al., 2017).

Genome‐wide genotyping was performed using the Infinium PsychArray‐24 v1.1 BeadChip, containing 265,000 tag SNPs, 245,000 exome markers, and 50,000 additional markers associated with common psychiatric disorders (http://www.illumina.com/products/psycharray.html). Data on the BeadChip were cleaned using the directions provided by the Broad Institute of MIT and Harvard (Purcell et al., 2009). To attain SNPs from the genome‐wide set that lie within the OXTR, OXT, AVPR1A, and AVPR1B genes, we used the UCSC genome browser (assembly GRCh37/hg19) to define the borders of each gene. All SNPs within these regions, including variants within 20 kilobase (kb) flanking region of each gene, were included in this study (see Table S1).

Quality control was performed using GenomeStudio and included several steps (see Appendix S1 for details). For this adolescent sample, N = 393 remained for analysis and 25 SNPs within the AVPR1B, AVPR1A, OXTR, and OXT genes were included (see Table S1 for the SNP list and see Table S2 in the online Supporting Information for identical SNPs in both samples).

Strategy of Analyses

The effects of different oxytocin‐ and vasopressin‐related genetic variants on loneliness measures were examined with SNP‐based, gene‐based, and polygenic risk score methods. Before considering the genetic effects, we first used latent growth curve modeling (LGCM) in Mplus (Muthén & Muthén, 1998–2007) to estimate both the baseline level of loneliness (i.e., intercept) and the change in loneliness over time (i.e., slope). With this approach, individual change across the three waves was determined for each adolescent (Duncan, Duncan, & Strycker, 2006). By default, Mplus uses a full information maximum likelihood (FIML) estimation approach to handle missing values (Muthén & Muthén, 1998–2007). To evaluate model fit, commonly used measures of fit were used, being the comparative fit index (CFI, with a cut‐off value of > .95) and the root mean square error of approximation (RMSEA, with a cut‐off value of < .06; Hu & Bentler, 1999). Analyses revealed that linear models described the data best. Therefore, the intercept and linear slope estimates were considered to describe baseline levels of loneliness and the development of loneliness over time, respectively.

In a second step, we examined the SNP‐based associations with both the intercept and slope of loneliness using association tests in PLINK (Purcell et al., 2007), with gender, depression (at W1), and social anxiety (at W1) as covariates (Lim et al., 2016; Vanhalst et al., 2012). We report adjusted p‐values as provided by PLINK to indicate significance levels after correction for multiple testing, using Bonferroni correction and Benjamini and Hochberg step‐up false discovery rate (FDR) control (p < .05) methods (Benjamini & Hochberg, 1995). Given that PLINK does not account for LD structures between SNPs in the multiple testing metrics, we also calculated the number of independent SNPs based on LD structures (with R threshold of 0.25). This approach provides a more accurate multiple testing corrected p‐value (corrected p‐value of .0025 in Sample 1, based on 20 independent SNPs and corrected p‐value of .0042 in Sample 2, based on 12 independent SNPs).

Third, we examined gene‐based associations with both loneliness growth parameters (i.e., intercept and slope) by using the set‐based option in PLINK (Purcell et al., 2007). This option allows for simultaneous testing of all independent SNPs (i.e., SNPs that are not in LD with other SNPs) within a gene, yielding empirical p‐values that account for the number of SNPs within the set and the number of permutations (in this case 10,000). We used default parameter values as defined by PLINK, being R 2 < 0.5, p‐value < .05, and the maximum number of SNPs = 5. The R 2 command describes the LD threshold. In this case, SNPs that are in LD (above a threshold of R 2 = 0.5) with previous SNPs are removed from the analysis. The last default command indicates that a maximum of five independent SNPs are included per set, under the notion that the most significant SNP is selected first and subsequent SNPs are selected on the basis of decreasing significance (and only if pSNP < .05). For both samples, we had four sets of genes (AVPR1A, AVPR1B, OXT, and OXTR).

Lastly, we calculated polygenic risk scores (PRS) based on a recent GWAS on loneliness (Gao et al., 2017) (see https://www.med.unc.edu/pgc/results-and-downloads). First, we identified which SNPs from the GWAS were present in our samples (54 for Sample 1 and 14 for Sample 2). Second, we estimated LD structures for both samples, using the clump command as implemented in PLINK (Purcell et al., 2007). We considered all SNPs and used an R threshold of 0.25 and a 250 kb window to cluster SNPs. With this procedure, SNPs in high LD with other SNPs were excluded from further analysis. After clumping, 20 SNPs remained in Sample 1, and 12 SNPs remained in Sample 2 (see Table S2). Third, the individual polygenic risk scores were calculated as the number of risk alleles (either 0, 1, or 2) multiplied by their respective unstandardized weights from the GWAS loneliness scores (Gao et al., 2017). The summed scores for each individual were used in linear regression models to estimate the effect of the polygenic risk scores on both the intercept and slope of loneliness.

Results

Descriptive Statistics

Participants of Sample 1 were somewhat older at Wave1 (M1 = 13.79 years; SD = .94) than participants of Sample 2 (M1 = 12.81 years; SD = .43), but they were of comparable age again at Wave 3 (Sample 1; M3 = 15.74 years; SD = .92, Sample 2; M3 = 15.14 years; SD = .46). Correlations between ages at Wave 1 and the psychological constructs were examined but they did not reveal significant associations (see Table S3 in the online Supporting Information; The same was true for correlations with ages at Waves 2 and 3, data not shown).

The annual loneliness scores showed a small, but significant, increase over the three consecutive years. In both samples and at all three waves, female participants scored significantly higher than males on loneliness, depressive symptoms, and social anxiety (see Table 1). There was one exception; boys and girls did not differ at the baseline measurement of loneliness in Sample 1. Correlations between the annual loneliness measures and the control measures (i.e., depressive symptoms and social anxiety at the first wave) were all positive and significant (see Table 2). Different social anxiety scales were administered to the two samples, but, correlations between social anxiety and the other measures were highly similar in both samples.

Table 1

Means (Standard Deviations) and Gender Differences for Age, Loneliness (Intercept and Slope), Depressive Symptoms, and Social Anxiety Symptoms in Both Samples at the Annual Waves

Variable	Sample 1 N = 1,030	Gender Difference (t)	Sample 2 N = 393	Gender Difference (t)
Gender (N, % boys)	557 (51.2%)		185 (47.1%)
Age W1	13.79 (0.94)	−2.87**	12.76 (0.39)	ns
Age W2	14.83 (0.92)	−2.64**	13.85 (0.43)	ns
Age W3	15.74 (0.92)	−2.19*	15.14 (0.45)	ns
Loneliness W1	1.54 (0.56)	ns	1.52 (0.48)	−2.22*
Loneliness W2	1.56 (0.54)	−3.44**	1.58 (0.52)	−3.37**
Loneliness W3	1.58 (0.55)	−5.15***	1.58 (0.53)	−3.12**
Intercept loneliness	1.54 (0.36)	−2.57*	1.52 (0.31)	−3.50**
Slope loneliness	0.02 (0.11)	−4.81***	0.02 (0.12)	−2.16*
Depressive symptoms W1	0.86 (0.34)	−4.17***	0.53 (0.44)	−3.61***
Depressive symptoms W2	0.57 (0.49)	−6.23***	0.52 (0.46)	−5.18***
Depressive symptoms W3	0.57 (0.51)	−5.78***	0.51 (0.42)	−5.72***
Social anxiety W1	2.40 (0.79)	−4.48***	1.62 (0.42)	−3.72***
Social anxiety W2	2.55 (0.81)	−5.35***	1.61 (0.44)	−7.20***
Social anxiety W3	2.53 (0.79)	−5.63***	1.61 (0.45)	−5.59***

W1 = Wave 1; W2 = Wave 2; W3 = Wave 3; ns, non significant.

*p < .05; **p < .01; ***p < .001.

Table 2

Correlations Between Annual Loneliness Scores and Depressive Symptoms (W1) and Social Anxiety Symptoms (W1) for Both Samples

Variable	1.	2.	3.	4.	5.
1. Loneliness W1		.50**	.46**	.50**	.48**
2. Loneliness W2	.54**		.62**	.35**	.40**
3. Loneliness W3	.47**	.62**		.26**	.40**
4. Depressive symptoms W1	.36**	.25**	.23**		.45**
5. Social anxiety W1	.53**	.40**	.44**	.35**

Below the diagonal the correlations for Sample 1, above the diagonal for Sample 2.

**p < .01.

Baseline Loneliness and the Development of Loneliness Over Time

Latent growth curve modeling revealed that both the intercept and linear slope of loneliness were significant in Sample 1 (β0 = 1.54, p < .001; β1 = .02, p < .05). These results indicated that adolescents on average scored 1.54 on the loneliness scale at baseline, and that the level of loneliness slightly increased over time (χ2 [df = 3, N = 1,108] = 779.841, CFI = 1.00, and RMSEA = 0.0). In Sample 2, a very similar pattern was observed. Again, both the intercept and slope of adolescents’ loneliness were significant (β0 = 1.53, p < .001; β1 = .03, p < .05). On average, participants scored 1.53 on the baseline measure and loneliness showed a slight increase over time (χ2 [df = 3, N = 418] = 275.385, CFI = 0.995, and RMSEA = 0.059).

SNP‐Based Associations With Adolescents’ Loneliness

Within Sample 1, the vast majority of oxytocin and vasopressin SNPs were unrelated to either the intercept or slope of the loneliness score. One SNP (rs918316) within the OXTR gene was marginally associated with the slope (see Table 3). This SNP remained significant after adding the covariates (i.e., gender, depressive symptoms, and social anxiety symptoms), but the association did not hold after strict correction for multiple testing (see Table S4 in the online Supporting Information for all SNP results).

Table 3

SNP‐Based Associations (if Unadjusted p < .05) for Development of Loneliness in the Two Samples, With and Without Covariates

Sample	SNP	Gene Name	Univariate (Unadjusted)		Multivariate (with Covariates)
			B slope	p‐value	B slope	p‐value
Sample 1	rs918316	OXTR	−.02 (.01)	.02	−.02	.03
Sample 2	rs6793234	OXTR	.03 (.01)	.03	.05	.003*
	rs75775	OXTR	−.02 (.01)	ns	−.03	.01

*SNP is significant after correction for multiple testing. Bonferroni p‐value < .004 (= alpha of .05 divided by 12 independent SNPs). Included covariates: Gender, depressive symptoms, and social anxiety symptoms.

In Sample 2, one SNP (rs6793234) within the OXTR gene was significantly associated with the slope of the loneliness score (see Table 3). Following inclusion of the covariates in the analyses, the association with this OXTR SNP remained significant and another OXTR SNP (rs75775) became significantly associated with the slope of loneliness. None of the SNPs remained significant after strict correction for multiple testing. However, using the less stringent multiple testing correction based on the number of independent SNPs (n = 12, p < .0042, as found in the clumping analyses) revealed that rs6793234 did pass this threshold and remained significant (see Table S5 in the online Supporting Information for all SNP results).

From the SNPs showing significance, only rs75775 was present in both samples. However, in Sample 1, this SNP was not associated with the development of loneliness over time (β = .0004, SE = .007, p = .95).

Gene‐Based Associations With Adolescents’ Loneliness

In Sample 1, the gene‐based analyses showed a significant association for the AVPR1A set with the intercept of loneliness when controlling for gender, depressive symptoms, and social anxiety, with a set‐based empirical p‐value of < .05. However, taking into account that multiple genes were analyzed simultaneously, the Bonferroni corrected p‐values should be set at p < .013. The AVPR1A gene did not pass that threshold (set‐based empirical p‐value = .036). The AVPR1A gene set was not associated with the slope of loneliness. The other three gene sets were unrelated to the intercept and slope of loneliness.

In Sample 2, the OXTR gene set was significantly associated with the slope of loneliness when controlling for gender, depressive symptoms, and social anxiety, with a set‐based empirical p‐value of .044. However, taking into account that multiple genes were analyzed simultaneously, the Bonferroni corrected p‐value threshold should be p < .013. The OXTR gene set was not associated with the intercept of loneliness, neither were the other gene sets.

Polygenic Risk Score Associations With Adolescents’ Loneliness

Finally, regression analyses with polygenic risk scores (PRS) for the four OT‐AVP related genes as predictor for adolescents’ loneliness did not reveal significant associations in one of the samples neither in the uncontrolled models nor in the models including gender, depressive symptoms, and social anxiety as covariates.

Discussion

In this study, involvement of OT‐AVP–related genes in adolescents’ loneliness was examined in two distinct samples, using SNP‐based and gene‐based tests and a polygenic risk score approach. Although not finding strong evidence for genetic main effects, we did find OXTR SNP‐ and gene‐set associations with the slope of loneliness in both samples, that did not survive correction for multiple testing. In addition, a significant association of the AVPR1A gene on the intercept of loneliness was found in Sample 1. Again, this finding did not survive correction for multiple testing. Further, the polygenic risk score approach did not provide strong evidence for involvement of the OT‐AVP pathway in adolescents’ loneliness either. Altogether, our approach has not consistently shown that genes within the OT‐AVP pathway modulate the intercept or slope of loneliness in adolescence.

In both samples of adolescents, loneliness slightly increased over three successive years. This finding is in contrast with other studies examining the development of loneliness over time, which mainly described decreasing levels of loneliness in adolescence over the years after an initial peak at age 13 (Harris, Qualter, & Robinson, 2013; Qualter et al., 2013; Vanhalst, Goossens, Luyckx, Scholte, & Engels, 2013). It is difficult to compare the different loneliness trajectories, because studies differed in ages at baseline, intervals between waves, and total study period. The only study also using annual assessments, with comparable age at baseline, found decreasing loneliness levels (van Roekel, Scholte, Verhagen, Goossens, & Engels, 2010). However, in the latter study adolescents were followed up for five consecutive years. Given that the loneliness scores in our samples were only slightly increasing over the first three assessments, and given that the steepest decrease in loneliness scores in van Roekel et al. (2010) were observed between the fourth and fifth measurement (which typically would be the year after our last assessment), an overall (nonlinear) decrease in loneliness may also have been found in our samples if followed in subsequent years.

Interestingly, our finding of involvement of the OXTR gene in the development of loneliness over time is consistent with a longitudinal study in adolescents also showing significant OXTR SNP effects for the development of loneliness, but not for baseline loneliness (van Roekel, Verhagen, Engels, et al., 2013). The previous study used a single‐SNP approach and examined rs53576 only. This SNP has also been described to be significantly associated with social and emotional loneliness (Lucht et al., 2009). This SNP was included in Sample 1 (but not in Sample 2) but did not show a significant association with loneliness in our sample (β = −0.003, SE = 0.01, p = .62).

The rs75775 SNP in the OXTR promoter region was measured in both samples, but only showed significance in Sample 2. In that same sample, the OXTR gene‐set test also provided evidence for significant association with the development of loneliness over time. Other significant associations for SNP rs75775 have been observed for autism and related social characteristics (Wang et al., 2009) and, more interestingly, this SNP has also been associated with prosocial behavior (Apicella et al., 2010), a trait that has been described in relation to loneliness (e.g., loneliness could either decrease or increase prosocial behavior due to enhanced self‐focus or the urge to restore relationships (Spithoven, Vanhalst, Lodder, Bijttebier, & Goossens, 2017; Twenge, Baumeister, DeWall, Ciarocco, & Bartels, 2007; Woodhouse, Dykas, & Cassidy, 2012). This SNP was not significantly associated with loneliness in the GWAS comprising adults aged 50 and over (Gao et al., 2017). For the other significant SNPs in our study (rs918316 and rs6793234), no previous associations have been documented.

To conclude, some OXTR SNPs and the OXTR gene set were marginally significant predictors of loneliness, but significance did not hold after correction for multiple testing. Our initial hope was, of course, that the inclusion of more OXTR SNPs than in previous studies and the examination of the combined effect of multiple SNPs would assist in unraveling OXTR gene effects in relation to loneliness. However, given previous studies describing associations with loneliness (Lucht et al., 2009; van Roekel, Verhagen, Engels, et al., 2013; van Roekel, Verhagen, Scholte, et al., 2013) and associated traits that could well be involved in the etiology or maintenance of loneliness (Feldman et al., 2016; Haram et al., 2015; Israel et al., 2009; Rodrigues, Saslow, Garcia, John, & Keltner, 2009), we consider this a trend finding worthy of further investigation.

The gene‐set tests also showed involvement of the AVPR1A gene on baseline loneliness (in Sample 1 only). This is the first study demonstrating this association. Previous studies have suggested that genetic variation in the promotor region of AVPR1A could play a role in impairments in social connectedness and interpersonal relationships (Walum et al., 2008). A possible factor that could lead to such impairments in social relations is the ability to correctly process and identify emotional stimuli. Indeed, a neuroimaging study showed that genetic variation in the AVPR1A promoter was associated with amygdala activation when viewing faces expressing negative affect (Meyer‐Lindenberg et al., 2009). Also, the ability to correctly identify static emotional expressions was affected by this promotor variant (Golimbet, Alfimova, Abramova, Kaleda, & Gritsenko, 2015), as was cognitive empathy (Uzefovsky et al., 2015). These findings illustrate its role in the processing and recognition of visual emotional stimuli. In addition, AVPR1A (Avinun et al., 2011) has also been associated with prosocial behavior. These previously associated traits are more or less distally related to loneliness, that is, a diminished capacity for emotion recognition could lead to less adequate behavioral reactions during social interactions, with the risk of increasing loneliness (Vanhalst et al., 2017; Woodhouse et al., 2012). Again, although showing only weak evidence for involvement of the AVPR1A gene in loneliness, there is suggestive evidence that this gene may modulate precursors of loneliness and could be explored in greater depth in the future.

Strengths and Limitations

A definite strength of this study is the possibility to examine the same biological pathway in two independent samples of developing adolescents, with similar demographic characteristics, similar measurements, and a similar development of loneliness over time. In addition, the combination of SNP‐based, gene‐based, and polygenic risk score approaches constitutes a powerful design, providing step‐by‐step insight in the respective roles of OT‐AVP pathway genes. Another strength is that we included several covariates, to ensure that the possible effects were uniquely associated with loneliness. Finally, the use of LGCM statistics to model loneliness over time provides a unique insight into the development of loneliness during adolescence that moves beyond the mere examination of cross‐sectional associations. Despite these strengths, an obvious limitation is the size of both samples. The samples could be regarded as too small to have sufficient power to detect relevant genetic effects, which is often the case when examining complex heterogeneous phenotypes. However, the finding that significant signals within the same gene were obtained in two independent samples is a clear suggestion that OXTR could be involved in loneliness. The assumed lack of power could also explain why the OXTR and AVPR1A gene findings did not reach statistical significance. Another limitation is that different genetic arrays were used to genotype the markers in the OT‐AVP genetic pathway in both samples. In the first sample, a high‐density chip was used to specifically genotype these markers, whereas the second sample provided genome‐wide genetic information in which these particular genes were less densely covered. This difference in overall approach may explain the absence of overlap in genetic variants across both samples. A last limitation concerns the PRS approach. The overall variance explained by the combined SNPs in the polygenic score was low (up to 1%). Polygenic scores rely heavily on the accuracy of the SNP effect estimates in the discovery sample (Wray et al., 2013); the summary statistics we used were based on a relatively small discovery sample (N = 10,760). Prediction accuracy will increase when larger samples sizes become available.

Implications for Future Research

Aside from sample size and replication considerations, there are several other avenues to explore. The first avenue concerns the level of phenotype measurement. Although the use of multiple measurements per individual as in this study leads to more intense phenotyping compared to cross‐sectional measures, and provides an estimate of the development of loneliness over time, this measure also comprises adolescents with stable low levels of loneliness. It would be interesting to focus only on adolescents reporting on the more severe ends of loneliness, that is, chronic or persistent loneliness, analogue to studies showing stronger effects for severe forms of psychopathology (versus mild forms: Uher et al., 2011). Next to this, though providing a completely different phenotype, an alternative measurement level enabling deep phenotyping of loneliness could be the so‐called moment‐to‐moment fluctuations in levels of loneliness. These consist of repeated loneliness assessments with very short time lags in between over a short period of time. This could for example provide insight in the relations between genes and the amount of variability in daily loneliness fluctuations or in the number of peaks above an individual's mean (van Roekel, Verhagen, Scholte, et al., 2013).

Secondly, given the OT‐AVP gene involvement in various traits that could play a role in both the etiology and maintenance of loneliness, it could be interesting to examine these as possible endophenotypes for loneliness. Examples of endophenotypes could be amygdala responsiveness to emotional faces to social integration possibilities and emotional withdrawal tendencies (Feldman et al., 2016; Haas et al., 2016; Meyer‐Lindenberg et al., 2009).

Aside from genetic factors, it is interesting for future research to incorporate social environmental factors such as social support or the perception of company, both in relation to trait levels and state levels of adolescents’ loneliness (van Roekel et al., 2010; van Roekel, Verhagen, Scholte, et al., 2013). Certain gene by environment interactions could further unravel the underlying dynamics in adolescents’ loneliness.

Another approach that could elucidate the genetic mechanisms underlying loneliness could be found in the field of epigenetics. The idea is that certain (environmental) influences (e.g., experienced trauma or nutrition) could affect gene activity by a process named methylation (Kumsta, Hummel, Chen, & Heinrichs, 2013). DNA methylation could increase the individual's vulnerability for psychopathology (McGowan & Szyf, 2010; McGowan et al., 2009). For example, epigenetic studies of the OXTR gene (Kumsta et al., 2013) have shown that OXTR gene alterations were associated with social perception processes and psychosocial stress sensitivity, both of which have been associated with loneliness (e.g., Hawkley, Browne, & Cacioppo, 2005; Vanhalst et al., 2017). In combination with an animal study showing that social isolation (which is on the extreme end of feeling lonely) led to both epigenetic and phenotypic changes in mice (Siuda et al., 2014), this indeed could be a relevant and interesting avenue to pursue.

Conclusion

We found involvement of the AVPR1A and OXTR genes in baseline levels and development of loneliness over time in two distinct adolescent samples using a variety of genetic approaches. After strict correction for multiple testing, genetic associations were no longer significant. This could reflect absence of genetic associations within the OT‐AVP pathway and loneliness. However, we argue that this absence of strong evidence for involvement of OT‐AVP pathway genes should not discourage researchers from further examining these interesting and plausible candidate genes in relation to adolescent loneliness. We recommend using deep phenotyping of loneliness, include environmental factors (especially factors indicating social resources), to consider epigenetic studies, and to examine possible endophenotypes.

Appendix S1 Genotyping Information for Sample 1 and Sample 2.

Table S1 SNP List and Numbers of SNPs Included in Each Gene in the Two Samples After Quality Control Steps.

Table S2 SNP Lists After Clumping for Both Samples.

Table S3 Correlations Between Ages at Wave 1 and Annual Loneliness Scores, Depressive Symptoms, and Social Anxiety Symptoms for Both Samples.

Table S4 All SNP‐Based Results for Sample 1.

Table S5 All SNP‐Based Results for Sample 2.

Click here for additional data file.