Author + information
- Received June 21, 2013
- Revision received August 30, 2013
- Accepted September 3, 2013
- Published online February 4, 2014.
- Kuixing Zhang, MD, PhD∗,
- Dekker C. Deacon, BS∗,
- Fangwen Rao, MD∗,
- Andrew J. Schork, BS∗,
- Maple M. Fung, MD∗,
- Jill Waalen, MD, MPH†,
- Nicholas J. Schork, PhD†,
- Caroline M. Nievergelt, PhD∗,
- Neil C. Chi, MD, PhD∗,‡ and
- Daniel T. O'Connor, MD∗,‡,§,‖∗ ()
- ∗Department of Medicine, University of California at San Diego, San Diego, California
- †Scripps Research Institute, La Jolla, California
- ‡Institute for Genomic Medicine, University of California at San Diego, San Diego, California
- §Department of Pharmacology, University of California at San Diego, San Diego, California
- ‖VA San Diego Healthcare System, San Diego, California
- ↵∗Reprint requests and correspondence:
Dr. Daniel T. O’Connor, Department of Medicine, UCSD School of Medicine, 9500 Gilman Drive, La Jolla, California 92093-0838.
Objectives The goal of this study was to understand the role of genetic variation in the catecholamine biosynthetic pathway for control of human heart rate (HR).
Background Human HR is an integrated cardiovascular trait predictive of morbidity and survival. Because the autonomic pathway exerts rapid control over the heart, we probed the role of heredity in the control of HR, focusing on a component of the autonomic sympathetic pathway already predictive of outflow responses: cytochrome b561 (CYB561), the electron shuttle in catecholamine vesicle membranes for transmitter biosynthesis.
Methods We studied hereditary control of HR with the twin pair design, at rest and during environmental (cold) stress. Single nucleotide polymorphism disruption of a microribonucleic acid (microRNA) recognition motif in the human CYB561 3′-UTR was identified computationally, and its differential effect on gene expression was demonstrated in a transfected luciferase reporter/3′-UTR variant. We exposed stem cell–derived human embryoid bodies to the microRNA mimic or antagomir oligonucleotides, and we observed the effects on contraction rate in proto-hearts.
Results Substantial heritability (h2) was demonstrated by using twin pair variance components for both basal/resting HR (h2 50.9 ± 6.4% of trait variation, p = 2.47 × 10−10) and stress-augmented HR (h2 55.1 ± 5.9%, p = 8.79 × 10−13), and the 2 HR traits shared genetic determination (genetic covariance ρG 0.747 ± 0.058, p = 2.85 × 10−9). CYB561 displayed 1 common genetic variant in the transcript region: A+1485G (rs3087776), in the 3′-UTR, 1485 bp downstream of the termination codon, in a conserved region, with the A-allele ancestral in primates. In a twin/sibling sample (n = 576), A+1485G influenced HR, both at rest (p = 0.010) and after environmental stress (p = 0.002), with the minor (A) allele displaying a recessive effect with lower HR. The effect of A+1485G on HR was extended by meta-analysis into 2 additional population samples (total n = 2,579), and the influence remained directionally consistent and significant (p = 0.007). A+1485G disrupted a microRNA (human microribonucleic acid-1294 [hsa-miR-1294]) recognition motif in the 3′-UTR, as demonstrated by a transfected luciferase reporter/human 3′-UTR variant system in 2 different neuronal/neuroendocrine cell types. The microRNA effect was further documented by cotransfection of an hsa-miR-1294 mimic, yielding an exaggerated decline in expression of the A-allele (better match) reporter (p = 4.3 × 10−5). Similar findings of differential 3′-UTR allelic susceptibility to hsa-miR-1294 were noted during expression of the full-length human CYB561 messenger ribonucleic acid with its cognate 3′-UTR. Finally, exposure of stem cell–derived human embryoid bodies to hsa-miR-1294 mimic or antagomir oligonucleotides yielded directionally opposite effects on contraction rate in proto-hearts.
Conclusions HR is a substantially heritable trait, with genetic influence by variation in the adrenergic pathway, here shown for messenger ribonucleic acid translational control at the CYB561 step of transmitter formation. The results have implications for potentially modifiable autonomic pathways that influence this risk trait in the population.
Human heart rate (HR) is an integrated cardiovascular trait that has assumed epidemiological importance because it may be predictive of premature mortality (both cardiovascular and noncardiovascular) and, consequently, lifespan (1,2). The HR trait is complex, reflecting not only cardiovascular but also metabolic adjustments (2). Rapid control of HR is achieved by the autonomic nervous system, in both its sympathetic (catecholaminergic, stimulatory) (3) and parasympathetic (vagal, inhibitory) branches. The role of heredity in control of HR has not been exhaustively examined, but genetic variation in the sympathetic pathway, including “tagging” (intronic) variation at the cytochrome b561 (CYB561) locus, reportedly influences cardiovascular responses to sympathetic activation (4). Indeed, naturally-occurring genetic variation at every other point in the adrenergic pathway (GCH1 , TH , DBH , PNMT , and ADRB1 ) has been associated with altered HR control, yet the interaction of CYB561 and HR is still unexplored.
CYB561 is an electron transfer protein unique to catecholamine and neuropeptide secretory vesicles of the adrenal medulla, pituitary gland, and other neuroendocrine tissues (10,11). The 30-kDa protein may comprise as much as ∼15% of the hormone storage vesicle membrane protein (12), and its role is to supply reducing equivalents to 2 monooxygenases, dopamine beta-hydroxylase (DBH) in chromaffin granules, and peptidylglycine alpha-amidating monooxygenase in neurosecretory vesicles (13). The cytochrome fulfills this role by catalyzing the transfer of electrons from a cytoplasmic donor, ascorbate, across a phospholipid bilayer to the luminal acceptor, semidehydroascorbate, in the interior of the vesicles. The continuously regenerated ascorbate within these vesicles is the immediate donor for the monooxygenases within the neuroendocrine secretory vesicles. Thus, cytochrome b561 is a transmembrane electron channel.
Because the sympathetic system exerts substantial control over HR, we first used the twin pair approach to determine whether HR responses were heritable; twin studies have successfully revealed hereditary contributions to cardiovascular stress traits, including HR (14). We then probed the role of genetic variation at the CYB561 locus in such responses. We next used a transfected reporter system to examine the role of common variation in the 3′-UTR of the gene, which disrupted a microribonucleic acid (microRNA) recognition motif, a potential control point for messenger ribonucleic acid (mRNA) translation. Finally, we explored whether microRNA influenced the contraction rate in cultured human embryoid bodies (EBs).
MicroRNAs are emerging as a widespread endogenous mechanism for control of gene expression at the posttranslational level (15), wherein the ∼22-nucleotide mature/processed microRNAs bind to particular motifs in the 3′-UTRs of target mRNAs, inhibiting translation by catalyzing transcript scission/degradation or by steric hindrance. Indeed, the human genome harbors more than 1,000 microRNA-encoding loci (http://www.microrna.org).
Subjects and characterization
Subjects were volunteers drawn largely from southern California, and each subject provided written informed consent; the protocol was approved by the institutional review board of the University of California at San Diego. Recruitment procedures, definitions, and confirmation of subject diagnoses are according to previous reports from our group.
Initial: twins and siblings
From 235 nuclear families, 576 individuals from twin and sibling pairs were recruited to conduct the following study. Zygosity was confirmed by extensive microsatellite and single nucleotide polymorphism (SNP) genotyping, as described previously (16). Twins ranged in age from 15 to 84 years. Individuals of white (European-American, 87%) or Hispanic (Mexican-American, 13%) biogeographic ancestry/ethnicity (according to self-identification) were included.
Physiological phenotyping in vivo (twins and siblings): environmental (cold) stress
To probe the functional significance of common variation at CYB561, we examined the potential influence of 1 common CYB561 polymorphism on HR before and during environmental (cold) stress testing (17) on 576 twin (monozygotic or dizygotic) and sibling individuals. Resting HR was recorded continuously for ∼5 min before cold stress. During the stressor, the subject immersed the nondominant hand into ice (0°C) water for 1 min, with averaged measurements of HR stable over at least 3 beats before and then again toward the end of the 60-s procedure.
Extension: primary care population samples
From a pool of >53,000 subjects in a primary care (Kaiser Permanente) database in southern California, we ascertained 2 samples (Kaiser-1 and Kaiser-2) of European-ancestry individuals of both sexes (18). Evaluation included physical examination (with vital signs), blood chemistry screening, hemogram, and medical history questionnaire. HR was measured by manual palpation of a radial artery for 30 s (then multiplying by 2 to obtain beats/min) in seated, resting subjects just before measurement of blood pressure.
Genotyping CYB561 variants
Bioinformatics of Polymorphism
The extent of polymorphism at the CYB561 locus was visualized at the National Center for Biotechnology Information (locus ID 1534), through GeneView at the Single Nucleotide Polymorphism database (http://www.ncbi.nlm.nih.gov/SNP/). The 3003-bp CYB561 mRNA (RefSeq clone NM_001915.3), which encodes a 251–amino acid protein product (Uniprot P49447), exhibits only 1 common (minor allele frequency [MAF] >5%) variant: rs3087776 A/G, located in the 3′-UTR, 1485 bp downstream from the stop codon (A+1485G), at MAF = 49.9%, heterozygosity = 0.499 (Fig. 1A). Sequence diversity results emerged from inspection of the following: the 1000 Genomes Project (http://www.1000genomes.org/) (chromosomes represented: 120 CEU [at 5.1× sequence coverage], 118 YRI, 120 CHB+JPT) and the HapMap (http://hapmap.ncbi.nlm.nih.gov/) (chromosomes represented: 226 CEU, 226 YRI, 86 CHB, 172 JPT). No common nonsynonymous (amino acid replacement at >1% frequency) variation has been observed at CYB561.
Genomic deoxyribonucleic acid (DNA) of each individual was prepared from leukocytes in ethylenediaminetetraacetic acid–anticoagulated blood by using PureGene extraction columns (Gentra Biosystems, Minneapolis, Minnesota). The 3′-UTR polymorphism rs3087776 A+1485G diploid genotypes were scored by using TaqMan technology (Applied Biosystems, Foster City, California) with a 5-μl reaction using 384 well plates: DNA 5 ng; probe (40×) 0.0625 μl; and master mix (2×) 2.5 μl. The polymerase chain reaction (PCR) profile was: 95°C (strand separation) for 10 min, followed by 40 cycles of 92°C for 15 s and 60°C for 1 min. Twins and siblings were genotyped with the Illumina610-Quad array (Illumina, Inc., San Diego, California), encompassing ∼590K SNP genotypes.
Computation and statistics
Heritability of Phenotype Expression in Vivo
In twin pairs, heritability (or the fraction of trait variance accounted for by genetic variance, heritability [h2] = VG/VP), as well as the genetic covariance (shared environmental determination or pleiotropy; parameter Rho_G or rG) and environmental covariance (shared environmental determination; Rho_E or rE), were estimated by variance components in SOLAR (Sequential Oligogenic Linkage Routines) (19), available at http://www.txbiomed.org/departments/genetics.
To test SNP on phenotype effects with explicit accounting for family structure, MERLIN version 1.1.2 (http://www.sph.umich.edu/csg/abecasis/merlin/) was used. As an additional quality control step, unlikely genotypes based on expected inheritance patterns were removed by using MERLIN's Pedwipe procedure. During the SNP association, subjects were categorized according to diploid genotype at the bi-allelic SNP locus. To control for additional genetic background heterogeneity in this predominantly white cohort, we performed a multidimensional scaling analysis by using PLINK (20), including all autosomal SNPs. We then included the first multidimensional scaling analysis dimension as covariate in the association analysis, which corresponded to the Native American admixture of Hispanic subjects (21,22).
Descriptive and inferential statistics
For twin/sibling analyses, descriptive (genotype-specific mean ± SEM) and inferential (chi-square, p value) statistics were computed across all of the twins using generalized estimating equations in SPSS version 17 (IBM SPSS Statistics, IBM Corporation, Armonk, New York) to account for correlated trait values within each sibship, using an exchangeable correlation matrix (23). After inspection of clustering of the descriptive statistics (mean ± SEM), heterozygotes and minor/major allele homozygotes could be grouped together to test dominant/recessive models. Estimates are stated as mean ± 1 SE. For unrelated individuals, 2-way analysis of variance or multivariable general linear modeling, as well as post-hoc corrections, were performed in SPSS version 17 to evaluate the significance of a single variant.
The coefficient of genotype effect (beta), its SE, and p values were obtained by using regression analysis in SPSS version 17. Meta-analyses were conducted with the command META, testing fixed effect (i.e., genotype) models in Stata version 12 (Stata Corp., College Station, Texas), after individual study regression analysis in SPSS, incorporating individual study data to derive significance as well as pooled genotype effect size (beta or slope per allele) and its SE. Analyses in twins were adjusted for age and sex (as well as ancestry); results in the Kaiser Permanente subjects were adjusted for age and sex.
We queried human mRNA expression data from the GTEx (Genotype Tissue Expression) program at its eQTL browser (http://www.ncbi.nlm.nih.gov/gtex/GTEX2/gtex.cgi), focusing on the effects of CYB561 genetic variation on CYB561 mRNA expression in neuronal tissues, originally available from Gibbs et al. (24). CYB561 3′-UTR variant A+1485G (rs3087776; chr-17, position 61,510,277) was tested by linear regression (assoc function within PLINK, with an additive genotype model) as a determinant of CYB561 mRNA abundance in a series of 143 human brain (frontal cortex) samples from neurologically normal white subjects, with available transcriptome (Illumina Human Ref-8 Expression bead-chips with probes for 22,184 mRNA transcripts; CYB561 probe ILMN_1771179, chr-17 position 60,214,268) and genome-wide association study data (Illumina Infinium HumanHap 550 bead-chips; 561,466 SNPs). Before regression, raw intensity values for mRNA expression were transformed by using the rank invariant normalization method.
CYB561 3′-UTR/luciferase reporter and activity assays
The CYP561 3′-UTR (2079 bp) was PCR-amplified from genomic DNA of known homozygotes and ligated into the unique XbaI site just downstream (3′) of the firefly luciferase reporter in plasmid pGL3-Promoter (Promega Corporation, Madison, Wisconsin), in which eukaryotic transcription is driven by the SV40 early promoter (Online Fig. 1). The variant 3′-UTR mutant G/A was created by site-directed mutagenesis (QuikChange, Stratagene, La Jolla, California) and verified by using dideoxy sequencing. Cells were transfected (at ∼50% to 60% confluence in 15.6-mm polystyrene dishes) with 500 ng of supercoiled plasmid by using the liposome method (TransFectin, Bio-Rad Laboratories, Hercules, California). After transfection and cell growth over a 12- to 48-h time course, cells were treated with passive lysis buffer (catalog no. E194A, Promega Corporation) for sequential measurement of luciferase enzymatic activity and protein concentration, and the results were expressed as the ratio of firefly luciferase/protein, as described previously (25,26). Each cellular experiment was repeated a minimum of 4 times. Results were expressed as mean ± SEM. Statistical significance (p < 0.05) was calculated by using an analysis of variance.
Cell and molecular biology
Neuronal or neuroendocrine cells, including rat pheochromocytoma cells (PC12), human HEK-293T cells, or human neuroblastoma cells (SH-SY5Y), were grown at ∼50% to 60% confluence in 15.6-mm polystyrene dishes before transfection. HEK-293T cells were chosen because of their prominent neuronal gene expression phenotypes (27). mRNA (e.g., CYB561) abundance was quantified by using quantitative real-time PCR and normalized to expression of endogenous “housekeeping” RNAs (glyceraldehyde-3-phosphate dehydrogenase, beta-actin, or 18S rRNA) in the same sample, according to the ΔΔCt method (28). Endogenous microRNA (i.e., human microribonucleic acid-1294 [hsa-miR-1294]) abundance was also quantified by using quantitative real-time PCR, but normalized to an endogenous small RNA (SNORD61). The full-length human CYB561 complementary DNA (NM_001915, including the 753 bp [251 amino acid] ORF with 2079 bp 3′-UTR) was obtained from Open Biosystems (Huntsville, Alabama) in the eukaryotic (cytomegalovirus promoter) expression plasmid pSPORT6; site-directed mutagenesis (QuikChange, Stratagene) created the 3′-UTR G/A variant, verified by using dideoxy sequencing. These plasmids were then transfected into PC12 or HEK-293T cells, with or without cotransfected hsa-miR-1294 oligonucleotides (double-stranded mimic or control), to test whether expression of the human CYB561 mRNA was influenced by the A+1485G variant in its cognate 3′-UTR.
Evaluation of 3′-UTR microRNA (hsa-miR-1294) motif and function
Interspecies sequence alignments were done by using the ClustalW program. 3′-UTR microRNA motifs were predicted at RegRNA (29) (http://regrna.mbc.nctu.edu.tw/). The effects of the variant on ribonucleic acid (RNA) hybrid structures and predicted minimum folding energies were analyzed by using BiBiServ (30) (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/submission.html).
An initial strategy of microRNA overexpression was adopted to explore microRNA function in cells. A Dharmacon (Lafayette, Colorado) miRIDIAN 22-mer (UGUGAGGUUGGCAUUGUUGUCU) double-stranded RNA mimic for hsa-miR-1294 (catalog no. C-301348-00) was cotransfected (TransFectin, Bio-Rad Laboratories) to increase that miRNA abundance, while a pre-designed negative control (No. 1, catalog no. CN-001000-01-05) was cotransfected to control for off-target effects. We then performed hsa-miR-1294 microRNA antagonist (antagomir, antisense version) studies, with a single-stranded oligonucleotide (Dharmacon IH-301348-01-0005). The identities of synthetic oligonucleotides were verified by using matrix-assisted laser desorption-ionization time-of-flight mass spectrometry, with confirmation of duplex integrity by PAGE. Synthetic oligonucleotides were applied to a final concentration of 15 nM, and cotransfected into PC12 or HEK-293T cells, along with wild-type or variant CYB561 3′-UTR/luciferase reporters, by using the liposome method described earlier. For studies in human EBs, we also transfected (TransFectin, Bio-Rad Laboratories) a synthetic, single-stranded hsa-miR-1294 antagomir (antisense version, or inhibitor; Dharmacon IH-301348-01-0005).
Human EB cultures (for contraction rate effects of hsa-miR-1294 in vivo)
Human embryonic stem cells (ESCs) (line H9) were grown on mouse embryonic fibroblasts in WiCell (Madison, Wisconsin) hESC medium (with 12 ng/ml basic fibroblast growth factor). The ESCs were differentiated to cardiomyocytes in EBs as described previously (31). After 6 weeks, cells were transfected for 24 h with either 20-nM hsa-miR-1294 double-stranded mimic, single-stranded antagomir (inhibitor), or control oligonucleotide (8 replicate wells per condition) by using the liposome method (TransFectin, Bio-Rad Laboratories). The contraction rate for each embryoid body was counted immediately before transfection and again 24 h after transfection. Due to the usual heterogeneity in contraction rate among EBs in all conditions, the contraction rate was compared for each EB before and after transfection. Because contraction rate values were not normally distributed, results were evaluated by using the nonparametric Mann-Whitney test.
Located on chromosome 17q23, the CYB561 locus spans ∼14 kbp, with only 1 common (MAF >5%) variant in the entire encoded transcript: 3′-UTR polymorphism A+1485G (rs3087776), with an MAF ∼49% (Fig. 1A). Linkage disequilibrium across the locus, reflecting substantial marker-on-marker correlations, indicates that the entire gene is contained within a single block (Online Fig. 1).
CYB561 3′-UTR Variant A+1485G (rs3087776) Conserved Across Primate Species
The A+1485G variant is located in a relatively conserved 3′-UTR region across primate species (human, chimpanzee, orangutan, rhesus, and marmoset) (Fig. 1B). The minor (less frequent) A-allele seems to be ancestral because it is present in other primates, including chimpanzee. The A-allele displays a better match for the hsa-miR-1294 motif than does the G-allele (Fig. 1C).
Human HR: heritability, stress responses, and influence of CYB561 3′-UTR variant
HR h2 in Twins
In the overall group, HR rose by a mean of 3.2 ± 0.39 beats/min (p = 1.92 × 10−13) (Table 1). Variance component analyses (Table 2) indicated that basal (resting) HR was under substantial genetic control, with h2 = 50.9 ± 6.4% of trait variance (p = 2.47 × 10−10), as was HR post-stress, at h2 = 55.1 ± 5.9% (p = 8.71 × 10−13) (Fig. 2A). The change during stress was somewhat less heritable, although still highly significant (h2 = 40.7 ± 6.9%, p = 1.00 × 10−7). Basal and post-stress HR were directly correlated (Spearman r = 0.682, p = 2.73 × 10−75) and shared substantial genetic codetermination (or pleiotropy; rG = 0.747 ± 0.058, p = 2.85 × 10−9). By contrast, change in HR was inversely correlated with basal HR (Spearman r = –0.299, p = 1.94 × 10−13), and the 2 shared principally environmental codetermination (rE = –0.443 ± 0.068, p = 2.46 × 10−8).
CYB561 3′-UTR Polymorphism A+1485G (rs3087776): Effects on Heritable HR Traits
In the twin/sibling sample, the 3′-UTR A+1485G polymorphism predicted both basal (p = 0.010) and post-stress (p = 0.002) HR (Table 3, Fig. 2B). Minor allele (A/A) homozygotes displayed clearly lower values for both basal and stressed HR, suggesting a recessive mode of inheritance of the A-allele on HR traits. Change in HR (delta, pre-→post-stress) was not influenced by A+1485G (p = 0.752).
We previously reported (4) that a “tagging” genetic variant at CYB561 (rs2058203, intron-1) predicted the human vasoconstrictive response to endogenous in vivo catecholamine release. This CYB561 tagging variant (rs2058203, intron-1) also predicted HR in the current study, with an identical p value to rs3087776 (3′-UTR), which perhaps not surprising in that the 2 SNPs are separated by only 12,726 bp and are in perfect (100%, R2 = 1.0) linkage disequilibrium in the current subjects.
Extension Studies of CYB561 rs3087776 on Basal HR in Additional Population Samples
Meta-analysis combining the twin/sibling pairs and 2 additional independent population samples (Kaiser-1 and Kaiser-2) indicated allelic effects consistent in direction (negative sign on slope) across groups; the overall slope of the meta-analysis regression (beta) was –0.899, with an SE (of beta) of 0.327 (significant at p = 0.007) (Table 4). Congruent with the consistent negative slope in each group, there was no demonstrable heterogeneity across the groups (Q = 2.212, df = 2 [i.e., n – 1], p = 0.331); nonetheless, inspection of individual study p values suggests that the twin/sibling study is the most important contributor to the overall meta-analysis results.
Functional studies of the trait-associated variant: human CYB561 3′-UTR variant and expression in transfected cells
CYB561 3′-UTR Reporter Activity Assays
To probe the functional significance of common variant rs3087776, we inserted each of the 2 versions (G-allele vs. A-allele) into the luciferase reporter plasmid pGL3-Promoter, just downstream (3′) from the luciferase open reading frame (Online Fig. 2). After transfection into rat chromaffin (PC12) cells, the 2 alleles yielded substantially different luciferase reporter activities (A<G) at each time point (Fig. 3A).
CYB561 3′-UTR Variant A+1485G: Coordinate Directional Functions in Cells and in Vivo (Twins/Sibling Pairs)
Taking advantage of in vivo (Fig. 2B) and cellular (Fig. 3A) data, we observed that the A/A genotype, which predicted lower HR in the twins/sibling pairs, also displayed lower gene expression (Fig. 3B).
MicroRNA Effects: hsa-miR-1294 Mimicry and Inhibition on Luciferase Reporter/Human CYB561 3′-UTR Transfections
Because A+1485G disrupted an hsa-miR-1294 recognition motif with a superior match for the A-allele (Fig. 1C), we probed the functional significance of the match with a specific microRNA mimic and inhibitor (compared with control oligonucleotide) cotransfected with the luciferase/human 3′-UTR allele reporter plasmids into HEK-293T or PC12 cells. In HEK-293T cells (Fig. 3C), the alleles were expressed differentially (G>A, p < 0.001), and exogenous/cotransfected hsa-miR-1294 mimic decreased reporter expression more effectively on the A-allele (p = 0.007). Likewise, PC12 chromaffin cells (Fig. 3D) displayed an expression difference between alleles (G>A, p = 0.028), and the difference was amplified on the A-allele (57.1% decline [p = 0.001] vs. 51.5% for the G-allele [p = 0.006]). The substantial difference in expression after exposure of each allele to the microRNA (p = 4.3 × 10−5) further documents the role of hsa-miR-1294. Upon exposure of transfected HEK-293T cells to an miR-1294 inhibitor (antagomir) (Fig. 3E), reporter expression was increased more effectively on the A-allele (by 11%; p = 0.004).
MicroRNA Effects: hsa-miR-1294 Mimicry on Human CYB561 mRNA Expression With Its Cognate 3′-UTR
When the full-length human CYB561 with its cognate 3′-UTR was expressed in HEK-293T cells (Fig. 4A), the 3′-UTR alleles were expressed differentially (G>A, p = 0.045), and the miR-1294 mimic preferentially inhibited the A-allele (p = 0.036). Likewise, upon expression in PC12 chromaffin cells (Fig. 4B), the miR-1294 mimic preferentially inhibited the 3′-UTR A-allele (p = 0.034), although basal expression of the 2 alleles was similar.
Endogenous RNAs in Model Neuronal/Neuroendocrine Cells: CYB561 mRNA and miR-1294 Small RNA
We quantified endogenous abundance of the pertinent RNAs by using quantitative real-time PCR in untransfected neural/neuroendocrine cells (PC12, SH-SY5Y, and HEK-293T), normalizing mRNA results to 18S rRNA and microRNA results to small RNA SNORD61 (Fig. 4C). In general, there was an inverse proportionality for the mRNA and microRNA (R = 0.998, p = 0.0001), with the mRNA reaching a plateau at low microRNA abundance, suggesting the expected dependence of target mRNA expression on the specific microRNA.
A+1485G (rs3087776) predicted CYB561 mRNA expression in human brain frontal cortex, at p = 3.14 × 10−8, accounting for ∼20% of the variance of mRNA abundance in the 143 samples (R2 = 0.1958). No trans-QTLs were detected for CYB561, even down to a permissive p < 0.01 threshold.
Contraction rate effects of hsa-miR-1294 in vivo: human embryoid body cultures
When exposed to 20-nM miR-1294 mimic oligonucleotide for 24 h, the proto-hearts decreased the contraction rate by 5.6 ± 1.4%, whereas the miR-1294 antagomir (inhibitor, again at 20 nM) elevated the contraction rate by 8.2 ± 0.7%. The result was a significant difference between treatments (p = 0.014) (Fig. 5).
HR not only predicts survival (1,2) and reflects cardiovascular performance, but it also serves as an integrated indicator of autonomic nervous (sympathetic and parasympathetic/vagal) system activity and metabolic rate (2).
Genome-wide association studies have now identified up to 21 genetic loci that contribute significantly to resting HR (32). Although 2 of these loci are in the classical parasympathetic pathway (ACHE, CHRM2), none of the genes encode catecholamine biosynthetic enzymes or adrenergic receptors, and the identified loci so far account for only ∼0.9% of trait variance. By contrast, twin studies indicate that heredity contributes up to ∼50% of HR variability; indeed, our twin results suggest HR heritability of 50.9 ± 6.4% of trait variance (p = 2.47 × 10−10) (Fig. 2A). The discrepancy between the 2 approaches is sometimes referred to as the “missing heritability problem” (33) and prompts a continuing search for responsible loci, in this case exploring the role of genetic variation in a contributory physiological pathway.
The autonomic system serves as a “master regulator” of diverse organ function, and autonomic dysregulation can manifest in a variety of ways. Because increased sympathetic tone may contribute to HR elevation (3), and “tagging” genetic variation at the CYB561 locus may predict responses to pre-synaptic stimulation in humans (4), we examined the association of sympathetic pathway gene CYB561 with HR traits, focusing on functional (rather than simply “tagging”) genetic variation. It was known that genetic variation across the adrenergic pathway (GCH1 , TH , DBH , PNMT , and ADRB1 ) predicted control of HR; thus, we focused on CYB561. We found that common genetic variation in the CYB561 was confined to the 3′-UTR at A+1485G, which disrupted a microRNA recognition motif. In twin pairs, A+1485G predicted both basal and stress-augmented HR, 2 highly heritable traits (Fig. 2).
CYB561 (cytochrome b561)
Cytochrome b561 (Uniprot P49447; named because of its optical absorbance at 561 nm) is a monomeric, 251-amino acid, ∼28 kDa protein with 6 transmembrane-spanning domains that functions as an electron shuttle for DBH in the catecholamine vesicle (chromaffin granule) membrane, wherein it is essential to catecholamine biosynthesis and, hence, the effects of these amines (34). CYB561 subserves catecholamine biosynthesis by accepting electrons from semidehydroascorbate and thus regenerating ascorbate for its role as electron acceptor (as enzyme cofactor) during oxidative conversion of dopamine to norepinephrine by DBH within catecholamine storage vesicles. Because the ascorbate/dehydroascorbate shuttle also subserves the action of the endocrine secretory vesicle enzyme peptidyl-alpha-amidating monooxygenase, the synthesis of carboxy-terminal amidated peptide transmitters (e.g., neuropeptide Y, pancreastatin, peptide histidine isoleucine, peptide tyrosine tyrosine, and neuropeptide K) is also dependent on CYB561 (35). CYB561 has also been described as an electron shuttle in macrophage lysosomal membranes (36).
The di-heme-binding CYB561 amino acid sequence is highly conserved across species (37,38) with 6 transmembrane helices, of which helices 2 through 6 are located closely in the same regions of the 26 sequences in the alignment. Indeed, we observed that the local region of the CYB561 3′-UTR harboring the hsa-miR-1294 motif is also substantially conserved (Fig. 1B).
Role of microRNA miR-1294 motif in the CYB561 3′-UTR
Isolation of the CYB561 3′-UTR variant onto a luciferase reporter allowed functional studies, confirming the role of A+1485G in differential responses to hsa-miR-1294 (Fig. 3), thus providing a mechanistic basis for the clinical findings. The differential allelic effects were confirmed during eukaryotic expression of full-length human CYB561 mRNA with its cognate 3′-UTR (Fig. 4). The apparent recessive effect of A/A homozygosity on HR traits (Fig. 2B) suggests that the A-allele may be rate limiting in determining the formation or effect of CYB561_mRNA:miR-1294 complexes, a proposal consistent with the superior match of the A-allele motif with miR-1294.
Differential 3′-UTR allelic (A/G) effects on CYB561 mRNA abundance (Figs. 4A and 4B), despite the use of identical promoters (cytomegalovirus promoter) for cDNA/mRNA expression, suggest that a decline in mRNA stability is likely for the A-allele, in the face of the enhanced microRNA:mRNA match for the A-allele; alterations in mRNA translation by the microRNA are also likely as a result of steric hindrance of the translational machinery by the mRNA-bound microRNA. The RNA components of this system (CYB561 mRNA, and miR-1294 small RNA) were endogenously coexpressed in model cell lines. Finally, exposure of human EBs to miR-1294 mimic or antagomir oligonucleotides yielded directionally opposite effects on contraction rate in beating proto-hearts (Fig. 5).
Study strengths and limitations
We took advantage of the twin pair approach to establish heritability of HR traits (Tables 1 to 3, Fig. 2) and then characterized the effect of sympathetic pathway variation at the CYB561 locus on such traits. We then observed directionally similar CYB561-on-trait effects in 2 additional population samples, with preserved significance (Table 4). We were later able to confirm the functional effects of 3′-UTR variation on responses to a particular microRNA by using a transfected chimeric 3′-UTR reporter system (Fig. 3), as well as the cognate (CYB561) mRNA (Fig. 4), and finally we observed reciprocal effects of miR-1294 mimic and antagomir on contraction rate in human EBs in vivo (Fig. 5). Nonetheless, additional questions arise during this work. For example, CYB561 acts as an electron shuttle in service to not only DBH but also peptidyl alpha-amidating monooxygenase and, hence, might exert effects through peptidergic in addition to catecholaminergic pathways. In addition, we have not yet documented the effects of miR-1294 on HR in freely living organisms. Finally, hsa-miR-1294 likely has additional mRNA targets in the transcriptome; indeed, algorithms such as TargetScan (http://www.targetscan.org/), requiring only very short microRNA:mRNA matches at 7 to 8 bp microRNA “seed sequences” (39), identify many potential mRNA targets for hsa-miR-1294; the functional importance of such very short, partial matches at these loci has not been explored.
A very common (MAF ∼45% to 49%) genetic variant in the CYB561 3′-UTR alters responses to a specific microRNA (hsa-miR-1294), ultimately leading to heritable changes in both basal and stress-augmented HR. Online Figure 3 presents a hypothetical schema outlining a step-by-step pathway in series that may unify our observations, as well as suggest avenues for further experimentation.
For supplemental figures, please see the online version of this article.
This research was funded by the National Institutes of Health (HL58120; DK094894; UL1RR031980 [UCSD Clinical and Translational Research Institute]; MD000220 [UCSD Comprehensive Research Center in Health Disparities]; 5U01HL107442) and the Department of Veterans Affairs. Dr. Fung is employed by and owns stock in Amgen Inc. All other authors have reported that they have no relationships relevant to the contents of this paper to disclose. Ali J. Marian, MD, served as Guest Editor for this paper.
- Abbreviations and Acronyms
- single nucleotide polymorphism (rs3087776) in the human CYB561 3’-UTR
- cytochrome b561 (electron shuttle for dopamine beta-hydroxylase)
- dopamine beta-hydroxylase
- deoxyribonucleic acid
- embryoid body
- embryonic stem cells
- human micro-ribonucleic acid-1294
- heart rate
- minor allele frequency
- microribonucleic acid
- messenger ribonucleic acid
- polymerase chain reaction
- single nucleotide polymorphism
- 3′ (downstream of the open reading frame) untranslated region (of the human CYB561 messenger ribonucleic acid)
- Received June 21, 2013.
- Revision received August 30, 2013.
- Accepted September 3, 2013.
- American College of Cardiology Foundation
- Greenland P.,
- Daviglus M.L.,
- Dyer A.R.,
- et al.
- Fung M.M.,
- Nguyen C.,
- Mehtani P.,
- et al.
- Zhang L.,
- Rao F.,
- Wessel J.,
- et al.
- Garland E.M.,
- Black B.K.,
- Harris P.A.,
- Robertson D.
- Jirout M.L.,
- Friese R.S.,
- Mahapatra N.R.,
- et al.
- Duong L.T.,
- Fleming P.J.
- Srivastava M.,
- Duong L.T.,
- Fleming P.J.
- Wu T.,
- Treiber F.A.,
- Snieder H.
- Rao F.,
- Zhang L.,
- Wessel J.,
- et al.
- Shaw G.,
- Morse S.,
- Ararat M.,
- Graham F.L.
- Huang H.Y.,
- Chien C.H.,
- Jen K.H.,
- Huang H.D.
- Kruger J.,
- Rehmsmeier M.