Colocalization analysis of 3′ UTR alternative polyadenylation quantitative trait loci reveals novel mechanisms underlying associations with lung function

Abstract

While many disease-associated single nucleotide polymorphisms (SNPs) are expression quantitative trait loci (eQTLs), a large proportion of genome-wide association study (GWAS) variants are of unknown function. Alternative polyadenylation (APA) plays an important role in posttranscriptional regulation by allowing genes to shorten or extend 3′ untranslated regions (UTRs). We hypothesized that genetic variants that affect APA in lung tissue may lend insight into the function of respiratory associated GWAS loci. We generated alternative polyadenylation (apa) QTLs using RNA sequencing and whole genome sequencing on 1241 subjects from the Lung Tissue Research Consortium (LTRC) as part of the NHLBI TOPMed project. We identified 56 179 APA sites corresponding to 13 582 unique genes after filtering out APA sites with low usage. We found that a total of 8831 APA sites were associated with at least one SNP with q-value < 0.05. The genomic distribution of lead APA SNPs indicated that the majority are intronic variants (33%), followed by downstream gene variants (26%), 3′ UTR variants (17%), and upstream gene variants (within 1 kb region upstream of transcriptional start site, 10%). APA sites in 193 genes colocalized with GWAS data for at least one phenotype. Genes containing the top APA sites associated with GWAS variants include membrane associated ring-CH-type finger 2 (MARCHF2), nectin cell adhesion molecule 2 (NECTIN2), and butyrophilin subfamily 3 member A2 (BTN3A2). Overall, these findings suggest that APA may be an important mechanism for genetic variants in lung function and chronic obstructive pulmonary disease (COPD).

alternative polyadenylation, GWAS, 3′ UTR length QTL, colocalization

Introduction

Chronic obstructive pulmonary disease (COPD) is a complex disease characterized by irreversible airflow obstruction on lung function testing. The leading environmental risk factor for COPD is cigarette smoking, however genetic factors have also been shown to play a role in disease susceptibility [1–4]. Through the use of genome-wide genotyping data, the COPDGene Study estimated the heritability of COPD as 37.7% in non-Hispanic white subjects and 37.9% in African Americans, while the heritability of lung function (forced expiratory volume in 1 s, FEV₁) is 38.4% in white subjects and 50.9% in African Americans [5]. Genome-wide association studies (GWAS) have been used to identify genetic variants associated with COPD and lung function [6–9]. A 2019 GWAS including 257 811 subjects from studies including the International COPD Genetics Consortium and UK Biobank identified 82 loci with genome wide significance [7], while a study of 400 102 individuals from UK Biobank and the SpiroMeta consortium identified 279 loci associated with lung function [8]. However, like most GWAS, the majority of identified variants lie in non-coding regions. To identify mechanisms connecting GWAS non-coding variants and human disease, the effect of genetic variants on gene expression (expression quantitative trait loci, or eQTL) have been successful in identifying target genes and functional mechanisms underlying disease associations. However, a large proportion of GWAS risk loci underlying complex traits such as COPD, remain unexplained [10].

Alternative polyadenylation (APA), which occurs in > 70% of human genes, plays an important role in posttranscriptional regulation by allowing genes to shorten or extend 3′ untranslated regions (UTRs) which may contain cis-regulatory elements including microRNA or RNA-binding protein binding sites. Thus APA can alter mRNA stability and translational efficiency, as well as protein localization [11]. The primary regulatory elements involved in APA are the polyadenylation signal (PAS), as well as other-cis elements located in the 3′ UTR, however other regulatory mechanisms underlying APA have been described [12].

There is growing evidence that genetic variants that alter 3′ UTR length, known as 3′ UTR APA quantitative trait loci (apaQTLs), can contribute to human disease, and several studies have demonstrated an association between genetic variants that alter 3′ UTR length and disease states [13–19]. However, to date, there has not been a comprehensive analysis of apaQTLs in lung tissue, or the role of apaQTLs in COPD and lung function, with prior analyses of apaQTLs including only 515 lung tissue samples from GTEx [20]. Here, we characterize apaQTLs in lung tissue samples from 1241 subjects and investigate their overlap with GWAS results for COPD case-control status and lung function.

Results

Discovery of alternative polyadenylation sites in lung tissue

