Gene Expression Networks in the Murine Pulmonary Myocardium Provide Insight into the Pathobiology of Atrial Fibrillation

The pulmonary myocardium is a muscular coat surrounding the pulmonary and caval veins. Although its definitive physiological function is unknown, it may have a pathological role as the source of ectopic beats initiating atrial fibrillation. How the pulmonary myocardium gains pacemaker function is not clearly defined, although recent evidence indicates that changed transcriptional gene expression networks are at fault. The gene expression profile of this distinct cell type in situ was examined to investigate underlying molecular events that might contribute to atrial fibrillation. Via systems genetics, a whole-lung transcriptome data set from the BXD recombinant inbred mouse resource was analyzed, uncovering a pulmonary cardiomyocyte gene network of 24 transcripts, coordinately regulated by chromosome 1 and 2 loci. Promoter enrichment analysis and interrogation of publicly available ChIP-seq data suggested that transcription of this gene network may be regulated by the concerted activity of NKX2-5, serum response factor, myocyte enhancer factor 2, and also, at a post-transcriptional level, by RNA binding protein motif 20. Gene ontology terms indicate that this gene network overlaps with molecular markers of the stressed heart. Therefore, we propose that perturbed regulation of this gene network might lead to altered calcium handling, myocyte growth, and contractile force contributing to the aberrant electrophysiological properties observed in atrial fibrillation. We reveal novel molecular interactions and pathways representing possible therapeutic targets for atrial fibrillation. In addition, we highlight the utility of recombinant inbred mouse resources in detecting and characterizing gene expression networks of relatively small populations of cells that have a pathological significance.

cells are difficult to analyze in situ and have not been extensively studied in the endogenous context.
Analyses of the prevalence, length, and thickness of myocardial extensions has resulted in conflicting reports as to whether various parameters do (Hassink et al. 2003;Kholová and Kautzner 2003) or do not (Kholová and Kautzner 2004;Saito et al. 2000) correlate with atrial fibrillation. However, mutations in genes encoding the transcription factors Nkx2-5 (Huang et al. 2013;Xie et al. 2013), Pitx2c (Qiu et al. 2014;Wang et al. 2014), and GATA6 ) are all associated with atrial fibrillation; these genes are known to be required for cardiovascular development. These data suggest that perturbation of the transcriptional program associated with cardiac muscle (but not necessarily disrupting morphology) may underlie the association of the pulmonary myocardium with atrial fibrillation. Consistent with this, both PITX2C and NKX2-5 are required to direct and maintain the identity of the pulmonary myocardium (Mommersteeg et al. 2007). Lost or reduced expression of NKX2-5 causes pulmonary myocardium cells to adopt a pacemaker-like phenotype similar to cells of the sinus horn myocardium, without altering cellular morphology (Mommersteeg et al. 2007;Ye et al. 2015). This implies that healthy pulmonary myocardium could gain pacemaker function from which an ectopic beat might originate via a shift in the gene expression program away from that of working myocardium and toward a pacemaker program. Thus, molecular analysis of pulmonary myocardium cells in vivo is warranted and could provide a powerful tool for understanding the pathobiology of atrial fibrillation and identifying novel therapeutic targets. This led us to investigate gene regulatory mechanisms that might drive the expression of cardiac genes in this tissue.
In this study, we used a systems-genetics approach to conduct complex trait analysis and uncover gene networks in the pulmonary myocardium in situ (Chesler et al. 2005;Sieberts and Schadt 2007). Systems-genetics approaches use genetic reference populations such as the BXD resource, a set of recombinant inbred mice which are the progeny of an intercross between the two parental mouse strains C57BL/6J and DBA/2J (Andreux et al. 2012;Peirce et al. 2004). Recombinant inbred strains like the BXD and the Collaborative Cross (Ram et al. 2014;Threadgill et al. 2011) are especially suited for systems-genetics studies because each strain is effectively immortal (Andreux et al. 2012). They provide a powerful test bed for linking genetic loci to variation, including mapping loci controlling gene expression levels. Such loci are termed expression quantitative trait loci (eQTL) (Chesler et al. 2005). Importantly, analyses of eQTL provide a means for detecting specific regulators of a gene of interest as well as identification of networks of coregulated transcripts (Andreux et al. 2012;Hall et al. 2014).
We have analyzed a publicly available BXD steady-state, whole-lung transcriptome data set generated from 51 BXD strains (Alberts et al. 2011) to interrogate the molecular phenotype of the pulmonary myocardium. Based on the assumption that genes with highly correlated expression in a tissue are likely to act in a common network or biological process, Alberts et al. (2011) previously identified gene networks associated with specific B-and T-cell populations from this same whole-lung data set. We reasoned that we could use a similar approach with the lung data set to identify and characterize pulmonary myocardium cells using the cardiac troponin I type 3 (Tnni3) gene (previously shown to be expressed in pulmonary myocardium) (Millino et al. 2000) by investigating highly correlated transcripts associated with cardiac muscle. This yielded 24 transcripts significantly correlated with the expression of Tnni3, which form the basis of a pulmonary myocardium gene network.

