The oxytocin signalling gene pathway contributes to the association between loneliness and cardiometabolic health

Increasing evidence has shown adverse effects of loneliness on cardiometabolic health. The neuromodulator and hormone oxytocin has traditionally been linked with social cognition and behaviour. However, recent implications of the oxytocin system in energy metabolism and the overrepresentation of metabolic issues in psychiatric illness suggests that oxytocin may represent a mechanism bridging mental and somatic traits. To clarify the role of the oxytocin signalling system in the link between cardiometabolic risk factors and loneliness, we calculated the contribution of single nucleotide polymorphisms (SNPs) in the oxytocin signalling pathway gene-set (154 genes) to the polygenic architecture of loneliness and body mass index (BMI). We investigated the associations of these oxytocin signalling pathway polygenic scores with body composition measured using body magnetic resonance imaging (MRI), bone mineral density (BMD), haematological markers, and blood pressure in a sample of just under half a million adults from the UK Biobank (BMD subsample n = 274,457; body MRI subsample n = 9796). Our analysis revealed significant associations of the oxytocin signalling pathway polygenic score for BMI with abdominal subcutaneous fat tissue, HDL cholesterol, lipoprotein(a), triglycerides, and BMD. We also found an association between the oxytocin signalling pathway polygenic score for loneliness and apolipoprotein A1, the major protein component of HDL. Altogether, these results provide additional evidence for the oxytocin signalling pathway's role in energy metabolism, lipid homoeostasis, and bone density, and support oxytocin's complex pleiotropic effects.


Introduction
Loneliness is a common human experience, described as the aversive state that arises from having fewer and less satisfying social relationships than desired (Peplau, 1982). Loneliness is a persistent state for 15-30 % of the general population (Heinrich and Gullone, 2006). The pervasive adverse effects of loneliness became particularly relevant during the COVID-19 pandemic as physical distancing rules resulted in a decline of in-person social contact, leading to an increase in the prevalence of loneliness (Groarke et al., 2020). Along with causing psychological distress, persistent loneliness is also strongly associated with negative somatic health outcomes and mortality (Holt-Lunstad et al., 2015). In particular, research has linked loneliness with cardiovascular mortality (Valtorta et al., 2016) and metabolic syndrome (Henriksen et al., 2019), a condition commonly thought to comprise four disease risk factors: abdominal obesity, dyslipidaemia, high blood pressure, and impaired glucose homoeostasis. A recent phenome-wide investigation of the association between a genetic predisposition to loneliness and health outcomes found that a genetic propensity to loneliness is associated with increased risk for cardiovascular disease and type-2 diabetes. Elevated triglycerides and reduced HDL, two elements of dyslipidaemia, were also associated with a predisposition to loneliness, along with elevated ☆ Dirk Hellhammer Award Paper 2022. body mass index (BMI) and body fat (Abdellaoui et al., 2019).
The oxytocin system is well known for its putative association with social dysfunction (Alvares et al., 2016), but growing evidence also points to its role in metabolic regulation due to its reported effects on appetite, caloric intake, lipolysis, insulin sensitivity, and body weight (McCormack et al., 2020). Oxytocin administration has been shown to reduce caloric intake in humans (Lawson, 2017;McCormack et al., 2020;Ott et al., 2013;Striepens et al., 2016;Thienel et al., 2016) and to reduce body weight in obese rhesus monkeys (Blevins et al., 2015) and rats (Deblon et al., 2011) via decreased food intake and increased energy expenditure. Oxytocin's reported effects on weight loss are also thought to occur via its effects on adipocytes, which specifically target adipose tissue leaving lean mass undiminished (Altirriba et al., 2014;Boland and Goren, 1987). Adipose tissue also helps regulate bone integrity (Amri and Pisani, 2016;Lee and Karsenty, 2008), which influences body weight control, glucose homoeostasis, and fat metabolism through the production of several endocrine factors (Amri and Pisani, 2016;Ferron et al., 2008;Lee and Karsenty, 2008). Notably, both osteoblasts and osteoclasts (respectively involved in bone synthesis and remodelling) express oxytocin receptors (Elabd et al., 2008), and bone mineral density (BMD) has been linked to cardiovascular disease (CVD) (Park et al., 2021;Veronese et al., 2017). This emerging body of research has contributed to the suggestion that oxytocin system dysregulation may be a mechanism common to both social and metabolic dysfunction (Quintana et al., 2017).
Candidate gene studies have been a popular approach for investigating oxytocin's role in social behaviour and loneliness (Lucht et al., 2009;Roekel et al., 2013), with research primarily focusing on variants in the oxytocin receptor gene (OXTR). Genetic studies can be a useful way to investigate oxytocin's influence on multiple phenotypes. However, given the small influence single gene variants are likely to exert on a given trait, conventional candidate gene studies do not typically have sufficient statistical power to identify susceptibility with small effects, which has led to difficulties replicating many candidate gene study findings (Harden, 2021). As there are over 150 genes in the oxytocin signalling pathway, a polygenic approach can increase sensitivity for the analysis of the role of the oxytocin signalling system as a possible common denominator linking loneliness, body composition, and metabolic syndrome risk factors.
Fat deposit location influences metabolic risk, with fat surrounding organs in the abdominal cavity (i.e., visceral adipose tissue), along with adipose infiltration (i.e. ectopic fat tissue) most associated with metabolic disturbances (Bergman et al., 2006;Després and Lemieux, 2006). As BMI is not particularly sensitive to the distribution of body fat, alternative techniques, such as body magnetic resonance imaging (MRI) segmentation, are required to measure the individual contributions of different fat depots (Beck et al., 2021;Gurholt et al., 2021;Linge et al., 2019Linge et al., , 2018). In the current study, we use an oxytocin pathway polygenic score (PGS oxt )  to evaluate the genetic contribution of SNPs belonging to 154 genes attributed to the oxytocin signalling pathway on cardiometabolic risk factors, including body fat distribution and infiltration indexed via body MRI, haematological markers, blood pressure, and bone mineral density in a large sample of adults from the UK Biobank. Oxytocin-specific polygenic scores for loneliness and BMI were calculated and then associated with these cardiometabolic risk variables and the differential expression of oxytocin-pathway genes was also evaluated to corroborate findings from these associations.