Using RNA sequencing (RNAseq) data from 1241 lung tissue samples from subjects in the Lung Tissue Research Consortium (LTRC) (Table 1), we identified APA sites using DaPars2, which identifies and quantifies the usage of one proximal APA site in relation to the distal APA site for each annotated transcript meeting expression thresholds. We identified a total of 178 012 APA sites, and after filtering out transcripts with low coverage of the UTR (< 30 reads per sample), we were left with 56 179 transcripts corresponding to 13 583 genes. On average, there were 3.9 transcripts containing proximal APA sites per gene. The mean percentage distal usage index (PDUI) per sample was 0.57, and we found a roughly equal distribution of the distal and proximal APA sites, with 10 715 transcripts primarily utilizing the distal site (mean PDUI>0.9) and 9971 transcripts primarily utilizing the proximal site (mean PDUI<0.1). The mean distance between the proximal and distal APA site was 561 bases.

Table 1

Clinical characteristics of LTRC study individuals included in the analysis.

	Overall (n = 1241)
Gender, male, n (%)	52.4
Age, mean (SD)	63.3 (10.6)
Race, n (%)  White  Asian  Black  Hispanic  Other	1118 (90.1) 4 (0.3) 81 (6.5) 25 (2.0) 13 (1.0)
Current Smokers, n (% smokers)	73 (5.9)
Pack-Years Smoked, mean (SD)	31.8 (33.2)

Generation of APA-QTLs

For the purposes of this analysis, we define an apaTranscript as a transcript in which we were able to identify a proximal APA site whose usage in relation to the distal APA site was significantly associated with at least one SNP (i.e. had an apaQTL) (see Supplementary Fig. 1). We identified a total of 8831 apaTranscripts, corresponding to 4216 unique genes, with APA sites that were associated with at least one SNP with q-value < 0.05. We found that the majority of lead SNPs associated with apaTranscripts (apaVariants) were located in introns (33.0%), followed by downstream gene variants (25.8%), 3′ UTR variants (16.9%), and upstream gene variants (10.0%) (Table 2). The majority of apaTranscripts were protein coding (62.6%), followed by retained introns transcripts (lacking open reading frames) (18%) and non-coding processed transcripts (7.7%) (Table 3). Out of the 5896 unique apaVariants, 5083 (86%) were also associated with gene expression with q value < 0.05. We next investigated what proportion of the apaQTLs found in our data were previously identified in GTEx lung (n = 515). Data was available for 6545 significant apaTranscripts in GTEx, and we detected 4085 of these transcripts in our dataset. Out of these, we found that 1219 were affected by significant apaQTLs. Furthermore, 209 apaQTLs shared the same SNP-transcript pair. We found that 65 apaVariants disrupted known PAS motifs (Supplementary Table 1), however only 36 of these fully abrogated the PAS motif, while the other 29 variants created a new PAS motif. Pathway analysis of genes containing apaQTLs revealed overrepresentation of pathways including “Neutrophil degranulation”, “Major pathway of rRNA processing in the cytosol”, and “mRNA splicing-Major Pathway”(Supplementary Table 2).

Table 2

apaQTL variant annotation.

	Number	Percentage
intronvariant	1310	33.0
downstreamgene_variant	1023	25.8
3primeUTRvariant	669	16.9
upstreamgenevariant	396	10.0
noncodingtranscriptexonvariant	348	8.8
missensevariant	72	1.8
synonymousvariant	62	1.6
5 prime UTR variant	30	0.8
Sequence feature	16	0.4
Splice region variant & intron variant	14	0.4
5 prime UTR premature start codon gain variant	5	0.1
Splice region variant & non coding transcript exon variant	5	0.1
Splice acceptor variant & intron variant	4	0.1
Splice region variant	3	0.08
Splice donor variant & intron variant	2	0.05
Stop gained	2	0.05
Stop lost	2	0.05
Missense variant & splice region variant	1	0.03

Table 3

apaQTL transcript function.

	Number	Percentage
Protein coding	2483	62.6
Retained intron	715	18.0
Processed transcript	307	7.7
Nonsense mediated decay	202	5.1
antisense	120	3.0
lincRNA	97	2.4
IG V genea	13	0.3
Sense overlapping	5	0.1
Transcribed processed pseudogene	5	0.1
Processed pseudogene	4	0.1
Transcribed unprocessed pseudogene	4	0.1
TECb	3	0.08
IG C genea	2	0.05
IG V pseudogenec	1	0.03
Non stop decay	1	0.03
prime3 overlapping ncrna	1	0.03
Unprocessed pseudogene	1	0.03

^aImmunoglobulin (Ig) variable chain and T-cell receptor (TcR) genes annotated according to the international ImMunoGeneTics information system.

^bTEC (To be Experimentally Confirmed): Regions with EST clusters that have polyA features that could indicate the presence of protein coding genes. These require experimental validation, either by 5′ RACE or RT-PCR to extend the transcripts, or by confirming expression of the putatively-encoded peptide with specific antibodies.

^cInactivated immunoglobulin gene.

Colocalization analysis

