A lead candidate functional single nucleotide polymorphism within the WARS2 gene associated with waist-hip-ratio does not alter RNA stability.

We have prioritised a single nucleotide polymorphism (SNP) rs2645294 as one candidate functional SNP in the TBX15-WARS2 waist-hip-ratio locus using posterior probability analysis. This SNP is located in the 3' untranslated region of the WARS2 (tryptophanyl tRNA synthetase 2, mitochondrial) gene with which it has an expression quantitative trait in subcutaneous white adipose tissue. We show that transcripts of the WARS2 gene in a human white adipose cell line, heterozygous for the rs2645294 SNP, showed allelic imbalance. We tested whether the rs2645294 SNP altered WARS2 RNA stability using three different methods: actinomycin-D inhibition and RNA decay, mature and nascent RNA analysis and luciferase reporter assays. We found no evidence of a difference in RNA stability between the rs2645294 alleles indicating that the allelic expression imbalance was likely due to transcriptional regulation.

Introduction

Body fat distribution is a risk factor for disease, independent of obesity. It can be easily assessed by waist-hip ratio (WHR), the ratio of waist circumference to hip circumference. Higher WHR, suggesting greater visceral fat accumulation, is associated with increased mortality and risk of coronary heart disease, myocardial infarction and type 2 diabetes (T2D) [[1], [2], [3], [4], [5]]. Furthermore, WHR has been established as a causal risk factor for T2D and CVD risk by Mendelian Randomization analyses [[6], [7], [8], [9]]. Twin studies have estimated that the variation in WHR is 31–61% heritable and the most recent genome-wide association study (GWAS) has discovered 346 different genetic loci associated with WHR adjusted for body mass index (WHRadjBMI) [[10], [11], [12], [13]]. The majority of these loci are located in non-coding regions, which complicates the identification of target genes and mechanisms underlying fat distribution [14,15].

The TBX15-WARS2 locus on chromosome 1 has been consistently associated with WHRadjBMI across multiple meta-analyses and up to four potentially independent association signals were discovered in the region, suggesting multiple functional SNPs and genes may be involved [13,[16], [17], [18], [19]]. The WHR-association signal spans ~1 Mb and includes genes T-box 15 (TBX15), mitochondrial tryptophanyl (W) tRNA synthetase 2 (WARS2), and regions downstream of SPAG17. The expression of both TBX15 and WARS2 has been associated with metabolic traits in humans [20].

Several studies have linked the TBX15-WARS2 risk SNPs to TBX15 expression in adipose tissue. However, data from the GTEx database links the locus predominantly to reduced expression of WARS2 across multiple tissues [16,20,21]. WARS2 is a nuclear-encoded mitochondrial tryptophanyl-tRNA synthetase, recently associated with angiogenesis, adiposity and brown adipose tissue metabolism [[22], [23], [24]]. Variants in another nuclear-encoded mitochondrial tRNA synthetase, aspartyl-tRNA synthetase (DARS2) has been associated with WHR and 79 additional mitochondrial - or mitochondrial nuclear – encoded variants associated with adipose measurements, pointing to genes involved in mitochondrial function as strong candidates [18,25]. These findings make WARS2 a good potential functional WHR gene within the locus.

Shungin et al. reported four index SNPs near TBX15-WARS2-SPAG17 each denoted as defining regions D, E, F and G within the locus (Supplementary Fig. 1). These contain correlated but different lead SNPs near WARS2 in men (region F, rs1106529) and women/sex combined (region G, rs2645294), an independent sex-combined signal in TBX15 (region E, rs12143789) and an independent sex-combined signal near SPAG17 (region D, rs12731372) [17].

We have performed Posterior Probability Analysis (PPA) of chromosome 1 SNPs for WHR in UK Biobank and GIANT datasets, to identify a total of 5 SNPs with PPA scores >0.2, which were localised in region G and D (3 and 2 SNPs respectively, including both the index SNPs). In this manuscript we consider the region G (WARS2) SNPs. Further, as rs2645294 was within the 3′UTR of WARS2 we hypothesised that it influenced RNA stability and post-transcriptional regulation of WARS2. To test this, we analysed the WARS2 3′UTR, used nascent RNA qPCR, actinomycin D-mediated transcription inhibition and 3′UTR luciferase assays in differentiating human white adipocyte cells. We found allelic imbalance in the transcription of WARS2, which we believe to be due to transcriptional regulation rather than due to alterations in RNA stability.

Materials and methods

Bioinformatic prioritisation of SNPs

We prioritised SNPs using posterior probability analysis (PPA) of each SNP in the TBX15-WARS2 locus, as described previously [26]. The analysis was performed independently for UK Biobank (UKBB) and The Genetic Investigation of Anthropometric Traits (GIANT) Consortium data in male-only, female-only or sex-combined groups. All SNPs with PPA scores ≥ 0.2 in any of the sub-groups were chosen for further analysis.

Bioinformatic analysis of risk block G SNPs