MATERIALS AND METHODS
Analysis of public data sets The GeneNetwork suite of online programs, incorporating WebQTL (Chesler et al. 2004;Wang et al. 2003), was used to analyze the publicly available BXD lung transcriptome data set (GN160) (www.genenetwork. org/webqtl/main.py) generated by Alberts et al. (2011). Genome-wide association analyses were performed using gene expression microarray values measured on a log 2 scale. Correlation analysis was performed with Pearson's correlation, and outliers were removed. For interval mapping, likelihood ratio statistic (LRS) scores were plotted on the y-axis vs. chromosomal location on the x-axis (autosomes and X chromosome). An LRS score .17 was significant and an LRS score .10 was suggestive at a genome-wide P value of ,0.05 using 5000 permutations. LRS scores were computed using Haley-Knott regression (Haley and Knott 1992). The genomic intervals analyzed were selected based on the location at which the eQTL peak crossed the suggested or significant LRS threshold. Principal component analysis (PCA) was used to merge multiple gene expression traits into a single, synthetic trait that represented a cellspecific gene expression signature.

Functional annotation and enrichment analysis
First-pass assessment of gene function was obtained from the GeneCards database (www.genecards.org) (Safran et al. 2010) and a tissue expression profile was analyzed using the Fantom5 database (http://fantom. gsc.riken.jp/zenbu/) (Forrest et al. 2014;Severin et al. 2014). The top 100 transcripts (mapping to 78 unique genes) that correlated with the pulmonary myocardium PCA trait were submitted for enrichment analyses and gene ontology (GO) using the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) . This resource incorporates information from different public databases including the Kyoto Encyclopedia of Genes and Genomes (KEGG), WikiPathways, and human phenotype ontology. We also used the oPOSSUM-3 Webbased tool (Kwon et al. 2012) to identify overrepresentation of transcription factor binding sites and integrated ChIP-seq data from two previously published data sets (Dupays et al. 2015;He et al. 2011).

Immunohistochemistry
Mouse experiments were approved by the University of Western Australia's Animal Ethics Committee. Four $8-wk-old male BALB/c mice were killed with an overdose of pentobarbitone (160 mg/kg via intraperitoneal injection) and fixed by cardiac perfusion with heparinized phosphate buffered saline (PBS) (15 IU/ml; heparin) followed by fixative (2% paraformaldehyde with 0.2% picric acid in PBS, pH 7.4). Heart and lungs were excised and postfixed for 48 hr at 4°in the same fixative. Tissue was washed, dehydrated, and embedded in paraffin wax. Five-micrometerthick sections through lungs or ventricular myocardium (control) were cut and mounted onto Superfrost Plus microscope slides (Sigma-Aldrich).
The distribution of TNNI3 in BALB/c lung tissue was detected with a mouse monoclonal IgG 2a antibody (Abcam, clone 4C2). Briefly, sections were dewaxed, rehydrated, and subjected to 20 min of heatinduced epitope retrieval (HIER) (pressure cooker) in Tris buffered saline (TBS)/EDTA buffer, pH 9.0 (Vector Laboratories). After 15 min cooling at room temperature, sections were washed in TBS (pH 7.4) and incubated in anti-TNNI3 primary antibody or mouse IgG isotype control (2 mg/ml in 3% fish skin gelatin in TBS) overnight at 4°. Sections were then washed in TBS for 1 hr before incubating in Mouse on Mouse Polymer (Abcam) for 30 min, washed in TBS, and incubated with diaminobenzidine (0.4 mg/ml, 3 min). Sections were counterstained with Mayers hematoxylin, blued with Scotts Tap Water Substitute, washed and dehydrated through graded alcohols to xylene, then coverslipped with Depex mounting medium. Cardiac a-actin (ACTC1) was detected with a mouse monoclonal IgG 1 antibody (Simga-Aldrich, clone Ac1-20.4.2), or mouse isotype control at 5 mg/ml, with the same protocol as above but with HIER in citrate buffer with pH 6.0 (Vector Laboratories). The presence of smooth muscle a-actin (ACTA2) was detected in a similar fashion, using mouse monoclonal IgG 2a antibody (Sigma-Aldrich, clone 1A4) at 9 mg/ml, and detection as above, but without HIER. Tissue was permeablized with 1% Triton X-100 for 15 min prior to the application of the primary antibody. Digital images were acquired on an Aperio ScanScope XT digital slide scanner (Leica Technologies). Images are typically of the left lobe and right middle lobe of four mice.

Data availability
All data used for this project are publicly available and accessible online (GN accession no. GN160).

RESULTS
Tnni3 can be used to identify a pulmonary myocardium gene network Since the expression of TNNI3 has previously been used as a marker for pulmonary myocardium (Millino et al. 2000;Mommersteeg et al. 2007;Ye et al. 2015), we reasoned that genes with expression that correlated with Tnni3 expression (trait ID 1422536_at, mean expression 11.552) Figure 1 Identification of a pulmonary myocardium gene network. The top 100 transcripts that correlated with Tnni3 in the BXD whole-lung transcriptome data set were used to generate a network graph. This resulted in a network of 24 transcripts with a Pearson correlation .0.8. (A) These 24 transcripts were subsequently regraphed to interrogate relationships at Pearson correlation values between 0.7 and 1.0 (red lines) and between 0.5 and 0.7 (orange lines). (B) The functional relationships of the 21 proteins translated from these 24 transcripts were also identified. All proteins except DOC2G and MYBPHL had known roles in cardiac function including structural sarcomeric proteins, regulation of sarcomere assembly, ion transport, transcriptional and post-transcriptional regulation, and hormone signaling.
in the whole-lung data set would also be expressed in the pulmonary myocardium. We performed sample correlation in the BXD data set and returned the top 100 significantly correlated transcripts (Supplemental Material, Table S1 in File S2; P , 0.0005). This set of 100 transcripts was enriched for genes known to be associated with cardiac muscle function, including structural proteins (actin and myosin), ion transporters (ryanodine receptors, calcium, and potassium transporters), and transcription factors (Hand2 and Nkx2-5). To define genes within this list that were likely to act in a concerted manner in the pulmonary myocardium, we generated a network graph and found 24 transcripts that were significantly associated at a Pearson correlation .0.8. These 24 transcripts were subsequently regraphed with lower stringency to visualize potential relationships between them ( Figure 1A). Three transcripts (Tnnt2, Pln, and Kcnj3) were detected by two separate probes and thus the 24 transcripts corresponded to only 21 unique genes. All 21 genes have previously been reported to be expressed in cardiac muscle (Chopra and Knollmann 2013;Ding et al. 2010;Fukuda and Mikoshiba 2000;Greulich et al. 2011;Guo et al. 2013;Hennessey et al. 2013;Holmegard et al. 2010;Huh et al. 2012;Marx et al. 2000;Reiser et al. 2001;Sequeira et al. 2015;Shaikh et al. 2016;Stobdan et al. 2015;Takeda et al. 2003;Warren et al. 2012;Yamashita et al. 2003;Zhou and Wu 2014) ( Figure 1B and Table 1); all were highly significantly correlated with Tnni3 (P , 3.02e 28 ). These genes can be loosely clustered into functional categories including structural sarcomeric proteins, regulation of sarcomere assembly, ion transport, transcriptional and post-transcriptional regulation, and hormone signaling (Table 1).
To ensure that the muscle-related transcripts that we had identified were indeed associated with the pulmonary myocardium and not the smooth muscle present in small arteries and surrounding the bronchus, we performed correlation analysis using Acta2. The top 100 transcripts correlated with Acta2 were enriched for smooth muscle-associated genes [including gamma-actin 2 (Actg2), leiomodin (Lmod2), transgelin (Tagln), and calponin 1 (Cnn1)], and did not significantly overlap with Tnni3-correlated transcripts (Table S2 in File S2).
n Table 1 Genes coexpressed with Tnni3 in the BXD whole-lung transcriptome have known roles in the function and regulation of the heart Functional Category Gene Name Function Sarcomeric structural protein Tnni3 Cardiac troponin I The troponin complex couples calcium availability with muscle contraction (Takeda et al. 2003) Tnnt2 Cardiac troponin T2 Myh6 Cardiac myosin heavy chain 6 (a-MHC) Two myosin heavy chains, two essential myosin light chains and two regulatory myosin light chains interact to form the thick filament of the sarcomere (Reiser et al. 2001;Yamashita et al. 2003)

Myh7
Cardiac myosin heavy chain 7 (b-MHC) Myl4 Atrial essential myosin light chain (ALC1) Myl7 Atrial regulatory myosin light chain (MYLC2a) Mybpc3 Cardiac myosin binding protein C Modulates the interaction between actin and myosin to regulate power output (Sequeira et al. 2015 In sum, the identity of the Tnni3-correlated genes themselves, and their distinct expression compared to smooth muscle-associated genes, strongly supports that we have identified a network of genes expressed in the pulmonary myocardium. The pulmonary myocardium gene network is coregulated from loci on Chr1 and Chr2 We next investigated whether transcripts within the pulmonary myocardium gene network were coregulated. Pair-wise comparisons were performed using a matrix function. All probe pairs within this set of 24 transcripts were positively correlated ( Figure S1 in File S1). All pairs had a Pearson correlation .0.5 with the exception of Myl4 and Doc2G (0.486), and Myl4 and Rbm20 (0.495). Further, heat-map analysis of eQTL for all 24 transcripts highlighted prominent regions of shared covariance on Chr1 and 2 with possible minor regions on Chr3, Chr4, Chr13, and Chr17 ( Figure 2A). This suggested that expression levels of multiple transcripts in this gene network were coordinately regulated by a shared set of cis and trans eQTL.
To identify the genomic intervals harboring genes that mediated variable expression, we performed genome-wide interval mapping and marker regression analyses for each transcript. The expression level of all transcripts varied by greater than twofold between BXD strains ( Table  2). All transcripts except one (Sln, sarcolipin) had either a suggestive (LRS .10.4) or significant (LRS .17) eQTL that mapped to the hotspot regions on Chr1 and Chr2 identified in the cluster map (Table  2). One transcript, Myh6 (myosin heavy chain 6, Chr14), varied by  2.62-fold between the highest and lowest expressing BXD strains (Figure 2B), such that a significant trans eQTL could be mapped to Chr1: 33 6 1 Mb ( Figure 2C). Of the remaining transcripts, 11 also had a suggestive eQTL overlapping the Chr1: 33 6 1 Mb region (Table 2). Two additional transcripts had significant eQTL mapping to Chr1 (Doc2g, Chr19; Trnd, Chr10) but these did not overlap and were not replicated in more than one other transcript (Table 2). Two transcripts, Kcnj3 and Tnni3 (Chr2 and Chr7, respectively), had overlapping significant eQTL on Chr2. For example, Tnni3 expression varied by 2.35fold between strains ( Figure 2D) and a trans eQTL was mapped to Chr2: 68.8 6 0.4 Mb ( Figure 2E). Of the remaining transcripts, 12 also had a suggestive eQTL overlapping the Chr2: 68.8 6 0.4 Mb region (Table 2). A second eQTL peak on Chr2: 77.5 6 1 Mb also reached significance for one gene (Mybpc3, Chr2) and was a suggestive eQTL for seven additional transcripts (Table 2). Using the Tnni3 eQTL plot as an example, this Chr2: 77.5 6 1 Mb region had a suggestive LRS ( Figure  2F). Since three significant eQTL intervals (Chr1: 33 6 1 Mb, Chr2: 68.8 6 0.4 Mb, and Chr2: 77.5 6 1 Mb) were mapped with suggestive LRS in multiple pulmonary myocardium transcripts, the genes within these regions were interrogated to identify putative upstream regulatory genes.
Rab23, Zfp451, Nostrin, Lrp2, and Ttn are candidate regulators of the pulmonary myocardium gene network The genes within the peak linkage regions were subsequently analyzed for known biological functions (GeneCards) (Safran et al. 2010) and expression in lung tissue [Fantom5 (Forrest et al. 2014;Severin et al. 2014); Table 3]. The Chr1: 33 6 1 Mb interval contained nine genes, with five of these (Prim2, Rab23, Bag2, Zfp451, and Dst) expressed in adult lung (Table 3). Both Rab23 and Zfp451 have known functions that could implicate them in regulating gene expression of the pulmonary myocardium network. RAB23, a GTPase of the RAB family involved in vesicular trafficking, also exerts downstream transcriptional effects by antagonizing the action of the transcription factor GLI1 (Sun et al. 2012). However, the function, localization, and mechanisms of action of RAB23 are still largely undefined. On the other hand, ZFP451 (zinc finger protein 451) is a promyelocytic leukemia nuclear body-associated transcriptional coregulator, and modulates transcription by regulating protein subnuclear localization (Karvonen et al. 2008).
The Chr2: 68.8 6 0.4 Mb interval contained nine genes, three of which were expressed at appreciable levels in adult lung (Lass6, Nostrin, and Lrp2) but only two had known roles in gene regulation (Nostrin) or cell signaling (Lrp2) ( Table 3). Nostrin has been shown to be required for vascular development and angiogenesis by acting as a molecular conduit for fibroblast growth factor receptor 1 (FGFR1)-mediated induction of Ras-related C3 botulinum toxin substrate 1 (RAC1) activity (Kovacevic et al. 2012). Interestingly, Nostrin has been shown in both mice (Kim et al. 2005) and humans (Wiesenthal et al. 2009) to bind to and negatively regulate transcription of its own promoter via a centrally located bZIP DNA binding motif (Bae et al. 2014).
However, transcriptional targets apart from Nostrin itself have not yet been found. LRP2 (low density lipoprotein receptor-related protein 2) is a critical component of the sonic hedgehog (SHH) signaling pathway (Christ et al. 2012). There is some evidence that LRP2 can shed its intracellular domain, which then shuttles to the nucleus where it may interact with as-yet-unidentified transcription factors (Li et al. 2008). However, a role for LRP2 has not yet been identified in either heart or lung. The Chr2: 77.5 6 1 Mb region contained six genes, only two of which were expressed in adult lung (Ttn and Ccdc141; Table 3). Of these, titin (Ttn) was a likely candidate modifier based on its role as a key integrator of myocyte signaling pathways (Linke 2008;Miller et al. 2004). Therefore, Rab23, Zfp451, Nostrin, Lrp2, and Ttn were analyzed further.
The transcription factors NKX2-5, SRF, MEF2A, TEAD1, and NR2F1 potentially cooperate to specify and/or maintain pulmonary myocardium identity Since NKX2-5 has previously been shown to be important for the maintenance of pulmonary myocardium cell identity, we investigated the possibility that genes within the pulmonary myocardium, and the candidate upstream regulators, might be direct targets of NKX2-5. To this end, we interrogated previously published NKX2-5 ChIP-seq data derived from two sources. He et al. (2011) used a doxycycline-inducible dual adenovirus system to express biotinylated NKX2-5 in the HL1 cardiomyocyte cell line, while Dupays et al. (2015) examined endogenous NKX2-5 binding in mouse hearts at embryonic day 11.5.
Of the 21 genes in our gene network, NKX2-5 binding sites were detected within 10 kb of seven genes (Rbm20, Corin, Myl7, Tnnt2, Ryr2, Fgf12, and Myh7) in the He et al. data set, two genes (Mybpc3 and Trdn) in the Dupays et al. data set, and three genes (Tbx20, Mybphl, and Atp2a2) in both data sets (Dupays et al. 2015;He et al. 2011) (Table  4). Collectively, this indicates that at least 50% of the genes in the pulmonary myocardium gene network are bound by NKX2-5 in heart tissue at some point during development. We also examined NKX2-5 enrichment at our candidate pulmonary myocardium control loci on Chr1 and Chr2 and found NKX2-5 enrichment at Rab23, Zfp451, and Nostrin in the He et al. data set. Although tissue-specific regulatory mechanisms may differ between the heart and the pulmonary myocardium, this is strong evidence to support the theory that NKX2-5 could directly regulate expression of the pulmonary myocardium gene network.
To further define transcriptional programs that might specify pulmonary myocardium identity, we performed transcription factor binding site enrichment analysis. We used the oPOSSUM-3 Web-based tool (Kwon et al. 2012) to identify overrepresented, conserved transcription factor binding sites in our set of 21 pulmonary myocardium genes and five candidate modifier genes. The pulmonary myocardium gene network was enriched in binding site sequences for SRF (serum response n and Dupays et al. (2015) data sets, we also investigated their enrichment in the genomic regions proximal to our genes of interest. Either one or both factors were enriched in 17 out of 21 genes (Table 4). Therefore, these transcription factors could form part of a regulatory network that directs or maintains the identity of the pulmonary myocardium.
We subsequently investigated the expression of Nkx2-5, Srf, Mef2a, Tead1, and Nr2f1 in the lung using the Fantom5 mouse promoterome database (data not shown). At the transcript level, the expression of all factors except Nkx2-5 could be detected in the lung (lung, adult: CNhs10474 ctss). The absence of Nkx2-5 expression in adult lung is an apparent paradox since this factor has previously been shown to be required for pulmonary myocardium identity (Millino et al. 2000;Mommersteeg et al. 2007). However, whether the lung sections sampled in the Fantom5 project actually contain pulmonary myocardium cannot be determined. In contrast, Nkx2-5 low-level expression in the BXD lung data set could be detected with a mean of 7.798 log 2 units [background is generally considered to be ,7 log 2 units (Geisert et al. 2009)]. Further, Nkx2-5 expression in the lung data set was significantly positively correlated with almost all the members of the pulmonary myocardium gene network (Table 4).
The pulmonary myocardium gene network and coexpressed genes are associated with cardiac phenotypes overlapping those of the stressed heart A limitation to single trait/transcript analysis is that it does not allow for combinatorial interactions with additional genetic loci, because each trait is individually regressed against every marker, i.e., only one locus is considered per trait (Michaelson et al. 2009). Therefore, all n 24 pulmonary myocardium gene network transcripts were merged into a single trait using PCA. A network of 100 probes that correlated (r . 0.553, P , 2.40e 205 ) with the synthetic PCA trait was identified.
Of the 100 probes, 76 had unique Entrez Gene IDs (Table S3 in File S2). These genes were considered to be part of a pulmonary myocardium gene coexpression network. This extended gene network contained a number of transcription factors (Tbx5, Hand2, Fhl2, and Gata4) and ion/calcium handling genes (Ank2, Cacna2d2, Scn5a, Kcnd2, Kcnj5, Myoz2, Casq2, Csrp3, and Gja3) of interest in the context of atrial fibrillation. We used this pulmonary myocardium gene coexpression network to investigate enriched GO terms, signaling pathways, and human phenotypes related to the coexpression of these genes in the lung using WebGestalt . The pulmonary myocardium gene coexpression network was significantly enriched with GO terms involving biological processes related to both striated and cardiac muscle function and development (Table 5). Enriched molecular functions included cytoskeletal protein binding and cation transmembrane transporter activity (Table 5), and enriched cellular components included contractile fiber and sarcomeric genes (Table 5). Collectively, these data support our approach for identifying genes that are likely to be expressed in the pulmonary myocardium, as our coexpressed genes are indeed enriched for cardiac muscle function.
Interrogation of related signaling pathways with the KEGG pathway tool revealed enrichment for "dilated cardiomyopathy," "cardiac muscle contraction," "hypertrophic cardiomyopathy," and "calcium signaling pathway" genes (Table 5). Similarly, analyses with the WikiPathway tool demonstrated enrichment for "striated muscle contraction" and "calcium regulation in the cardiac cell" (Table 5). In agreement with the GO and pathway results, the pulmonary myocardium gene coexpression network was enriched with higher order muscle phenotypes such as "abnormal cardiac muscle contractility," "abnormal cardiac tissue morphology," and "abnormal heart size" (Figure 3). Enrichment for pathways related to both abnormal contractility and calcium signaling/regulation led us to investigate the ion channels and regulatory proteins expressed in the extended gene network.
The expression of calcium handling networks in the pulmonary myocardium are dysregulated in atrial fibrillation We chose to investigate the potential for genes involved in ion transport and calcium handling to be dysregulated in atrial fibrillation. We selected 16 genes which are known to be involved in regulation of ion transport, almost all of which are also known susceptibility genes for arrhythmia (Table 6) Wang et al. 1995), and analyzed their expression in a publicly available mRNA microarray data set generated by Deshmukh et al. (2015). This study used RNA extracted from left atrial appendage tissue obtained from cardiac surgery patients who were divided into three categories: those with atrial fibrillation who either were or were not in sinus rhythm at the time tissue samples were collected, and those with no history of atrial fibrillation. The patient cohorts represent susceptibility to atrial fibrillation and persistent atrial fibrillation, respectively, when compared to control samples with no history of atrial fibrillation (Deshmukh et al. 2015). We compared our genes of interest (Table 6) to expression changes in this data set to determine if genes in our network are likely to be n Table 4 ChIP-seq enrichment of Nkx2-5, Srf, and Mef2a in genomic regions proximal to genes within the pulmonary myocardium gene network and candidate upstream regulators dysregulated in atrial fibrillation. Of these 16 genes, four were significantly downregulated (Cacna2d2, 0.5-fold; Kcnj5, 0.7-fold; Myoz2, 0.7fold; and Sln, 0.7-fold) and three were significantly upregulated (Atp2a2, 1.6-fold; Csrp3, 1.8-fold; and Gja3, 1.6-fold) in the persistent atrial fibrillation cohort compared to either the susceptibility or the control cohorts (Table 6).
The dysregulation of these genes in atrial tissue of patients with persistent atrial fibrillation, especially Ttn which we propose acts upstream of many of the genes in our network, highlights the need for similar gene expression studies to be performed with pulmonary myocardium tissue.
Gene network analysis effectively identifies markers expressed in the pulmonary myocardium from wholelung transcriptome data To confirm that our analyses have, in fact, identified genes expressed in the pulmonary myocardium, lung tissue samples were analyzed from four BALB/c mice using immunohistochemistry to detect the cellular localization of TNNI3 and ACTC1. ACTC1 was selected as its transcript was highly expressed in the lung data set (12.422 units) and was part of n our pulmonary myocardium gene coexpression network (Table S3 in File S2). Further, Actc1 has previously been shown to be a direct target of NKX2-5 (Chen et al. 1996). To differentiate between the pulmonary myocardium and typical smooth muscle, we also examined the localization of ACTA2. Cardiac ventricular tissue was used as a positive control for cardiomyocyte (TNNI3 and ACTC1) and smooth muscle (ACTA2) immunostaining ( Figure S2 in File S1). Distinct cell-type localization of TNNI3 and ACTA2 was observed in lung samples from four BALB/c mice (Figure 4; representative image from one mouse shown). A thin muscle coat staining positively for ACTA2 was situated beneath the columnar epithelium ( Figure 4A, black filled arrow) typical of bronchi and bronchioles in BALB/c mice (Rockx et al. 2007). Consistent with previous reports describing the pulmonary vein myocardium (Kracklauer et al. 2013;Mueller-Hoecker et al. 2008), a thick muscle coat positively staining for TNNI3 was situated beneath a monolayer of endothelium lining the vein lumen (Mueller-Hoecker et al. 2008) ( Figure 4B, white filled arrow). This myocardial layer also stained positively for ACTC1 ( Figure 4C, white filled arrow), and contained visible striations ( Figure 4C, insert with Ã ) previously reported to be evident in pulmonary vein myocardium (Kracklauer et al. 2013;Mueller-Hoecker et al. 2008). Weak staining of ACTC1 but not TNNI3 was detected in bronchial smooth muscle cells (Figure 4, B and C, black filled arrow). No positive staining was detected with an isotype control antibody ( Figure 4D).

DISCUSSION
Numerous studies now indicate that polymorphisms in genes encoding the cardiac transcription factors NKX2-5 and PITX2C increase susceptibility to atrial fibrillation (Huang et al. 2013;Qiu et al. 2014;Wang et al. 2014;Xie et al. 2013), suggesting the potential involvement of disrupted gene networks in this disorder. Both factors are required for development of the pulmonary myocardium (Mommersteeg et al. 2007;Ye et al. 2015), which is also thought to have a pathological role in atrial fibrillation (Chen et al. 1999). A large body of literature has also been amassed to show that variants in genes regulating calcium handling lead to arrhythmia (Table 6). However, how the pulmonary myocardium gains pacemaker activity in atrial fibrillation patients is still unknown. Further, the underlying molecular identity of the pulmonary myocardium, including which ion channels and transcription factors are expressed in this tissue, has not been described in detail. We have used a systems-genetics approach to examine a large network of coexpressed genes simultaneously without prior knowledge of their identity, and have attempted to shed light on the relationships between transcriptional regulation, ion/calcium handling, contractile phenotypes, and the gene networks involved. By correlating variation in gene expression and known BXD genotypes in a genome-wide eQTL analysis we identified three regulatory loci associated with the expression of a cardiac gene network in the pulmonary myocardium.
The best candidate gene in the Chr2: 68.8 6 0.4 Mb interval was Nostrin which is expressed in endothelial cells, heart, and highly n Table 6 The pulmonary myocardium gene coexpression network contains genes involved in ion transport that are associated with arrhythmia and dysregulated in atrial vascularized tissues including the lung (Zimmermann et al. 2002), and has clear potential as a regulatory molecule. A primary function of NOSTRIN is to modulate the activity of endothelial nitric oxide synthase (eNOS) by sequestering eNOS away from the plasma membrane (Zimmermann et al. 2002). In turn, this prevents the calcium-mediated release of nitric oxide (NO) (Zimmermann et al. 2002), a potent modifier of cardiac function, including regulation of contractility (Massion et al. 2003;Petroff et al. 2001). Alternative splicing of Nostrin produces a truncated isoform (NOSTRIN-b) which is primarily localized to the nucleus (Wiesenthal et al. 2009) and represses transcription from its own promoter (Kim et al. 2005;Wiesenthal et al. 2009). Therefore NOSTRIN-b could also influence transcription of other genes within the pulmonary myocardium gene network. A model can be envisaged where NOSTRIN-b alters expression of RBM20 (the linkage peak for the Rbm20 eQTL contains Nostrin), leading to altered splicing of Nostrin (as well as other key myocyte genes), to provide feedback from stretch-induced NO release and result in altered contractile properties. Interestingly, RBM20 has been shown to directly influence alternative splicing of Zfp451 and Ttn (Guo et al. 2013), both of which were identified as candidate upstream regulators in the other two eQTL peaks identified in this study. The eQTL interval located at Chr1: 33 6 1 Mb contained two possible candidates for regulating expression of the pulmonary myocardium gene network: Rab23 and Zfp451. Since a defined relationship exists between ZFP451 and cardiac function, we suggest that ZFP451 is the better candidate. ZFP451 has been shown to bind directly to SMAD3/4 which prevents recruitment of p300 to target promoters and results in the downregulation of transforming growth factor-b (TGFb) target genes in A549 lung epithelial cells (Feng et al. 2014). This is relevant since transgenic mice expressing a constitutively active form of TGFb 1 were shown to develop atrial but not ventricular fibrosis and had increased susceptibility to atrial fibrillation (Nakajima et al. 2000;Verheule et al. 2004). Upregulation of TGFb 1 is also known to induce myocardial hypertrophy (Parker et al. 1990) and is evident in animal models of heart failure (Dobaczewski et al. 2011). However, therapies which block TGFb would be expected to have adverse effects due to the pleiotropic nature of TGFb signaling (Dobaczewski et al. 2011). In future studies, it will be important to determine if ZFP451 expression in the pulmonary myocardium and heart of atrial fibrillation patients is altered in parallel with increased TGFb expression; if so, this transcription factor might be a viable therapeutic target.
A likely candidate underlying the Chr2: 77.5 6 1 Mb eQTL is the Ttn gene, encoding TTN, a giant muscle structural protein which acts as a myofibrillar backbone attached to the Z-disk, the thin filament, the thick filament, and the M-band (Linke 2008). TTN regulates the passive stretch of cardiomyocytes via the I-band region of the protein which is comprised of sequential Ig, PEVK, and N2B domains (Linke 2008). In the heart, alternative splicing of Ttn results in a shorter ventricular isoform (N2B) and a longer atrial isoform (N2BA) . These two isoforms are coexpressed at different ratios during development as a mechanism for fine-tuning the passive stiffness of the cardiomyocyte . TTN acts as a signaling hub from which direct and indirect interactions with multiple proteins can lead to diverse cellular responses, including changes in gene expression (Linke 2008).
The pulmonary myocardium gene coexpression network (Table S3 in File S2) contains a number of genes coding for TTN-interacting proteins, including cardiac myosin-binding protein C (Mybpc3; one of the original 24 transcripts), a-actinin (Actn2), obscurin (Obscn), Actc1, cardiac ankyrin protein [Ankrd1; also known as cardiac ankyrin repeat protein (CARP)], myomesin (Myom2), and four and a half LIM domains 2 (Fhl2). Of these TTN-interacting proteins, FHL2 interacts with the N2B domain (present in both isoforms), while CARP interacts with the N2A domain (present only in the longer N2BA isoform) (Linke 2008). Both FHL2 and CARP have been shown to act as transcription factors in response to specific stimuli. For example, CARP is upregulated at the myofibril and in the nucleus in response to stretch (Miller et al. 2003), while FHL2 is upregulated in response to RhoA signaling (Philippar et al. 2004). FHL2 has also been shown to interact with SRF to antagonize expression of smooth muscle genes in differentiating embryonic stem cells (Philippar et al. 2004). Similarly, CARP overexpression in murine neonatal cardiomyocytes represses expression of Anf, Actc1, Acta1, bMHC, ventricular myosin light chain 2 (Myl2), and cardiac troponin C (Tnnc) (Mikhailov and Torrado Figure 5 A model for pulmonary myocardium gene regulation in response to hemodynamic stress. Based on our gene network analysis we propose that hemodynamic changes trigger mechanosensory (stretch-stress) signaling from the sarcomere via TTN, and initiate NO signaling cascades which are transmitted to the myocyte nucleus via NOSTRIN. These signals initiate dynamic regulation of gene expression to enable finetuning of the contractile properties of the pulmonary myocardium. Gene expression changes may also be regulated at the post-transcriptional level by RBM20 which could feedback to regulate Nostrin, Ttn, and Zfp451. Collectively, these effects are predicted to result in altered contractile properties, ion handling, and hormonal signaling in the pulmonary myocardium to direct an effective response to the initial hemodynamic cue. 2008); genes which are associated with our pulmonary myocardium gene network. We therefore propose a model whereby mechanosensing by TTN could lead to altered gene expression of the pulmonary myocardium gene network via FHL2 or CARP ( Figure 5), which might exert opposing transcriptional effects depending on the TTN isoform expressed (as CARP signaling would only be activated by the N2BA isoform). Interestingly, Fhl2 had one of the highest fold-changes in expression in the left atria of patients with atrial fibrillation compared to healthy controls (Deshmukh et al. 2015), supporting a role for dysregulated TTN signaling in the pathobiology of atrial fibrillation.
Our eQTL analysis indicates that one of the genes potentially regulated by TTN is Rbm20, which in turn regulates Ttn isoform expression. RBM20 mediates alternative splicing of Ttn, resulting in increased expression of the shorter N2B isoform (increased stiffness) at the expense of the longer N2BA isoform (Beraldi et al. 2014;Guo et al. 2013;Wyles et al. 2016). RBM20 is also involved in alternative splicing of Mef2a, Znf451, and Trdn (Beraldi et al. 2014;Guo et al. 2013), and is associated with altered expression of Nkx2-5, Myl4, Myl7, Tnnt2, Mybpc3, and Actc1 (Beraldi et al. 2014). These changes in gene expression could be due to alternative splicing of MEF2A, which is likely to regulate the pulmonary myocardium gene network (Table 4). In knockout/-down studies, loss of RBM20 is associated with altered splicing and expression levels of .200 genes, and the authors conclude that RBM20 underpins structural and functional development of cardiomyocytes (Beraldi et al. 2014). RBM20 is also likely to be an important regulator of the pulmonary myocardium molecular phenotype, is associated with NKX2-5 expression, and could be an important therapeutic target in atrial fibrillation.
As a whole, our data imply that regulatory polymorphisms which alter expression levels of upstream regulators of this gene network (such as Nkx2-5, Mef2a, SRF, Nostrin, Znf451, and Ttn) might predispose to atrial fibrillation. This hypothesis is supported by the observation that expression of Ttn, one of our candidate upstream regulators, is increased in atrial fibrillation patients with persistent disease (Deshmukh et al. 2015) and a recent GWAS that has for the first time identified the TTN locus as a susceptibility locus for atrial fibrillation in humans (Christophersen et al. 2017). Further, genes within our gene network that are regulated by the Chr2 eQTL containing Ttn are also dysregulated in atrial fibrillation (Rbm20, Atp2a2, Fndc5, and Mybphl) (Deshmukh et al. 2015) and dysregulation of transcription factors and ion/calcium handling genes could drive the switch from working myocardium to the pacemaker phenotype observed in the pulmonary myocardium in atrial fibrillation.
Of particular interest in our gene network is the expression of Cacna2d2 which codes for the a-2/d subunit of the voltage-dependent calcium channel complex, a subunit that is highly and predominantly expressed in the sinoatrial and atrioventricular nodes (Marionneau et al. 2005). This gene is expressed at low levels in our gene network but is regulated from a suggestive eQTL on Chr1: 33 (one of the three loci from which our gene network is coregulated). Since expression of this subunit is associated with nodal tissues, upregulation of Cacna2d2, subsequent to dysregulation of our gene network, could predispose to gain of pacemaker activity. In support of this, Cacna2d2 expression is reported to be highly sensitive to dosage of the transcription factor TBX5 (Mori et al. 2006), another gene in our extended gene network. TBX5 is an essential factor for development of the cardiac conduction system including the atrioventricular bundle and bundle branches (Arnolds et al. 2012;Mori et al. 2006;Moskowitz et al. 2004). TBX5 has also been shown to interact and cooperate with both NKX2-5 (Hiroi et al. 2001;Puskaric et al. 2010) and GATA4 (Munshi et al. 2009) (also in our gene network) to drive commitment to a nodal cell phenotype (Stefanovic et al. 2014). Since NKX2-5 is known to suppress the nodal phenotype (Ye et al. 2015), a fine balance between TBX5 interactions with NKX2-5 vs. GATA4 may underlie the maintenance of the working myocardium in the pulmonary myocardium.
Taken together, a model is proposed whereby hemodynamic stress triggers mechanosensory (stretch-stress) signaling from the sarcomere and initiates NO signaling cascades from the adjacent pulmonary epithelia, which are transmitted back to the pulmonary myocardium myocyte nuclei ( Figure 5). This leads to dynamic regulation of gene expression, enabling fine-tuning of the contractile properties of the pulmonary myocardium, such as the rate and force of contraction, as well as influencing blood volume, ion transport, and hormone signaling. Our findings indicate that genetic variability at these loci can modify levels of gene expression in the pulmonary myocardium, and could modify venous return to the left atrium as well as pathological features such as hypertrophy, fibrosis, and initiation of ectopic beats. Hence the genes/loci identified in this study could be linked to atrial fibrillation in humans triggered by pulmonary events. Indeed, these loci regulate genes that in humans are known to be associated with various heart diseases (dilated cardiomyopathy, hypertrophic cardiomyopathy, and arrhythmia) and are therefore viable candidates for involvement in the pathobiology of atrial fibrillation.