To determine whether GWAS associations with COPD and lung function may be attributed to alternative polyadenylation, colocalization analysis was performed between apaQTL data and published GWAS data from COPD case control status [7] and lung function measures FEV₁, forced vital capacity (FVC), FEV₁/FVC ratio, and peak expiratory flow (PEF) [9]. We tested a total of 2104 unique genomic windows, including 569 windows for FEV₁/FVC, 656 windows for FEV₁, 638 windows for FVC, 467 windows for PEF and 21 windows for COPD case–control status containing both significant GWAS and apaQTLs. We identified APA sites in 193 genes that colocalized with at least one phenotype; specifically, APA sites in 63 genes colocalized with PEF, sites in 62 genes colocalized with FVC, sites in 27 genes colocalized with FEV₁, sites in 33 genes colocalized with FEV₁/FVC and sites in 8 genes colocalized with COPD (Supplementary Table 3). To select for APA sites with the strongest effect, we focused on apaQTLs with a beta score (slope from eQTL linear regression) of greater than 0.1, resulting in 27 loci for follow up (Table 4). Of these, four APA sites were located in the vicinity of the APA variant (within 60 bp), or within the 5′ UTR, and these were selected for further analysis.

Table 4

Top loci colocalized between GWAS and apaQTL data.

Locus	GWAS Phenotypes	Lead apaQTL SNP	apaQTL slopea	apaQTL   p-valuea	Gene Nameb	apaQTL transcript(s)	Predicted Proximal APA Site & Motifc	GWAS implicated gened
6:32645386:C:A 6:33228793:G:A	FEV1, PEF	6:32646995:G:A 6:32661737:T:A	0.41	< 1×10−317	HLA-DQA1	ENST00000343139.11 ENST00000443574.1	32 643 217 32 659 946	AGER
19:44892009:G:A 19:44922203:A:G	FEV, FVC	19:44878660:C:T	0.32	< 1×10−317	NECTIN2	ENST00000252485.8 ENST00000585601.1	44 878 606 44 878 613 No known motif
9:133274084:C:T	PEF	9:133273813:C:T	−0.32	1.21E-78	ABO	ENST00000611156.4 ENST00000651471.1 ENST00000453660.4 ENST00000647353.1	133 256 205 133 256 205 133 254 367 133 250 897	ABO
6:27607489:T:A	FVC	6:26378060:T:A 6:26376800:G:A 6:25384133:C:T	−0.31	1.78E-49	BTN3A2	ENST00000527422.5 ENST00000396948.5 ENST00000528222.1 ENST00000377708.7 ENST00000508906.6	26 375 666 26 375 666 26 368 267 26 375 948 26 375 948
19:3385231:G:A	FEV1	19:2936537:G:A	0.28	2.59E-10	ZNF77	ENST00000314531.5	2 934 664
13:78267917:G:A	FVC	13:79306066:C:A	−0.24	1.95E-42	RBM26-AS1	ENST00000456602.5 ENST00000656910.1 ENST00000606584.6 ENST00000654024.1	79 422 214 79 422 214 79 422 213 79 422 401
11:2001961:G:A	FEV1/FVC	11:1219991:G:T	0.25	1.45E-65	MUC5B	ENST00000529681.5	1 261 708
10:79946568:A:G	COPD	10:79747240:A:T 10:79753551:T:C 10:79758006:G:A	0.10	1.77E-103	NUTM2B-AS1	ENST00000665716.1 ENST00000610681.1 ENST00000483080.6 ENST00000661316.1 ENST00000658528.1	79 766 173 79 766 173 79 766 173 79 766 173 79 723 162	DYDC2   FAM213A   TMEM254
2:229904475:A:C	PEF	2:229948378:C:G	−0.20	0.000734	FBXO36	ENST00000283946.8	230 011 082
19:45790878:G:A	COPD	19:45721396:C:T	−0.16	7.36E-82	DMPK	ENST00000683086.1	45 770 270	DMPK   SNRPD2
6:33745319:C:T	FEV1	6:32658788:C:T	−0.15	3.49E-259	HLA-DQB1	ENST00000460185.1	32 660 068
19:3471605:A:G 19:4977049:T:G	PEF	19:4437453:A:G 19:4449290:T:C 19:4385662:A:G	−0.14	8.07E-09	CHAF1A	ENST00000585854.1 ENST00000587368.1 ENST00000301280.10	4 409 054 4 433 221 4 443 163
1:161545536:C:T	FEV1/FVC	1:161523491:A:G 1:161593265:G:A	0.04	2.76E-10	FCGR2C	ENST00000543859.5 ENST00000496692.6	161 599 806 161 589 809
6:28333322:A:G	PEF	6:29947986:C:T	−0.13	1.20E-17	AL645929.1	ENST00000420084.1	29 888 117
14:68978668:G:T 14:70045724:C:T	FVC	14:69487244:G:A	−0.13	2.06E-46	GALNT16	ENST00000448469.8	69 352 182
1:24555524:T:G	PEF	1:24555051:C:G	0.13	3.67E-10	NCMAP	ENST00000374392.3	24 606 527
15:41548040:A:G	PEF, FEV1/FVC	15:41554978:T:G	−0.13	6.46E-17	LTK	ENST00000453182.2 ENST00000563518.5	41 504 034 41 503 904
9:133270015:A:C 9:133271182:T:C	FVC, FEV1	9:133271182:T:C 9:133249012:G:A 9:133279786:G:A	−0.12	1.28E-57	ABO	ENST00000453660.4 ENST00000680600.1 ENST00000647353.1	133 254 367 133 250 524 133 250 897
19:8595399:C:T	FEV1/FVC	19:8438737:T:G	0.12	1.98E-112	MARCHF2	ENST00000381035.8 ENST00000215555.7	8 438 712 8 438 712 TTTTTT
11:68849554:A:G 11:68853972:C:T	FVC, FEV1	11:68857818:T:G	0.12	4.22E-09	IGHMBP2	ENST00000675684.1	68 929 301	IGHMBP2
16:19715459:A:G	PEF	16:19708624:C:A	0.11	7.37E-63	KNOP1	ENST00000565844.1	19 710 435
1:36322287:G:A	FEV1/FVC	1:36323945:A:G 1:36233260:C:T	0.10	1.39E-10	SH3D21	ENST00000672313.1 ENST00000480549.6 ENST00000373137.2 ENST00000672976.1	36 327 871 36 324 091 36 327 871 36 327 872 AATAAA
1:39594353:G:A	COPD	1:39757456:T:C	0.10	1.10E-11	PPIEL	ENST00000649325.1	39 522 354	PABPC4   OXCT2   MYCL
5:140937818:A:T	FVC	5:140735758:C:T 5:140650658:A:T	−0.10	2.66E-14	WDR55	ENST00000504897.2 ENST00000358337.10	140 672 529 140 669 690
16:30108980:C:T	PEF	16:30312489:C:T 16:30327682:T:C	−0.10	3.77E-46	NPIPB11	ENST00000524087.5 ENST00000614927.1	29 383 688 29 385 534
18:47098808:C:G	FVC	18:47108385:A:G	0.10	1.46E-54	AC012254.2	ENST00000602459.6	47 108 439
6:29176755:C:T 6:31598391:C:G	FVC FEV1	6:30726039:A:G	0.02	1.52E-25	HCG18	ENST00000681435.1 ENST00000680530.1 ENST00000681421.1 ENST00000656629.1 ENST00000657595.1 ENST00000670071.1 ENST00000659563.1 ENST00000665953.1 ENST00000656629.1	30 724 468 30 724 468 30 724 468 30 289 146 30 288 780 30 288 780 30 288 780 30 288 780 30 289 146	AGER