To assess altered transcription-factor binding sites (TFBSs), we interrogated risk block G SNPs using the Haploreg v4.1 web server with default settings [27]. We then tested all the SNPs in the TBX15-WARS2 locus, defined by Shungin et al. [17], for the likelihood of effect on transcriptional regulation using Phylogenetic Module Complexity Analysis (PMCA), exactly as described previously [28].

Cell culture and differentiation

HEK293T cells were purchased from Public Health England ECACC General Cell Collection and used between passages P14 and P17 (ECACC 12022001, lot 16G020, Sigma-Aldrich). The hWAT cells were provided by Prof Yu-Hua Tseng, PhD at Harvard Medical School, Joslin Diabetes Center, Boston, MA 02215. The cells were previously isolated from the white fat in the neck of a female, aged 56 with a BMI of 30.8, and immortalised using telomerase reverse transcriptase (TERT)-expressing plasmid [29]. hWAT cells were used between passages P17 and P23. Both cell lines were cultured at 37 °C, 5% CO2 in DMEM GlutaMAX (# 10569010 DMEM, high glucose, GlutaMAX™ supplement, pyruvate, Thermo Fisher) with 10% fetal bovine serum (FBS, Thermo Fisher) and 1% Penicillin-Streptomycin (P/S, Thermo Fisher). For hWAT differentiation in T75 flasks, 2-days post-confluent cells were induced using a freshly prepared DMEM medium with 10% FBS, 1% P/S, 0.5 μM insulin, 500 μM isobutyl methylxanthine (IBMX), 2 nM 3,3′,5-triiodo-L-thyronine (T3), 0.1 μM dexamethasone, 17 μM panthoneate, 30 μM indomethacin, 33 μM biotin (all from Sigma-Aldrich). The differentiation media was replaced after three days and on day 4 cells were either harvested for nuclear protein isolation or transfected with luciferase plasmids.

DNA isolation and sequencing

The genomic DNA of hWAT, HEK293T and Simpson Golabi Behmel Syndrome (SGBS) cells was isolated from ~3 million cells using DNeasy Blood and Tissue Kit (Quiagen), according to manufacturer′s protocol. The 379 bp sequence surrounding rs2645294 was amplified using the primers rs264-FW IV (5′-ATGTGACCACGGTTCTGTGA-3′), rs264-RV III/IV (5′-AAGAGCCCAAGTCCCTGAAT-3′) and the Q5® High-Fidelity DNA Polymerase (NEB). The PCR product was then sequenced using the same primers and the Source BioScience (Nottingham, UK) Sanger sequencing service.

Adipocyte nuclear protein extraction

Nuclear protein was extracted from Day 4 differentiated hWAT adipocytes using the NE-PER™ Nuclear and Cytoplasmic Extraction Reagents (Thermo Fisher). The quantity of protein was measured with the DC™ Protein Assay (Biorad).

EMSA

The EMSA was performed according to the LightShift Chemiluminescent EMSA Kit protocol (ThermoFisher). Probes of 41 base-pairs were designed for each SNP, this was selected to provide a target longer than a typical transcription factor nominal binding site without introducing more binding sites unrelated to the SNP site, giving larger mobility shifts than longer fragments, with the caveat that this was at the expense of larger complex binding. Briefly, the 3′ biotinylated double stranded probes (20fmol) of either allele were incubated with day 4-differentiated hWAT nuclear protein (4.5 μg) in final concentration of 10 mM Tris, 50 mM KCl, 1 mM DTT and 50 ng/μl Poly (dI·dC) (ThermoFisher). The reactions were incubated at room temperature for 20 min and resolved on a 6% DNA retardation gel (Thermo Fisher) in 0.5 M TBE buffer. The gel contents were transferred onto a nylon membrane, UV-crosslinked, then blocked and visualised using the Chemiluminescent Nucleic Acid Detection Module (Thermo Fisher). The forward biotinylated sequences are shown in (Supplementary Table 1).

RNA isolation and cDNA synthesis

Total RNA was extracted using the TRIzol reagent (Thermo Fisher) according to the manufacturer's instructions, including the DNase I treatment step. 2 μg of RNA was reverse transcribed to cDNA in a 20 μl reaction using the SuperScript™ III Reverse Transcriptase Kit (Thermo Fisher). Quantitative PCR (qPCR) was carried out with the TaqMan™ Universal PCR Master Mix (Thermo Fisher) and FAM-labelled probes (ThermoFisher) on the ABIPRISM 7500 Fast Real-Time PCR System (Applied Biosystems) using the ‘Fast’ protocol with 40 cycles and quantitation – comparative CT (∆∆CT) analysis.

Allele-specific qPCR

The relative expression of the two rs2645294 alleles in hWAT cells was assessed using a TaqMAN SNP Genotyping Assay (Assay ID: C_15913675_20, Thermo Fisher), which contains VIC or FAM labelled allele-specific probes, following the standard TaqMAN Assay protocol. The Ct values of each allele were normalised by β-actin expression fold change ratio of C:T determined using the Comparative CT Method (ΔΔCT Method).

Violin eQTL plots were downloaded from GTEx Portal (https://www.gtexportal.org/home/). The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from: the GTEx Portal on 10/12/19.

Transcription inhibition

hWAT cells were plated at 120K/well (Scepter™ 2.0 Cell Counter, gating: 11 μm–24 μm) in 6-well plates to reach 60% confluence next day. After 24 h, cells were treated either with Actinomycin D (Sigma Aldrich, final concentration: 2 μg/ml) or corresponding volume of DMSO for 0, 4, 8, 12 and 24 h. At each time-point, cells were lysed and frozen in TRIzol (Thermo Fisher). Expression levels were assessed using the WARS2 genotyping probe (above) and MYC TaqMAN probe (Hs00153408_m1, Thermo Fisher). House-keeping genes were assessed using Double-dye Hydrolysis geNorm probes (PrimerDesign) for beta actin (ACTB), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and 18S rRNA (18S). Transcript half-life was determined using t1/2 = ln(2) / −slope, where the slope was calculated from a plot of the log(2) relative expression (WARS2 or MYC transcript normalised to ACTB) using the Comparative CT Method (ΔΔCT Method).

Bioinformatic analysis of WARS2 3′UTR

To extract the WARS2 3′UTR sequence, we used the ENST00000369426.9 transcript which has the longest (2144 bp) 3′UTR sequence of all the WARS2 isoforms annotated in GENCODE v24. We analysed this sequence using RegRNA 2.0 using default settings for all RNA motifs, except TRANSAFAC TFBSs (search all motifs and sites only for ‘Homo sapiens’, long stems ≥ 40 bp, functional RNA sequences of similarity ≥ 0.9 or match_length ≥ 30 bp, miRNA target sites with score ≥ 170 & free_energy ≤ −25, non-coding RNA-hybridization sites with length ≥ 20 & free_energy ≤ −20, GC-content ratio for 100 bp window size, RNA accessibility with maximum pair distance of 100 bp and consecutive unpair size of 6 bp, open reading frames (ORFs) predicted with “AUG” start codon only) [30]. To test the effect of rs2645294 on RNA-binding protein (RBP) binding, we used ‘a database of RNA-binding proteins and associated motifs’ AtTRACT to scan a sequence 20 bp surrounding the SNP for motifs of length between 4 and 8 bases [31]. For the protein motifs that were specific only to one or the other allele, the position-weight matrices (PWMs) were downloaded from AtTRACT and their occurrence probability further quantitatively assessed for each allelic sequence by Find Individual Motif Occurrences (FIMO, http://meme-suite.org/tools/fimo) with a cut-off p-value < 0.01 [32]. For each protein, only the motif with lowest p-value was shown. To assess the presence of known miRNA-binding sites, we searched 100 bp in the proximity of the SNP for known motifs using the ‘Scan for Motif’ tool in Transterm and selecting to ‘Show targets of conserved microRNA families as predicted by Targetscan’ [33,34]. To study the effect on secondary RNA structure, we tested each allelic 3′UTR sequence (full 2144 bp) using the “Predict a Secondary Structure Web Server” function of the RNA Structure web server with default settings (temperature 310.15 K, maximum loop size: 30 bases, maximum % energy difference: 10, maximum number of structures: 20, window size: 3, gamma (MEA): 1, iterations (pseudoknot prediction): 1, minimum helix length (pseudoknot prediction): 3. We show the most energetically-favourable structure for each of the algorithms – Fold, MaxExpect, ProbKnot) [35].

Nascent RNA isolation & validation

PolyA− and PolyA+ RNA were isolated using the PolyATtract® mRNA Isolation Systems (Promega). To ensure improved clearance of polyadenylated mRNA, the PolyA− fraction was passed through a fresh PolyA+ binding column for a second time. WARS2 allelic levels were assessed by allele-specific qPCR described above. The method was validated by enrichment of spliced over unspliced WARS2 in the polyA+ mRNA fraction, assessed by specifically designed qPCR primers (LGC Biosearch) targeted to regions between exons 2 and 3 of WARS2, as listed in Supplementary Table 2 and depicted in Fig. 6A.

Fig. 6

ASE is present both at the mature and nascent RNA level in hWAT cells. Diagram illustrating putative RNA processing of the WARS2 transcript and the location of primers for both the spliced and unspliced WARS2 transcripts (A). RNA from hWAT cells was separated into polyA+ and PolyA− fractions and assessed by qPCR. The ratio of spliced to unspliced WARS2, calculated by 2^(Ct(unspliced) − Ct(spliced)) without further normalisation. The polyA+ fraction showed a clear enrichment of spliced WARS2 RNA (B). Analysis of the C to T ratio ((2^−(deltaCt)) results showed there was no statistically significant difference between any of the RNA fractions (C). The raw Ct values used to calculate the spliced to unspliced WARS2 ratios in (B) and the C:T ratio in (C) are shown in the table (D). Comparison is by one-way ANOVA. **** for Total cDNA vs poly A+, p ≤ 0.0001; #### for PolyA− vs poly A+, p ≤ 0.0001.

3′UTR luciferase assays

The 3′UTR of WARS2 was amplified from hWAT cell genomic DNA using Q5® High-Fidelity DNA Polymerase (NEB) and cloned using SalI-HF and XhoI downstream of the luciferase gene in the pmirGLO Dual-Luciferase miRNA Target Expression Vector (Promega, referred to as empty vector). Cloning was verified by Sanger sequencing. The 3′UTR sequence was chosen from ENST00000369426.9, the WARS2 transcript with the longest annotated 3′UTR in GENCODE v24. The insertion was confirmed by sequencing and restriction enzymes digest. The plasmids were then mutated to generate the alternative alleles using the Q5® Site-Directed Mutagenesis Kit (NEB) and the inserts sequenced to ensure no additional mutations were present. The primers used are listed in Supplementary Table 3. On the day 4 of differentiation, 4000 hWAT cells/well were plated into a white solid bottom 96-well Greiner Bio-One CELLSTAR plate. 8 technical replicates (wells) per plasmid were included. 24 h later, the cells were transfected with 0.2 μg of pmirGLO with or without the 3′UTR sequence and 0.2 μl Lipofectamine® 3000 per well according to the commercial protocol. Each experiment also included pEGFP_C1 (BD Biosciences Clontech, cat. no. 6084-1) as a transfection control. If ~10% transfection was observed, the cells were lysed and the luciferase activity assayed using the Dual Luciferase Assay Kit (Promega) and the Varioskan® Flash microplate reader (ThermoFisher). Fold change luciferase values were relative to empty vector which was the pMIR-GLO vector that contained a moderate-strength PGK promoter driving the firefly luciferase luc2 reporter without the 3′UTR test sequence, in addition to the SV40 early enhancer/promoter driven Renilla luciferase (hRluc-neo) used as the control reporter for normalisation.

Statistical analysis

Statistical analyses were conducted with GraphPad PRISM 6 software package. The D'Agostino-Pearson Omibus test was used to evaluate data normality and the appropriate parametric or non-parametric tests were used. Unpaired two-tailed Student's t-tests, one-way ANOVA or Kruskal-Wallis tests were applied to compare two or multiple groups in the qPCR and luciferase assay data. Rates of RNA degradation were calculated using a linear regression, as described above.

Results

Prioritising potentially functional SNPs

The TBX15-WARS2 locus contains four independent haplotype regions defined by different lead SNPs (regions D, E, F and G) associated with WHRadjBMI [17]. Sixty-two SNPs were identified in close LD (r2 ≥ 0.8), with the four lead SNPs for all regions stretching between WARS2, TBX15 and SPAG17. Except for rs10494217 (TBX15-H156N) in region E, all SNPs were non-coding [17].

Using posterior probability (PP) analysis on the UK Biobank and GIANT GWAS datasets we assigned a posterior probability for a causative effect to each SNP in the locus [26]. Selecting a cut-off PP of ≥0.2 we identified three SNPs in region G (WARS2) and two in region D (SPAG17) (Table 1). In this study, we further investigate the three region G SNPs as candidates for likely causal SNP(s) in the region. The rs2645294 index SNP showed the highest PPA value, but only in females in the GIANT cohort.

Table 1

Prioritising SNPs within risk block G using posterior probability. The three SNPs having PP > 0.2 in either GIANT of UKBB are shown with the rounded results for female, male and combined (All) in each dataset. Base pair position is shown in hg38.

SNP	Chromosome	Base pair hg38	PPA score
			UK BioBank			GIANT
			All	Female	Male	All	Female	Male
Region G, index SNP rs2645294 (WARS2)
rs2645294	1	119,031,964	0.05	0.00	0.01	0.33	0.60	0.01
rs10923724	1	119,004,219	0.53	0.00	0.07	0.29	0.16	0.01
rs6428789	1	119,007,857	0.27	0.00	0.05	N/A	N/A	N/A
Region D, index SNP rs12731372 (SPAG17)
rs12731372	1		0.01	0.00	0.01	0.68	0.56	0.43
rs7534091	1		0.00	0.00	0.05	0.25	0.13	0.33

PP values ≥ 0.2 are shown in bold text.

Initial analysis using HaploREg indicated that the region G SNPs in linkage disequilibrium (r2 ≥ 0.8) with the rs2645294 index SNP nominally have the potential to alter predicted transcription factor binding sites (Supplementary Table 4). We then further assessed these SNPs using Phylogenetic Module Complexity Analysis (PMCA), a method that uses conservation of sequence, order and distance of TFBS motifs, in conjunction with convolutional neural networks (CNN) that predict the regulatory activity of a given variant. This revealed low ranks for the three region G PPA shortlisted variants, indicating a lack of regulatory activity for those variants (rs2645294, 50th rank; rs10923724, 21st; rs6428789, 53rd; the two other region G SNPs rs7553422 and rs1886914 were ranked 30th and 39th - data not shown) [28,36,37]. Furthermore, when overlaid with human white adipocyte cell (hWAT [29]) ATAC-sequence data (Sinnott-Armstrong et al. accepted for publication), none of the three SNPs overlapped open chromatin in preadipocytes (Day zero, D0) or at any stage of differentiation (D3, 6 or 24 of adipocyte differentiation) (Fig. 1).

Fig. 1

Risk block G SNPs do not overlap open chromatin in hWAT cells. Alignment of the 3 short-listed SNPs in risk block G with the ATAC-Seq signal of differentiating hWAT preadipocytes (Day 0, Day 3, Day 6, Day 24; NCBI Sequence Read Archive (SRA) PRJNA664585) and the primary state ChromHMM annotation from Adipose Nuclei (FAT.ADIP.NUC), Mesenchymal Stem Cell Derived Adipocyte Cultured Cells (FAT.MSC.DR.ADIP) and Adipose Derived Mesenchymal Stem Cell Cultured Cells (FAT.ADIP.DR.MSC) [38] (Sinnott-Armstrong et al. accepted for publication). Region, GRCh37/hg19, chr1:119,531,043-119,576,244, visualised using UCSC genome browser.

Next, we tested the three SNPs using an electrophoretic mobility shift assay (EMSA) with nuclear protein from early differentiating hWAT cells (Day 4) and found no effect of these SNPs on differential transcription factor binding (Fig. 2 & Supplementary Fig. 2).

Fig. 2

Prioritised risk block G SNPs do not affect protein binding in differentiating hWAT cells. Double stranded biotin-labelled DNA probes, 38 bp in length surrounding each SNP and carrying either genotype were incubated with nuclear protein from day 4 differentiated hWAT cells and resolved on a 6% DNA-retardation gel. No reproducible difference, within the limited sensitivity of these assays, was observed in the protein bound bands between the alleles in four replicates (see Supplementary Fig. 2).

In summary, we found no evidence for an effect of the three prioritised SNPs on cis-regulatory element (CRE) activity or protein binding. We therefore considered alternative mechanisms of gene regulation in the risk block. The highest-ranking SNP in the PP analysis was rs2645294 with a PPA score of 0.6 in females of the GIANT Consortium (Table 1). Interestingly, this SNP overlaps a 3′ untranslated region (3′UTR) of the WARS2 gene that was annotated with transcribed histone marks in adipose tissue chromatin (Fig. 1). There was no indication of enhancer mark annotation (Fig. 1). Since rs2645294 is an eQTL for WARS2 in multiple tissues, including subcutaneous and visceral adipose tissue, we hypothesised that the rs2645294 SNP could act by altering WARS2 RNA stability, thus leading to reduced WARS2 levels and potentially contributing to the WHRadjBMI association through mitochondrial metabolism alterations.

Bioinformatic analysis of WARS2 3′UTR RNA regulatory elements

Elements within the 3′UTR of genes can recruit miRNAs, non-coding RNA or RNA-binding proteins (RBPs) thus affecting stability, translation or localisation of the RNA [39,40].

We used RegRNA 2.0 to scan the 2144 bp 3′UTR of WARS2, but found no direct overlap of rs2645294 with regulatory elements, such as splicing sites, splicing regulators, polyadenylation sites, structural sequences and miRNA binding sites (Supplementary Fig. 3) [30]. The lack of overlap with miRNA binding sites was confirmed by TargetScan and Transterm, although Transterm revealed that two miRNAs could potentially bind nearby to the SNP (Supplementary Fig. 4) [33,34].

To predict which RNA binding motifs could be affected by rs2645294, we used the dAtabase of RNA binding protein and AssoCiated moTifs (AtTRACT), a comprehensive database with information on 370 RBPs and 1583 RBP consensus binding motifs [31]. The G allele (RNA) created additional RNA binding sites for UGBP Elav-Like Family Member 1–4 (CELF1–4), Eukaryotic translation initiation factor 4B (EIF4B) and Fused in Sarcoma/Translocated in Sarcoma (FUS) (Table 2).

Table 2

The RNA binding motifs introduced by the G allele of rs2645294. The 100 bp surrounding rs2645294 with either allele were scanned using AtTRACT (https://attract.cnic.es/). Six novel motifs were introduced by the G allele. Since WARS2 is transcribed from the minus strand, SNP alleles C and T in the DNA are referred to as G and A, respectively, on the single stranded RNA level.

Protein	Gene ID	Motif	Experiment
CELF1	ENSG00000149187	UGUU	X-ray diffraction
CELF2	ENSG00000048740	UGUU	X-ray diffraction
CELF2	ENSG00000048740	UGUUG	SELEX
CELF4	ENSG00000101489	GGUGUUG	RNAcompete
EIF4B	ENSG00000063046	GUUGGAA	SELEX
FUS	ENSG00000089280	GGGUGU	SELEX. SDS-PAGE, EMSA, UV crosslink, competition and immunoprecipitation assays

To gain a more quantitative assessment of these protein binding RNA motifs we then ran a Find Individual Motif Occurrences (FIMO) analysis on the five proteins found by AtTRACT [32]. These analyses indicate the strongest likelihood for binding EIF4B (q-value ~0.002) and weaker support for the G allele binding FUS (q-value ~0.04) and CELF4 (q-value ~0.03) (Table 3). These analyses support the potential for altered regulation of the WARS2 transcript through alteration of protein binding.

Table 3

The RNA binding motifs introduced by the G allele of rs2645294. The 20 bp surrounding rs2645294 with either allele were scanned using AtTRACT (https://attract.cnic.es/) discovering 5 novel motifs introduced by the G allele. Their matching sequences and experiments that derived their respective position-weight-matrices (PWM) are listed. The probability of a random sequence occurring and matching the same PWM was obtained by Find Individual Motif Occurrences (FIMO, http://meme-suite.org/tools/fimo) with a p = 0.01 cut-off. For each protein, only the motif with lowest p-value is shown. Since WARS2 is transcribed from the minus strand, SNP alleles C and T in the DNA are referred to as G and A, respectively, on the single stranded RNA level.

Protein	Motif	Experiment	FIMO p-value		FIMO q-value
			C	G	C	G
EIF4B	GUUGGAA	SELEX	0.0009	6E−05	0.0128	0.0018
FUS	GGGUGU	SELEX	0.0186	0.0012	0.119	0.0365
CELF4	GGUGUUG	RNAcompete	0.0258	0.0019	0.129	0.0284
CELF1	UGUU	X-ray diffraction	0.0508	0.0039	0.914	0.141
CELF2	UGUU	X-ray diffraction	0.0508	0.0039	0.914	0.141

In addition, a change in RNA structure due to the SNP could indirectly alter accessibility of surrounding RBPs or miRNA binding motifs. Indeed, a major change in the stem loop architecture due to the SNP was predicted by Fold and ProbKnot algorithms of the RNA Structure Web server (Fig. 3) [35]. On the other hand, MaxExpect predicted RNA structures that were almost identical for the two alleles.

Fig. 3

The effect of rs2645294 on RNA structure predictions. Comparison of results of the 3 algorithms available at the RNA Structure server [35]. 100 bp sequence surrounding the SNP was used. Since WARS2 is transcribed from the minus strand, SNP alleles C and T in the DNA are referred to as G and A, respectively, on the single stranded RNA level. Structures with the lowest folding energies for Fold and MaxExpect and ProbKnot algorithms are shown, respectively. Probabilities of base positions in the structures are color-coded according to the legend.

In summary, the WARS2 3′UTR SNP does not directly overlap a miRNA binding sequence, but may directly affect RBP binding leading to a change in RNA stability. The SNP could also affect RNA structure and thus indirectly affect RBPs and miRNAs binding at surrounding sites.

Allele-specific expression (ASE) at rs2645294 for WARS2 gene expression, but no difference in RNA degradation

Sanger sequencing revealed that hWAT cells (a human white adipose cell line) are heterozygous for rs2645294 (Fig. 4A). In order to test for allele-specific expression (ASE) in this cell line we selected a TaqMAN allele-specific genotyping assay and confirmed its specificity using two cell lines one with a T/T and the other with C/C genotype (Supplementary Table 5). We then tested the hWAT cells and found that the T risk allele was expressed >2 fold lower than the C allele (p < 0.0001) in these cells which agrees with the direction in the GTEx database (Fig. 4B and C).

Fig. 4

hWAT cells are heterozygous for rs2645294 and show allele-specific expression of WARS2. (A) Sanger sequencing of rs2645294 in hWAT DNA revealed two peaks, C and T. (B) Allele-specific TaqMAN genotyping probes were used to analyse both genomic gDNA and cDNA from hWAT cells, showing that the C allele is expressed >2 fold higher than the T (WHRadjBMI-increasing) allele in the cDNA. Comparison is by unpaired t-test, **** for p < 0.0001. (C) The Genotype-Tissue Expression (GTEx) Project eQTL violin plot of WARS2 expression in subcutaneous adipose tissue for either the homozygous reference alleles (C/C), heterozygous (C/T) or homozygous for the alternative alleles (T/T) obtained from the GTEx Portal on 10/12/19 [41].

To compare the RNA stability of the two WARS2 alleles, we used actinomycin D, an RNA polymerase inhibitor, to inhibit transcription in hWAT cells. We first tested the effect of two different actinomycin D concentrations on the RNA stabilities of WARS2, a MYC positive control with a short half-life and housekeeping genes at different time-points. Degradation over time of both alleles of WARS2 was observed at both concentrations of actinomycin, compared to DMSO treated cells showing treatment was effective (unpublished data).

Since, 2 μg/ml actinomycin D treatment for 0–24 h did not interfere with expression of housekeeping genes, but rapidly reduced MYC gene expression as previously reported, we chose this condition for further experiments (unpublished data) [42]. We performed three independent replicates of actinomycin D inhibition in hWAT cells for 0–24 h. The half-life of WARS2 RNA was 6.3 h, but no differences between the degradation rate of the two alleles was observed (linear regression of log2(C) vs log2(T), p = 0.6550) showing that rs2645294 does not impact RNA stability of the WARS2 transcript (Fig. 5A and B). We confirmed rapid degradation of MYC indicating that the actinomycin D inhibition was successful (Supplementary Fig. 5).

Fig. 5

RNA stability analysis of the two WARS2 allele transcripts.

hWAT cells were treated with Actinomycin D (ActD, 2 μg/ml) or DMSO (control) for 0–24 h and allele-specific qPCR was used to assess differences in allele-specific RNA stability of WARS2. All data were normalised to highly stable β-actin and are shown as mean ± SD. (A) Log(2) normalised expression of WARS2 (n = 3 separate replicate experiments, each with 3 technical replicates, for each genotype/treatment). The linear regression slopes of T and C alleles upon actinomycin treatment were not significantly different. Half-life was calculated as ln(2) / −slope and was t1/2 = 6.3 h. (B) Representation of the same data as in A, showing linear regression slopes for the difference between the Ct values of the two alleles over time, which were not significantly different.

Allele-specific expression is present at both the mature and nascent RNA level in heterozygous adipocytes

We hypothesised that if the WARS2 adipocyte allele specific expression arose post-transcriptionally, it should be less pronounced at the nascent RNA level and only appear after RNA processing at the mature RNA level (Fig. 6A). It was previously shown that PolyA− RNA can be used as a nascent mRNA fraction [43]. We thus separated RNA from undifferentiated hWAT preadipocytes into polyA+ and PolyA− RNA. First, we used exon and intron-specific primers to show ~3-fold enrichment of mature spliced WARS2 RNA in the polyA+ fraction compared to total RNA and PolyA− RNA (p < 0.0001 for both comparisons) (Fig. 6B). We then used the rs2645294 allele-specific probes to assess ASE in the different fractions. The mean ∆CT between C and T alleles in total RNA, poly A+ and PolyA− RNA was −1.16, −1.29 and −1.45, respectively, without any statistically significant difference between the values (p = 0.1143) (Fig. 6C and D). Thus, rather than the expected decrease in the C:T ratio, we observed no difference between the C:T ratios for the PolyA− fraction compared to the poly A+ RNA fraction (Fig. 6C). Normalising to the spliced WARS2 probe (2^−(Ct(C or T) − Ct(Spliced WARS2))) indicates that the T and C allele account for 25% and 62% of the WARS2 mRNA respectively in these assays accounting for approximately 87% of the polyA+ mRNA (the other 13% may be accounted for by differences in probe efficiencies and their location in the transcript if there is degradation). This result supports the hypothesis that a transcriptional mechanism accounted for the observed ASE in hWAT cells rather than posttranscriptional regulation.

SNP rs2645294 does not alter the expression of a luciferase-WARS2-3′UTR reporter

Finally, to assess at the protein level whether the 3′UTR rs2645294 SNP in WARS2 altered transcript stability or translation we cloned the 3′UTR downstream of the luciferase gene in the pMIR-GLO vector. We used site-directed-mutagenesis to make constructs with either C or T alleles, and transfected these into hWAT or HEK293T cells and measured reporter luciferase activity (Fig. 7A). We used HEK293T cells in addition to hWAT cells in order to have a comparison biological replicate in an easily transfectable and widely used cell line, although these cells could lack the necessary transcription factors found in adipose cells, however the WARS2 eQTL is observed in many tissues in GTEx. The presence of the WARS2 3′UTR caused a ~3-fold reduction in reporter signal compared to the empty vector both in hWAT cells (C vs empty: p = 0.0056, T vs empty: p = 0.0007) and HEK293T cells (C vs empty: p = 0.0194, T vs empty: p = 0.0001). The reduction is independent of the allele inserted and reflects the lower expression rates often observed with longer 3′UTRs, presumably due to the presence of binding sites that inhibit translation [44]. In line with the computational predictions, no difference in luciferase activity between the two alleles was observed in either cell line, we did not find any evidence for rs2645294 regulating WARS2 protein or RNA levels (Fig. 7B and C).

Fig. 7

rs2645294 has no effect on expression of the luciferase fused to the 3′UTR of WARS2. (A) The 2144 bp 3′UTR sequence encompassing the rs2645294 SNP was cloned 3′ of the luciferase gene in the pMIR-GLO vector. (B and C) The vectors with either SNP (C or T) were then transfected into day 4-differentiated hWAT cells (B) or HEK293T cells (C) for 24 h. Fold change was relative to empty vector which is the pMIR-GLO vector containing a moderate-strength PGK promoter driving the firefly luciferase luc2 reporter without the 3′UTR test sequence. Data plotted as mean ± SD. Comparisons are by Kruskal–Wallis test. *p < 0.05 vs empty vector, **p < 0.01 vs empty vector, ***p < 0.001 vs empty vector.

Discussion

The TBX15-WARS2 locus contains multiple independent association signals for WHRadjBMI. Here, we focused on the Shungin et al. association region-G signal [17]. Posterior probability analysis of SNPs within this region indicated three potentially functional SNPs with rs2645294 attaining the highest score of 0.6 in females of the GIANT dataset. The 3 SNPs we focused on in our downstream analyses are those with the highest PPA in region-G, given the available data. PPA scores were markedly different between UKBB and GIANT datasets and between sexes and may be explained by the differing sample sizes between the two studies, the sexual dimorphism of fat distribution genetics, and by the different ancestral structures between the two meta-analyses [[45], [46], [47]]. The three SNPs shortlisted in this approach scored poorly in sequence-based computational analyses designed to identify regulatory SNPs, suggesting that the SNPs are unlikely to have a genomic regulatory function.

Using multiple online databases, we found that rs2645294 was not predicted to directly affect miRNA binding, but that the G allele, associated with higher WARS2 expression levels, introduces new binding motifs for six RNA Binding proteins (RBPs) with strongest support for EIF4B. Binding of some of these RBPs could lead to increased stability and levels of WARS2 RNA [40]. For example, the CELF family was shown to affect alternative splicing, C to U RNA editing, de-adenylation, mRNA decay, and translation of multiple transcripts [48]. In order to test the hypothesis that rs2645294 affects RNA stability and leads to allele specific differences in WARS2 RNA expression levels we used three different approaches.

Actinomycin D inhibition was previously used to detect a stabilising effect of a schizophrenia-linked variant on Kalirin (KALRN) mRNA [49]. Using the same method in a human white adipose cell line, we determined the half-life of WARS2 to be 6.3 h, similar to the 7.2 h reported for Wars2 in a transcriptome-wide study of mouse embryonic stem cells [42]. No difference in degradation of the two WARS2 alleles was detected after actinomycin D inhibition showing that the allele specific expression (ASE) in the hWAT cell line is not due to an allelic difference in RNA stability. However, it has been shown that transcription and RNA degradation are tightly linked and the effects of transcription inhibition are thus buffered by consequential downregulation of mRNA decay [50,51]. Indeed, using comparative dynamic transcriptome analysis, the median RNA half-life is shorter than reported by actinomycin D inhibition or pulse-chase [42,52]. To avoid this problem, we assessed the allelic expression at the level of nascent RNA. Previously, the PolyA− fraction was used to estimate the nascent RNA population [43]. We found no difference in ASE of WARS2 between the polyA+ and PolyA− fractions, further suggesting a transcriptional origin of the ASE. Although we observed an enrichment of spliced to unspliced WARS2 in the polyA+ fraction, no reduction of the ratio compared to the total RNA was observed in the PolyA− fraction. This could be explained by the co-transcriptional nature of splicing [53]. Since our qPCR probes for unspliced WARS2 targeted intron 2 and exon 2 boundary, it is possible that the majority of these early introns were already spliced before polyadenylation.

Finally, both the actinomycin experiment and nascent RNA qPCR were carried out using a human WAT cell line where other SNPs differing between the two WARS2 alleles could have concealed an effect of rs2645294 on RNA stability. Thus, to isolate the effect of the single SNP, we cloned the 3′UTR of WARS2 in a luciferase vector and found no differences between the expression of the two alleles. To prevent potential artefacts such as the effect of the luciferase gene on the 3′UTR structure, ideally, rs2645294-directed mutagenesis of endogenous WARS2 in hWAT cells would be performed in the future.

In conclusion, although each method had its limitations, none of the three different approaches showed any evidence for rs2645294 affecting RNA stability. The effect of region G on WARS2 RNA levels is therefore likely regulated through a transcriptional mechanism involving other SNPs in the risk block, although alterations in nuclear export cannot be ruled out. For example, rs6428789 which overlaps a putative bivalent enhancer in adipose nuclei, or rs10923724, with a high UK Biobank PPA and the highest PMCA and CNN-based scores in region G. As outlined in the introduction WARS2 function has been linked with adiposity and thus alteration in its expression could have a physiological effect leading to altered fat distribution [23,24].

CRediT authorship contribution statement

Milan Mušo: Conceptualization, Methodology, Formal analysis, Investigation, Writing - original draft, Writing - review and editing, Visualization.

Rebecca Dumbell: Conceptualization, Writing-review and editing, Supervision.

Sara Pulit: Methodology, Formal analysis, Writing - review and editing.

Nasa Sinnott-Armstrong: Methodology, Software, Formal analysis.

Samantha Laber: Conceptualization, Investigation.

Louisa Zolkiewski: Validation, Investigation, Visualization.

Liz Bentley: Conceptualization, Supervision, Project administration.

Melina Claussnitzer: Conceptualization, Writing-review and editing.

Roger D Cox: Conceptualization, Writing-original draft, Writing-review and editing, Visualization, Supervision, Project administration and Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.