Population
The UK Biobank is a prospective cohort including approximately half a million participants (54 % female) in the age range of 40-69 years (Sudlow et al., 2015), recruited between 2006 and 2010. All participants provided signed informed consent. The UK Biobank has Institutional Review Board approval from North West Multi-centre Research Ethics Committee and its Ethics Advisory Committee (https://www.ukbiob ank.ac.uk/ethics) oversees the UK Biobank Ethics & Governance Framework (Miller et al., 2016). Specific details regarding recruitment and data collection procedures have been previously published (UK Biobank, 2007). Genotyping and initial quality control (QC) of ~ 96 million variants were performed by Affymetrix (2014). Further QC procedures have been described previously (Bycroft et al., 2018). Only participants of European British descent were selected for the main analysis. The present research has been conducted using the UK Biobank Resource under Application no. 27412. An overview of the variables used in the main analysis is available in Table 1. Data were used from 488,377 genotyped individuals, including a subset of 9796 participants with available body composition profiles obtained from body MRI and 274,457 with heel bone mass density measurements.

Oxytocin pathway polygenic score (PGS oxt )
We obtained the genotypes of 458,775 individuals of European ancestry meeting QC requirements. We calculated two oxytocin pathway polygenic scores (PGS oxt ) using PRSice-2 (version 2.3.3) (Euesden et al., 2015), by limiting the calculation to SNPs belonging to the oxytocin signalling pathway, adopting a methodology similar to that of a previous report (Darst et al., 2017). Using BMI (Locke et al., 2015) and loneliness (Gao et al., 2017) GWAS summary statistics, we calculated 2000 PGS oxt scores with thresholds ranging from 5 × 10 − 8 to 1, in increments of 0.001. Using a permutation approach (calculating the empirical p-value by obtaining the p-value of association of the best p-value threshold and a distribution of best p-value thresholds under the null, shuffling the target phenotype 10,000 times) included in PRSice-2 (Euesden et al., 2015), we determined optimal (controlling for Type 1 error) PGS oxt thresholds of p < .035 for BMI (resulting in 295 SNPs) and p < .011 for loneliness (resulting in 92 SNPs), and used PGS oxt 's computed at these thresholds for the main analysis, after a transformation to z-scores. The 10 first principal components of the sample's variance-standardised relationship matrix were annotated alongside the PGSs to account for population stratification in the sample itself. We also mapped the SNPs included in each PGS oxt using the snp2gene component of FUMA (Watanabe et al., 2017) and estimated their enrichment in KEGG pathways using DAVID (Huang et al., 2009a(Huang et al., , 2009b, confirming enrichment for the oxytocin signalling pathway. This procedure helped make sure that the subsets of SNPs that are selected by thresholding the PGS oxt 's are specifically tied to the oxytocin pathway. The gene lists (using the standard gene symbol nomenclature) and the output of the enrichment analysis are available in the Supplementary materials (Supplementary Tables 1 and 2, Supplementary materials 1-2).

Body composition profiles
To analyse body composition, we used data derived from whole-body MRI supplied to the UK Biobank by AMRA (Advanced MR Analytics AB, AMRA, Sweden). Detailed methods on how these fields were derived are described in previous work (Linge et al., 2018). In short, body and liver MRI was acquired on a 1.5 T Siemens MAGNETOM Aera scanner using a body dual-echo Dixon Vibe protocol 4.5 and a single-slice multi-echo gradient Dixon acquisition 5, respectively (Borga et al., 2015;Wilman et al., 2017). Then, a supervised automated segmentation tool was used to segment different adipose tissue compartments Linge et al., 2019Linge et al., , 2018West et al., 2016).
We selected several variables describing the distribution of abdominal fat based on their metabolic significance (Britton and Fox, 2011;Linge et al., 2019Linge et al., , 2018: visceral adipose tissue volume [VATi -L, normalised by height squared following previous literature (Beck et al., 2022;Gurholt et al., 2021;Heymsfield et al., 2007)], abdominal subcutaneous adipose tissue volume (ASATi -L, also normalised by height squared), percentage of muscle fat infiltration (MFI) and liver proton density fat fraction (PDFF). In addition, we used the composite variables of abdominal fat ratio (total abdominal fat divided by total abdominal fat and thigh muscle volume) to assess the distribution of fat and muscle volume (Linge et al., 2018), and weight-to-muscle ratio (WMR; kg/Lbody weight divided by thigh muscle volume) as an index of the ability of individuals to carry their weight (Linge et al., 2018).

Bone mineral density (BMD) measurement
To evaluate BMD, we extracted the Heel BMD T-score calculated from the calcaneal BMD ultrasound measurement and the value normally expected in someone of the same sex. The units of the T-score are the number of standard deviations that the bone density is above or below the standard (Frost et al., 2001) and are used to assess and predict bone health, osteopenia and osteoporosis. We chose this measurement over raw BMD scores for two main reasons-its broad clinical application (Knapp et al., 2001) and the standardisation by sex, which better controls for potential sex differences that can influence results in oxytocin research (Winterton et al., 2020).

Blood biochemistry
We selected several cardio-metabolically relevant blood assays: total cholesterol (mmol/L), HDL cholesterol, LDL cholesterol (mmol/L), triglycerides (mmol/L), apolipoprotein A1 and B (g/L), and lipoprotein(a) (nmol/L) to evaluate dyslipidaemia and lipid homoeostasis; and glycated haemoglobin (HbA1c; mmol/mol) as an index for impaired glucose homoeostasis, commonly used in assessing glycaemic trends in a population over time (Beck et al., 2017). In addition, we selected pulse pressure (calculated as the difference between systolic and diastolic pressure, averaged from two automated measurements a few moments apart, including manual measurements when automatic measurement was not available; mmHg) as an index of cardiovascular health, as it appears to be a better predictor of both cardiovascular and total mortality in older adults than either mean, systolic or diastolic pressure (Glynn et al., 2000).

Statistical analysis
All statistical analyses were done using R version 4.0.0 (2020-04-24) (R Core Team, 2020). The code for the analysis performed is made available at https://osf.io/hy5aw/. First, to analyse the relationship between variables, we computed the pairwise correlations between phenotypical variables, and used complete linkage hierarchical clustering on the resulting correlation matrix to identify variable clusters. To analyse the relationships between continuous variables, we fitted linear regression models, controlling for sex, BMI (following Linge et al., 2018), age, age squared, and the 10 top principal components from the variance-standardised relationship matrix to account for population stratification. ASATi, VATi, MFI and liver PDFF were log-transformed following testing residual vs. fitted value and due to Q-Q plots showing a significant (and consistent) departure from normality (Lumley et al., 2002). Since the heel BMD T-score measure already accounts for sex, this was not included in the respective model. FDR correction across all tests (N = 32) was used to control for multiple testing at alpha = 0.05, as it does not assume independence (Yekutieli and Benjamini, 1999). Sex was determined from genotyping analysis (UKB data field 22001). Since PGSs calculated from any set of SNPs will have a non-zero effect (Yang and Zhou, 2020), we compared the effect sizes of the various PGS oxt models against a distribution of effect sizes from models using PGSs calculated from 100 random and unique gene-sets of equal size to evaluate the specificity of the results. If the effect sizes from the PGS oxt models fell in the lower or upper 5% of the effect size distribution, it was considered to be specific to the oxytocin pathway.

Tissue enrichment in differentially expressed oxytocin pathway genes
Analysis of the GTEx v8 30-tissue dataset revealed statistically significant enrichment of differentially expressed genes (DEG) belonging to the oxytocin pathway in muscular tissue, brain tissue, bladder, breast, kidney, pancreas, oesophagus (Bonferroni corrected p < .001), heart, salivary gland (Bonferroni corrected p < .01), and skin (Bonferroni corrected p < .05) (Fig. 1, Supplementary Table 3). Analysis of the 54tissue dataset from the same donors adds some nuance to these results. The difference between subcutaneous adipose tissue and visceral adipose tissue is particularly notable (Fig. 1, Supplementary Table 4).

Pairwise phenotypic correlations and clustering
Pairwise correlation analyses and complete linkage hierarchical clustering revealed weak associations within four main variable clusters: one including HDL and apolipoprotein A1, one with cholesterol, LDL, and apolipoprotein B, one with liver PDFF, VAT, triglycerides, and glycated haemoglobin, and a fourth with ASAT, abdominal fat ratio, weight to muscle ratio, and muscle fat infiltration (Fig. 2, Supplementary Table 7).

Linear models
Comparing the effect sizes of the BMI PGS oxt models (obtained from linear models as standardised estimates) to the random gene-set model distributions, the effect sizes for the models of ASATi (p fdr =0.04), lipoprotein(a) (p fdr <0.001), HDL (p fdr =0.008) and heel BMD T-score (p fdr <0.001) were among the top 5% of their respective distributions, and that for the model of triglycerides (p fdr <0.001) the bottom 5% (Fig. 3). For the loneliness PGS oxt models (Fig. 4), the apolipoprotein A1 model had a significantly (p fdr =0.046, Supplementary Table 6) negative. association that fell in the lower 5% of the distribution. Results from these five BMI PGS oxt models and one loneliness PGS oxt model survived FDR correction ( Fig. 5; Supplementary Table 5). Associations between BMI PGS oxt and apolipoprotein B, total cholesterol and glycated haemoglobin, and between loneliness PGS oxt and HDL were also statistically significant but showed effect sizes falling outside the top or bottom 5% of their respective random gene-set model distribution. Therefore, these associations could not be confidently considered specific to the oxytocin pathway. The rest of these models were neither statistically significant nor specific to the oxytocin signalling pathway.

Discussion
The present results, derived from a large sample of adults, provide evidence for an association between loneliness and BMI polygenic scores computed exclusively for oxytocin signalling pathway genes, and a series of traits related to lipid homoeostasis and body composition. Importantly, permutation analyses demonstrated specificity for many of these associations to the oxytocin signalling pathway. The associations between PGS oxt for BMI and abdominal subcutaneous adipose tissue and blood lipids (e.g., HDL and triglycerides) are congruous with the suggested link between oxytocin and energy homoeostasis, body fat  metabolism, and body fat distribution (Leng and Sabatier, 2017;Winterton et al., 2020). An increased genetic load for higher BMI has been associated with increased subcutaneous fat, which unlike visceral adipose tissue, is not directly associated with insulin resistance (Alvehus et al., 2010) and elevated triglycerides (Porter et al., 2009), two components of metabolic syndrome. Subcutaneous adipose tissue is a fat depot specialised to provide long-term nutrient storage (Tchkonia et al., 2013), acting as a sink to sequester fatty acids. When subcutaneous fat is removed in mice, there is an increase of visceral fat mass, insulin resistance, circulating insulin, and TNF-α; while reimplantation of subcutaneous fat reverses these effects (Ishikawa et al., 2006). This is in line with the interpretation of different adipose depots as metabolically different endocrine organs (Tchkonia et al., 2013). This finding is supported by the results of our tissue expression enrichment analysis, which demonstrates an enrichment of oxytocin pathway genes in subcutaneous adipose tissue (as opposed to visceral fat) and is consistent with previous findings in mice where oxytocin administration reduced visceral and liver fat (Maejima et al., 2011). The observed association between BMI PGS oxt and abdominal subcutaneous fat is also in line with the associations between HDL and BMI PGS oxt , with HDL considered protective for cardiovascular events (Ali et al., 2012). In addition, the BMI PGS oxt was inversely associated with triglycerides, which is in agreement with elevation in HDL and the metabolic role of subcutaneous fat. The link between oxytocin signalling and BMD is also consistent with the evolutionarily ancient role of osteocytic bone tissue in energy metabolism (Haridy et al., 2021). This relationship is maintained in humans and appears to be reciprocal (Dirckx et al., 2019;Lee and Karsenty, 2008), which points to common regulatory mechanisms, such as the oxytocin system.
The negative association between the loneliness PGS oxt and apolipoprotein A1, which exerts a protective effect on cardiovascular disease, provides further evidence for pleiotropic effects of oxytocin gene variants on both energy homoeostasis and social behaviour. In other words, genetic variants in the oxytocin signalling pathway that are associated with loneliness are also linked with energy homoeostasis. Given oxytocin's secretion into both the central nervous system and the periphery, the oxytocin system is well-placed to coordinate both behaviour and energy regulation (Quintana and Guastella, 2020). Loneliness is also linked to reduced physical activity (Hawkley et al., 2009;Philip et al., 2020), which in turn negatively affects BMD (Benedetti et al., 2018) and lipid homoeostasis (Association, 2003;Delavar et al., 2011). Altogether, these results suggest a potential mechanism for how the oxytocin system integrates physiological and behavioural processes to regulate energy requirements (Quintana and Guastella, 2020). Moreover, the results highlight the treatment potential for stimulating the oxytocin system via intranasal administration.
This study had some limitations worth noting. First, it is conceivable that part of the observed genetic pleiotropy could be explained by Fig. 3. Effect size distributions for the BMI PGS oxt models (red confidence intervals) and random gene-set PGS models (grey confidence intervals). Highlighted in red for each distribution are the upper and lower 5 % to identify the effects specific to the oxytocin pathway. HDL -High Density Lipoprotein, LDL -Low Density Lipoprotein, BMD -Bone Mineral Density. Fig. 4. Effect size distributions for the Loneliness PGS oxt models (blue confidence intervals) and random gene-set PGS models (grey confidence intervals). Highlighted in blue for each distribution are the upper and lower 5 % to identify the effects specific to the oxytocin pathway. HDL -High Density Lipoprotein, LDL -Low Density Lipoprotein, BMD -Bone Mineral Density. population stratification and non-random sample selection in the original GWAS' (e.g., due to varying exclusion criteria related to somatic conditions for patients and healthy controls or general phenotypic correlations in the population that may or may not be due to overlapping mechanisms). We addressed this potential issue by using a careful quality control procedure for the genetic analysis and controlling for population stratification in our models. Second, a sex-disaggregated PGS score might have accounted for sex-specific confounders (Flynn et al., 2020), which is relevant for oxytocin research (Winterton et al., 2020). To address this, we adjusted for sex in our models and used T-scores for our BMD measures. Third, polygenic approaches cannot explain the biological mechanisms that underlie the relationships found, as PGSs assume an additive effect of individual alleles and do not model higher order relationships between variants. Fourth, the various phenotypes have unequal sample-sizes: comparatively few participants of UK Biobank had body MRI measures compared to the other measures, leading to decreased statistical power, which is of particular importance in studies seeking to find small effect sizes. On the other hand, the large sample size for the other phenotypes is one of the strengths of the study, along with using standardised methods and the inclusion of novel body composition measures.
Altogether, these results provide evidence for the involvement of the oxytocin pathway in energy metabolism, lipid homoeostasis and storage, bone density, and mental states such as loneliness in humans, hinting at complex pleiotropy and interconnected effects. These findings are in line with the notion that oxytocin plays a role in both social behaviour and metabolic traits. Moreover, this research provides a starting point for further research to elucidate these mechanisms and potentially find therapeutic applications of oxytocin in a broad range of conditions characterised by social and metabolic dysfunction, such as schizophrenia and autism spectrum disorders.

Conflicts of Interest
OAA has received speaker's honorarium from Lundbeck and is a consultant to HealthLytix. The remaining authors have no conflicts of interest to disclose.

Data Availability
The UK Biobank resource is open for eligible researchers upon application (http://www.ukbiobank.ac.uk/register-apply/). The code for the statistical analysis is available at https://osf.io/hy5aw/. #2019101, #2022080), and EU H2020 #847776 CoMorMent, the European Research Council under the European Union's Horizon 2020 research and Innovation program (Grant 802998). This research has been conducted using data from the UK Biobank Resources (www.ukb iobank.ac.uk), a major biomedical database, under Application no. 27412. This work was performed on Services for sensitive data (TSD), University of Oslo, Norway, with resources provided by UNINETT Sigma2 -the National Infrastructure for High Performance Computing and Data Storage in Norway.

Appendix A. Supporting information
Supplementary data associated with this article can be found in the Fig. 5. Summary of the estimates of the models that showed significant associations specific to PGSs calculated from oxytocin pathway genes. All these models are statistically significant after correcting for multiple tests (FDR). ASATi -Abdominal Subcutaneous Adipose Tissue volume, HDL -High Density Lipoprotein. online version at doi:10.1016/j.psyneuen.2022.105875.