^aFor some loci, the lead apaQTL-variant was associated with APA sites in multiple transcripts, the p-value and slope for the association with the lowest P-value is reported here.

^bIf the lead apaQTL variant was associated with APA sites in multiple genes, only the gene containing associations with the lowest P-value are reported.

^cMotif information is only provided for loci where the variant is within the 5′ UTR.

^dGWAS implicated gene was ascertained from published GWAS findings for COPD case-control status [7] and lung function [9].

Colocalization of membrane associated ring-CH-type finger 2 (MARCHF2)

Variants in the 19p13.2 region have been associated with lung function in a GWAS [8]. Colocalization of GWAS results and apaQTL from a region flanking 2 MB of the lead variant revealed a shared signal between GWAS associations with FEV₁/FVC and apaQTLs, with a colocalization posterior probability (CPP) of 0.7. The top colocalized variant was rs1055554 (grch38, chr19:8438737_T/G) (Fig. 1a), and SuSie finemapping identified rs1055554 as one of five likely causal variants in the 95% credible set. The T allele of rs1055554 is associated with increased FEV₁/FVC (P = 5.8 × 10⁻⁶; ß = 0.012) as well as preferred usage of the proximal APA site (mean PDUI of TT subjects = 0.36). rs1055554 is located in the 3′ UTR of MARCHF2 and has a minor allele (G) frequency of 0.26 in LTRC. The G allele of rs1055554 disrupts the proximal PAS motif, changing the sequence of “TTTTTT” to ‘TTGTTT’, and resulting in increased usage of the distal APA site (mean PDUI of GG =0.58) (Fig. 1b and c). The T allele of rs1055554 is additionally associated with decreased gene expression of MARCHF2 (P = 4.1 × 10⁻¹⁴³; β = 4.0) (Fig. 1d).

Figure 1

apaQTLs for MARCHF2 colocalize with GWAS data for FEV₁/FVC. (a) Locus association plot for FEV₁/FVC GWAS and MARCHF2 apaQTL corresponding to ENST00000381035. The lead colocalized SNP, rs1055554, is highlighted with a diamond shape and used as the LD reference. (b) IGV coverage plot showing read pileup for the region spanning chr19:8,438,030-8,439,359 for 81 subjects from each genotype of rs1055554. (c) Boxplot of PDUI values for the ENST00000381035 transcript of MARCHF2. (d) Boxplot of RSEM gene expression values for MARCHF2.

Colocalization of NECTIN cell adhesion molecule 2 (NECTIN2)

Variants located in the 19q13.32 region have previously been associated with lung function [8] (rs3729640: FEV₁  P = 5.9 × 10⁻⁶, ß =0.013; FVC P = 7.7 × 10⁻⁶, ß = 0.013). Colocalization analysis of a 2 MB region flanking rs3729640 (grch38, chr19:44878660_C/T) revealed evidence of a shared signal between GWAS data for FEV₁ and FVC with apaQTL data with a CLPP of 0.8 and 0.6, respectively. The lead colocalized variant for both FEV₁ and FVC was rs3729640 (Fig. 2a); furthermore SuSiE finemapping revealed that the most likely causal variant for this apaQTL association is also rs3729640 (PIP = 1.00). The T allele of rs3729640 is associated with increased FEV₁ and FVC as well as decreased usage of the distal APA site. The mean PDUI of CC subjects is 0.43 compared to 0.97 in TT subjects (Fig. 2b and c). The variant did not result in a change in any known PAS motif, therefore the mechanism behind this association is unknown. We additionally found that the T allele is associated with a significant increase in NECTIN2 gene expression (P = 3.7 × 10⁻⁴²; β = 5.0) (Fig. 2d).

Figure 2

apaQTLs for NECTIN2 colocalize with GWAS data for FEV1 and FVC. (a) Locus association plots for FVC and FEV₁ GWAS and NECTIN2 apaQTL corresponding to ENST00000252485. (b) IGV coverage plot showing read pileup for the region spanning chr19:44,878,000-44,879,200 for 51 subjects from each genotype of rs3729640. (c) Boxplot of PDUI values for the ENST00000252485 transcript of NECTIN2. (d) Boxplot of RSEM gene expression values for NECTIN2.

Colocalization of butyrophilin subfamily 3 member A2 (BTN3A2)

A genome-wide significant association of variants in the 6p22.2 region with lung function has been identified in two large meta-analyses of lung function GWAS data [8, 9]. The minor allele (A) of the rs72841536 variant (grch38, chr6: 26378060:_T/A), is associated with decreased FVC (P = 2.98 × 10⁻²⁰, β = −0.032), however the definitive causal variant and mechanism has not been identified. We performed colocalization analysis of a 2 MB region in 6p22.2 with apaQTL data and identified a strong colocalization (PPA = 0.91) of signals from the two datasets (Fig. 3). The best colocalized SNP is rs72841536, which was also identified by SuSie as one of 55 variants in the 95% credible set. The A allele of rs72841536 was associated with decreased usage of the proximal site and therefore a shortened UTR. The mean PDUI for the TT subjects is 0.84 while the PDUI for AA is 0.79. The distance between the predicted proximal and distal APA sites is 222 bp. We further found that rs72841536 is an eQTL, with the A allele associated with decreased expression of BTN3A2 (P = 1.3 × 10⁻²⁷⁸, β = −57.7).

Figure 3

apaQTLs for BTN3A2 colocalize with GWAS data for FVC. (a) Locus association plot for FVC GWAS and BTN3A2 apaQTL corresponding to ENST00000527422 (BTN3A2-209). (b) IGV coverage plot showing read pileup for the region spanning chr6:26,375,333-26,376,136 (top) and chr6:26,374,161-26,378,479 for 17 subjects from each genotype of rs12116935. (c) Boxplot of PDUI values for the transcript ENST00000527422 (BTN3A2-204). (d) Boxplot of RSEM gene expression values for BTN3A2.

See supplementary text and Supplementary Fig. 3 for colocalization of SH3 domain containing 21(SH3D21).

Discussion

Although APA has been identified as an important aspect of gene regulation, affecting most human protein coding genes and contributing to many complex diseases and traits, the role of APA in COPD and lung function has not yet been characterized. In this analysis, which is the largest study of apaQTLs in human lung tissue, we identified a total of 8831 apaTranscripts from 4216 genes. Pathway analysis of genes containing significant apaQTLs identified pathways related to the inflammatory response such as “Neutrophil degranulation”, and “Interferon alpha/beta signaling”, as well as pathways related to gene expression and regulation such as “Major pathway of rRNA processing in the nucleolus and cytosol” and “mRNA splicing-Major Pathway”. Both inflammation and dysregulated expression as a result of cigarette smoke are known to play a role in the development and progression of lung disease [21, 22]. We identified substantial genetic colocalization between apaQTLs and genetic associations with COPD and lung function, with evidence of colocalization of APA sites in 193 genes with at least one GWAS phenotype. Out of the 27 top loci identified in our analysis, our study identified the same gene implicated by the GWAS for three loci, while it identified a different gene for four loci. In addition, we identified a target gene for 20 loci which previously had no known target. For three of these genes, MARCHF2, NECTIN2, and SH3D21, the apaVariant was located within 60 bp of the proximal APA site or within the 3′ UTR, and a variant in MARCHF2 was identified to modify known PAS motifs.

MARCHF2 is a member of the MARCH family of membrane bound E3 ubiquitin ligases, which function in the ubiquitination of target lysine protein to signal vesicular transport between membrane compartments [23, 24]. MARCHF2 has previously been shown to regulate the intracellular trafficking and secretion of α1-antitrypsin, a serine protease inhibitor in which mutations are associated with COPD, through the ubiquitination and degradation of the cargo receptor ERGIC3 [25]. In addition, MARCHF2 is involved in the ubiquitination and degradation of CFTR, the causal gene for cystic fibrosis, by associating with adaptor proteins CAL and STX6 [26], and can further play a role in CFTR-mediated autophagy [27]. Furthermore, MARCHF2 has been shown to negatively regulate the innate immune response by repressing NF-ΚB and type 1 IFN signaling pathways through the polyubiquitination and proteosomal degradation of IKBKG/NEMO [28]. We found that the lead colocalized variant for MARCHF2, rs1055554, results in the modification of a “TTTTTT” PAS signal to a “TTGTTT” sequence and decreases the usage of the proximal APA site resulting in a longer 3′ UTR. Previous studies have identified rs1055554 as an apaVariant for the same APA site in eleven cancer types [29]. Furthermore, a variant in LD with rs1055554 (rs59372292) has previously been associated with coronary artery disease [30], with the allele linked to a longer 5′ UTR (and decreased lung function) also associated with increased disease risk. The elongated 3′ UTR of MARCHF2 contains well conserved binding sites for several miRNAs including miR-383-5p.2, miR142-5p, miR-130-3 and miR-148-3p, which are not present in the shorter UTR sequence, suggesting that 3′ UTR length of MARCHF2 could trigger downstream effects, miRNA mediated transcript degradation. Furthermore we found that rs1055554 is located in a region containing H3K27Ac marks and overlapped by transcription factor binding sites (Supplementary Fig. 2a), further reinforcing the likely regulatory function of this variant.

NECTIN2 encodes a membrane glycoprotein which is one of the plasma membrane components of adherens junctions. It is mainly localized at the adherent junction in epithelial cells, and is ubiquitously expressed in a variety of cells including epithelial cells, endothelial cells, neurons and fibroblasts [31]. NECTIN2 can play a role in the initiation of intracellular signal transduction, tissue organization and cell polarity [32]. Genetic variants in NECTIN2 have been associated with many phenotypes including Alzheimer’s disease [33–35], coronary heart disease [36], multiple sclerosis [37], and cleft lip/palate [38]. The apaVariant, rs3729640, does not modify any known PAS motif near the predicted proximal APA site. This is consistent with our overall findings; we found that only a small number of apaVariants disrupted known PAS motifs. This is supported by two previous studies of apaQTLs [20, 39] which have shown that disruption of canonical signal motifs is not a major mechanism by which genetic variants influence APA. It has been proposed that other co-transcriptional mechanisms may play a role, such as competition between the spliceosome and polyadenylation factors mediated, for example, by spliceosomal RNA U1 [40] and RNAPII pausing [41]. Furthermore, RNA secondary structure may play a role, in addition to more complex mechanisms such as the recently described connection between DNA methylation, gene looping and APA regulation [42]. Investigation of regulatory and functional data (Supplementary Fig. 2b) revealed that rs3729640 is in a region containing H2K27Ac marks, contains transcription factor binding sites, and has a distal enhancer-like signature according to the ENCODE Registry of candidate cis-Regulatory Elements.

BTN3A2 encodes a transmembrane protein which is a member of the immunoglobulin superfamily [43]. The gene is located within the major histocompatibility class 1 locus, surrounded by the other family members on chromosome 6. Members of the BTN3A family are expressed in most types of immune cells including T cells, B cells, monocytes, dendritic cells and natural killer cells [44]. BTN3 molecules have been shown to inversely correlate with lymphocyte activity and are thus thought to contribute to regulation of the cellular immune response [45]. Mendelian randomization analysis has demonstrated a causal role between BTN3A2 variants and risk of COPD [46] as well as both childhood and adult onset asthma [47]. Therefore, BTN3A2 may contribute to lung function by altering the severity of the immune response. While rs72841536 has been found to be a significant eQTL in several tissues [48], it is located within the 3′ UTR of BTN3A2, and therefore the mechanism of action of this variant on gene expression was previously unknown. Here, we find that rs72841536 is associated with 3′ UTR length of BTN3A2, with the risk allele for decreased lung function (A), resulting in both a shortened UTR and significantly decreased gene expression. While we were not able to identify any specific regulatory elements within the lengthened UTR with may contribute to increased gene expression, the association between rs72841536 and BTN3A2 has previously been shown in five cancer cell lines [29], suggesting that this is a true and ubiquitous apaQTL. We also found that rs72841536 is located in a predicted transcription factor binding site for ARID5A, and is in a region containing K3K27Ac marks (Supplementary Fig. 2c).

The major strength of this study is the large sample size which allowed us to comprehensively characterize APA in lung tissue. While apaQTLs have previously been generated in several tissue types from GTEx, only 515 lung samples were available in that analysis [20], while our study contains over double that number (n = 1241). DaPars2 may have inflated Type 1 error in some scenarios, when compared to 3′ RNAseq [49], which may explain why only 1187 of the 6395 significant apaTranscripts identified in GTEx were also apaQTLs in our dataset. However, the large sample size of our analysis allowed us to be more stringent in selecting candidates for follow up, and we selected only apaQTLs with large effect sizes and low p-values. An additional explanation for the modest overlap between the apaQTLs identified in GTEx compared to our study may be that LTRC largely consists of current and former smokers, most with end stage lung disease. This may result in altered gene expression profiles compared to that in GTEx.

One potential weakness of this study is that we used standard RNAseq to estimate APA usage, which may have reduced accuracy compared to directly quantifying usage through 3′ RNAseq. In addition, using standard RNAseq, it is challenging to separate the contributions of co-transcriptional mechanisms resulting in APA site usage from post-transcriptional mechanisms such as isoform-specific decay [39]. Furthermore, our RNAseq was performed using whole lung tissue samples, which contain a variety of cell types. Therefore, it is impossible to disentangle which changes are occurring as a result of changes in transcriptional regulation from differences in inter-individual cell proportions. Additionally, it is important to note that when performing our colocalization analyses we included all GWAS variants with P < 1 × 10⁻⁵ instead of using the genome wide significance threshold of 1 × 10⁻⁸. The Bonferroni correction is known to be overly stringent in GWAS analysis due to the collinearity of variants from linkage disequilibrium. Therefore, a p-value of 1 × 10⁻⁵ is often used as a threshold for associations that are “suggestive” of significance, to identify subthreshold associations that can be validated using additional datasets. Finally, while we have identified changes in the 3′ UTR that are associated with genotype, using short read RNAseq alone, it is possible that some of the 3′ UTR changes observed are linked to mechanisms upstream of the 3′ UTR such as alternative promoter usage. Further work using long read RNAseq and protein characterization will provide more definitive evidence regarding the contribution of APA of identified genes to lung function.

In conclusion, there are many mechanisms by which genetic variation can contribute to gene regulation and subsequently modulate human disease. While many studies have demonstrated the role of gene expression and splicing in disease, few have characterized the role of 3′ UTR length on phenotype. Our study shows that apaQTLs contribute significantly to lung function, and identifies MARCHF2, NECTIN2 and BTN3A2 as novel candidate genes for lung function.

Materials and methods

Study population

All samples collected as part of the Lung Tissue Research Consortium were collected using a standardized protocol that has been described in the original study design [50]. A total of 1241 subjects had both RNAseq data and whole genome sequencing data (TOPMed Freeze 8) that passed QC and were thus included in the current analysis (Table 1). All samples were sequenced at the University of Washington Northwest Genomics Center (RNASeq) and Broad Genomics (whole genome sequencing), during phase 4 of the TOPMed program.

RNA extraction, quality control, sequencing and alignment

Details have been previously published [50]. Briefly, mRNA sequencing (RNAseq) was performed through the NHLBI TOPMed program at the University of Washington. Poly-A selection and cDNA synthesis was performed using the TruSeq Stranded mRNA kit (Illumina), and sequencing was performed on the NovaSeq6000 instrument. Sequences were aligned to GRCh38 using STAR (v2.6.1d) with the GENCODE release 29 reference. Gene-level expression quantification was performed using RSEM (v1.3.1).

Identification and quantification of alternative polyadenylation (APA)

DaPars2 [20] was used to calculate the percentage of distal poly(A) site usage index (PDUI), which is a measure of the expression of isoforms using the distal poly-A site relative to the total expression of isoforms using either the distal or proximal site. A PDUI closer to 1 indicates that there is an increased abundance of transcripts using the distal poly(A) site, resulting in a longer 3′ UTR. APA sites with less than 30 reads covering the UTR region in each subject were excluded, and transcripts with missing values in greater than 50% of subjects were removed. RSEM gene-level (v1.3.1) expression counts were generated using the Gencode GTF (v29). Genes with low expression were filtered out by keeping only genes with at least 1 count per million in at least 20% of subjects. Pathway analysis of genes containing APA sites was performed using Sigora [51].

3′ aQTL and eQTL analysis

PDUI ratios calculated using DaPars2 were quantile normalized, and RSEM gene levels counts were TMM normalized [52] prior to QTL testing. TensorQTL [53] was used to test for association between genotype of all SNPs within 1000 kb of a gene (cis-) and quantifications of PDUI or gene expression using linear models, adjusting for sex, batch, principal components of APA data and principal components of genetic ancestry. Calculation of principal components of genetic ancestry has been previously described [54]. A total of 8 792 206 variants (biallelic SNPs with MAF > 0.01) were tested for association with APA sites in 56 179 transcripts (corresponding too 13 582 unique genes) and expression of 16 266 genes. Results were annotated using ANNOVAR [55] with annotations derived from dbSNP build 150. Motif analysis has been previously described [56].

Colocalization analysis

Published GWAS data for COPD case-control status and lung function outcomes (FEV1, FVC, FEV1/FVC, PEF) [8, 9] were used for this analysis. Testing windows were generated by identifying all GWAS variants with P < 1 × 10⁻⁵ and calculating non-overlapping windows of 1 MB on either side of each SNP. Only windows containing apaVariants with FDR < 0.05 were tested. For each window, Bayesian colocalization tests were performed using the Moloc R package [57] to quantify the probability that the GWAS and apaQTL associations were due to a shared causal variant. Windows with a colocalization posterior probability (CPP) of greater than 0.6 were reported. Finemapping was performed using SuSie [58].

Acknowledgements

Molecular data for the Trans-Omics in Precision Medicine (TOPMed) program was supported by the National Heart, Lung and Blood Institute (NHLBI). Whole Genome Sequencing and RNASeq for “NHLBI TOPMed: The Lung Tissue Research Consortium (phs001662)” was performed at Northwest Genome Center (NWGC, HHSN268201600032I, RNASeq) and Broad Genomics (HHSN268201600034I, WGS). Core support including centralized genomic read mapping and genotype calling, along with variant quality metrics and filtering were provided by the TOPMed Informatics Research Center (3R01HL-117626-02S1; contract HHSN268201800002I). Core support including phenotype harmonization, data management, sample-identity QC, and general program coordination were provided by the TOPMed Data Coordinating Center (R01HL-120393; U01HL-120393; contract HHSN268201800001I). We gratefully acknowledge the studies and participants who provided biological samples and data for TOPMed.

Author contributions

Conceptualization: A.S., P.J.C., C.P.H.; Data Curation: M.H.C., P.J.C., C.P.H.; Formal Analysis: A.S.; W.K., Z.X., R.P.C.; Writing—Original Draft Preparation: A.S., C.P.H.; Writing—Review & Editing: A.S., W.K., Z.X., R.P.C., M.H.C., A.L., P.J.C., C.P.H.

Conflict of interest statement: PJC received grant support from Bayer and Sanofi and consulting fees from Verona Pharmaceuticals. CPH has received research grants from Bayer, Boehringer-Ingelheim, and Vertex, and consulting fees from Chiesi, Sanofi, and Takeda. MHC has received grant funding from Bayer.

Funding

This work was funded by K01HL157613, R01HL157879, P01HL114501, X01HL139404, R01HL124233 and R01HL126596. MHC was supported by R01HL153248, R01HL149861, R01 HL111527HL135142, and NIGMS R35 GM140844. AL was supported by R01HL147148.

Data availability

LTRC genotyping data are available on dbGaP with accession number phs001662.v1.p1.

LTRC RNA-seq data are available through TOPMed, https://topmed.nhlbi.nih.gov

Results of apaQTL analysis are available on dbGap (phs001662.v2.p1).