Pneumoconiosis is the most serious occupational disease in China and its leading cause is occupational silica exposure. Pneumoconiosis takes several years to develop depending on the exposure level of silica. However, individual variation in the susceptibility to pneumoconiosis has been observed among the subjects with similar exposure. We conducted a genome-wide screening with 710 999 single nucleotide polymorphisms (SNPs) in a cohort of 400 coal workers (202 cases and 198 exposed controls) for pneumoconiosis susceptible loci. Seven promising variants were evaluated in an independent cohort of 568 coal workers (323 cases and 245 exposed controls), followed by a second replication on 463 iron ore workers (167 cases and 296 exposed controls). By pooling all of the genome-wide association studies and replication stages together, we found a genome-wide significant (P < 5.0 × 10−8) association for rs73329476 (P = 1.74 × 10−8, OR = 2.17, 95% CI = 1.66–2.85) and two additional replicated associations for rs4320486 (P < 0.05) and rs117626015 (P < 0.05) with combined P-values of 4.29 × 10−6 and 5.05 × 10−6, respectively. In addition, the risk allele T of rs73329476 was significantly associated with lower mRNA expression levels of carboxypeptidase M (CPM) in total cellular RNA from whole blood of 156 healthy individuals (P = 0.0252). The identified pneumoconiosis susceptibility loci may provide new insights into the pathogenesis of pneumoconiosis, and may also have some clinical utility for risk prediction for pneumoconiosis and high-risk population screening for workers with occupational silica exposure.

## INTRODUCTION

Pneumoconiosis is the most serious occupational disease in China. In 2012, 88.28% of the reported cases were attributed to pneumoconiosis and of these, over 50% are directly related to underground mining (1). Underground mining has mixed exposures, such as silica, metal or coal (carbon) dust generated during tunnel drilling, mining, roof bolting and transportation. Therefore, most cases of coal workers' pneumoconiosis (CWP) and iron ore workers' pneumoconiosis (IWP) are accompanied by silicosis. It has been reported that silica dust exposure is the leading cause of preventable pneumoconiosis (2,3). However, in China, >23 million workers are still exposed in their occupational settings (4). Recently, a large prospective cohort study showed that long-term silica dust exposure is associated with increased mortality, including deaths from respiratory diseases, cardiovascular diseases and lung cancer (4).

Occupational exposure of silica dust is the major etiological factor of pneumoconiosis; however, individual variation in the susceptibility to pneumoconiosis was observed among the workers with similar exposure at a given time (5). The past several years have witnessed an explosion of the associations between single nucleotide polymorphisms (SNP) and pneumoconiosis (69). However, for the majority of these studies using candidate gene approaches, the sample sizes were limited and the exposure information was not well controlled between cases and controls, which may lead to inconsistent results.

Recently, genome-wide association studies (GWAS) designed by multi-stage replications and high-throughput genotyping technology have been widely accepted as an effective way to identify susceptibility loci associated with various complex diseases or traits (10). To systematically explore genetic variants associated with susceptibility to pneumoconiosis, we conducted a GWAS of pneumoconiosis in Han Chinese workers by genotyping 900 015 SNPs in 205 pneumoconiosis cases and 203 controls with similar exposure levels using IlluminaOmniExpressZhonghua chips, and further evaluated suggestive associations involving pneumoconiosis risk by a two-stage replication with a total of 490 pneumoconiosis cases and 541 controls with similar exposure levels in the Han Chinese populations. Demographic and clinical information of the subjects are summarized in Table 1. Notably, an m: n matched case–control design for each stage based on cumulative dust exposure (CDE) and age at first dust exposure (AFDE) was carried out in our study (11,12). Detailed information is summarized in Materials and Methods.

Table 1.

Characteristics of the subjects in this study

Variables GWAS

Replication I

Replication II

Case Control Case Control Case Control
(n = 202) (n = 198) (n = 323) (n = 245) (n = 167) (n = 296)
AFDE (mean ± SD) 22.82 ± 4.23 22.21 ± 4.15 21.85 ± 3.97 21.42 ± 3.84 23.25 ± 4.82 22.89 ± 4.93
CDE (mg/m3-y) (mean ± SD) 76.04 ± 37.44 79.03 ± 43.08 97.33 ± 47.94 90.71 ± 32.96 112.93 ± 47.42 106.30 ± 38.50
Smoking status
Never 68 (33.66) 62 (31.31) 122 (37.77) 114 (46.53) 77 (46.11) 70 (23.65)
Ever 134 (66.34) 136 (68.69) 201 (62.23) 131 (53.47) 90 (53.89) 226 (76.35)
Stage
I 132 (65.34)  237 (73.37)  161 (96.41)
II 43 (21.29)  60 (18.58)  6 (3.59)
III 27 (13.37)  26 (8.05)  0 (0.00)
Variables GWAS

Replication I

Replication II

Case Control Case Control Case Control
(n = 202) (n = 198) (n = 323) (n = 245) (n = 167) (n = 296)
AFDE (mean ± SD) 22.82 ± 4.23 22.21 ± 4.15 21.85 ± 3.97 21.42 ± 3.84 23.25 ± 4.82 22.89 ± 4.93
CDE (mg/m3-y) (mean ± SD) 76.04 ± 37.44 79.03 ± 43.08 97.33 ± 47.94 90.71 ± 32.96 112.93 ± 47.42 106.30 ± 38.50
Smoking status
Never 68 (33.66) 62 (31.31) 122 (37.77) 114 (46.53) 77 (46.11) 70 (23.65)
Ever 134 (66.34) 136 (68.69) 201 (62.23) 131 (53.47) 90 (53.89) 226 (76.35)
Stage
I 132 (65.34)  237 (73.37)  161 (96.41)
II 43 (21.29)  60 (18.58)  6 (3.59)
III 27 (13.37)  26 (8.05)  0 (0.00)

AFDE, age at first silica dust exposure; SD, standard deviation; CDE, cumulative silica dust exposure.

## RESULTS

After quality control (QC) procedures to filter qualified samples and SNPs (Supplementary Material, Fig. S1), 202 cases and 198 controls with 710 999 autosomal SNPs were retained for subsequent association analyses. In the GWAS scan, P-values presented in the scatter plot of Figure 1 are derived from the additive model in the conditional logistic regression model (CLRM) (full set of genome-wide summary association statistics for GWAS data is provided in Supplementary Material). A quantile–quantile plot derived from the CLRM revealed a good match between the distributions of the observed P-values and those expected by chance (Supplementary Material, Fig. S2). A small genomic control inflation factor (λ) of 1.037 indicates minor contribution of population stratification to the association results. Details of the study design are shown in Supplementary Material, Figure S3.

Figure 1.

Genome-wide association results on pneumoconiosis in Han Chinese (202 cases and 198 controls). Scatter plot of P-values in the –log10 scale from the additive model on 710 999 SNPs.

Figure 1.

Genome-wide association results on pneumoconiosis in Han Chinese (202 cases and 198 controls). Scatter plot of P-values in the –log10 scale from the additive model on 710 999 SNPs.

We determined promising SNPs associated with risk of pneumoconiosis based on P < 1 × 10−4 in the additive model from the GWAS scan. Seven SNPs were selected to be further evaluated in the first stage of replication (Replication I) including 323 CWP cases and 245 exposed controls (Table 2 and Supplementary Material, Table S1). Although the SNP rs6660355 also met the significance threshold, it was excluded from further analysis as it was in strong linkage disequilibrium (LD) with the selected SNP rs9970018 (Supplementary Material, Table S2). The results of the 7 SNPs had P < 1 × 10−4 in the GWAS scan and the first-stage replications are shown in Table 2 and Supplementary Material, Table S1. Among these 7 SNPs, 3 SNPs at 7q31.1 (rs4320486), 12q12 (rs117626015) and 12q15 (rs73329476) were successfully replicated in the first-stage replication (P = 3.13 × 10−2, 2.39 × 10−2 and 7.93 × 10−3, respectively). These three loci were further validated in the second-stage replication (Replication II) including 167 IWP cases and 296 exposed controls (P = 3.05 × 10−2 for rs4320486, P = 3.12 × 10−2 for rs117626015 and P = 6.23 × 10−3 for rs73329476; Table 2).

Table 2.

Results of GWAS, replication and combined analysis

SNP Chr. Allelea Study Casesb Controlsb MAFc

Cases Controls
rs73329476 12q15 C/T GWAS 143/53/6 163/33/1 0.161 0.089 4.69 (2.40–9.15) 5.90 × 10−6
Replication I 231/77/9 177/47/2 0.150 0.113 1.75 (1.16–2.64) 7.93 × 10−3
Replication II 121/43/3 226/56/1 0.147 0.102 1.97 (1.21–3.21) 6.23 × 10−3
Combined     2.17 (1.66–2.85) 1.74 × 10−8
rs4320486 7q31.1 C/T GWAS 141/55/6 105/85/8 0.166 0.255 0.34 (0.20–0.58) 7.06 × 10−5
Replication I 190/109/14 129/86/17 0.219 0.259 0.69 (0.49–0.97) 3.13 × 10−2
Replication II 107/41/7 166/112/16 0.177 0.245 0.65 (0.44–0.96) 3.05 × 10−2
Combined     0.59 (0.47–0.74) 4.29 × 10−6
rs117626015 12q12 T/C GWAS 168/31/3 183/14/0 0.092 0.036 5.39 (2.38–12.21) 5.48 × 10−5
Replication I 278/41/2 220/20/1 0.070 0.046 1.91 (1.09–3.36) 2.39 × 10−2
Replication II 148/18/0 264/23/0 0.054 0.040 2.39 (1.08–5.27) 3.12 × 10−2
Combined     2.43 (1.66–3.56) 5.05 × 10−6
SNP Chr. Allelea Study Casesb Controlsb MAFc

Cases Controls
rs73329476 12q15 C/T GWAS 143/53/6 163/33/1 0.161 0.089 4.69 (2.40–9.15) 5.90 × 10−6
Replication I 231/77/9 177/47/2 0.150 0.113 1.75 (1.16–2.64) 7.93 × 10−3
Replication II 121/43/3 226/56/1 0.147 0.102 1.97 (1.21–3.21) 6.23 × 10−3
Combined     2.17 (1.66–2.85) 1.74 × 10−8
rs4320486 7q31.1 C/T GWAS 141/55/6 105/85/8 0.166 0.255 0.34 (0.20–0.58) 7.06 × 10−5
Replication I 190/109/14 129/86/17 0.219 0.259 0.69 (0.49–0.97) 3.13 × 10−2
Replication II 107/41/7 166/112/16 0.177 0.245 0.65 (0.44–0.96) 3.05 × 10−2
Combined     0.59 (0.47–0.74) 4.29 × 10−6
rs117626015 12q12 T/C GWAS 168/31/3 183/14/0 0.092 0.036 5.39 (2.38–12.21) 5.48 × 10−5
Replication I 278/41/2 220/20/1 0.070 0.046 1.91 (1.09–3.36) 2.39 × 10−2
Replication II 148/18/0 264/23/0 0.054 0.040 2.39 (1.08–5.27) 3.12 × 10−2
Combined     2.43 (1.66–3.56) 5.05 × 10−6

aMajor allele/minor allele.

bWild-type homozygote/heterozygote/variant homozygote.

cMAF, minor allele frequency.

dAdjusted for smoking status and eigenvectors where is appropriate.

By pooling all of the GWAS and replication stages together, we found a genome-wide significant (P < 5.0 × 10−8) association for rs73329476 (P = 1.74 × 10−8, OR = 2.17, 95% CI = 1.66–2.85) and two additional replicated associations for rs4320486 (P < 0.05) and rs117626015 (P < 0.05) with combined P-values of 4.29 × 10−6 and 5.05 × 10−6, respectively (Table 2). In addition, we observed similar association patterns of the three SNPs in subpopulations stratified by CDE, smoking status and stage classifications (Supplementary Material, Table S3).

In view of the modest or small effect of each individual locus, we put the above three SNPs together to assess the joint effect based on the replication stages. We found that the more unfavorable alleles the subjects carried, the higher risk of developing pneumoconiosis they have, suggesting a allele-dosage effect (Ptrend = 2.13 × 10−6; Table 3). Individuals with ‘2’ unfavorable alleles had a significantly 68% increased risk of developing pneumoconiosis (OR = 1.68, 95% CI: 1.17–2.41) compared with those having ‘0–1’ unfavorable allele. The elevated risk was more evident among subjects having ‘3’ or ‘4–5’ unfavorable alleles (OR = 2.45, 95% CI: 1.54–3.89 and OR = 5.66, 95% CI: 2.24–14.32, respectively). When using ‘2’ unfavorable alleles as a reference group, as expected, we found those with ‘0–1’ risk allele had a decreased risk (OR = 0.60, 95% CI: 0.42–0.86) while those with ‘3’ or ‘4–5’ unfavorable alleles showed higher risk (OR = 1.46, 95% CI: 0.98–2.18 and OR = 3.38, 95% CI: 1.38–8.27, respectively). In addition, for weighted risk score analysis, the comparison of the high-risk group to the low-risk and median-risk group showed a significantly increased risk of developing pneumoconiosis (OR = 3.58, 95% CI: 2.51–5.10 and OR = 1.82, 95% CI: 1.33–2.49, respectively) (Supplementary Material, Table S4).

Table 3.

Joint effect of SNPs at 12q15 (rs73329476), 7q31.1 (rs4320486) and 12q12 (rs117626015) on pneumoconiosis risk

Number of risk allelea Case (n = 466) n (%) Control (n = 499) n (%) Adjusted OR (95% CI)b P-valueb
0–1 122 (26.18) 171 (34.27) Reference –
223 (47.85) 239 (47.90) 1.68 (1.17–2.41) 5.19 × 10−3
100 (21.46) 80 (16.03) 2.45 (1.54–3.89) 1.48 × 10−4
4–5 21 (4.51) 9 (1.80) 5.66 (2.24–14.32) 2.52 × 10−4
Trend test    2.13 × 10−6
0–1 122 (26.18) 171 (34.27) 0.60 (0.42–0.86) 5.19 × 10−3
223 (47.85) 239 (47.90) Reference –
100 (21.46) 80 (16.03) 1.46 (0.98–2.18) 6.42 × 10−2
4–5 21 (4.51) 9 (1.80) 3.38 (1.38–8.27) 7.71 × 10−3
Number of risk allelea Case (n = 466) n (%) Control (n = 499) n (%) Adjusted OR (95% CI)b P-valueb
0–1 122 (26.18) 171 (34.27) Reference –
223 (47.85) 239 (47.90) 1.68 (1.17–2.41) 5.19 × 10−3
100 (21.46) 80 (16.03) 2.45 (1.54–3.89) 1.48 × 10−4
4–5 21 (4.51) 9 (1.80) 5.66 (2.24–14.32) 2.52 × 10−4
Trend test    2.13 × 10−6
0–1 122 (26.18) 171 (34.27) 0.60 (0.42–0.86) 5.19 × 10−3
223 (47.85) 239 (47.90) Reference –
100 (21.46) 80 (16.03) 1.46 (0.98–2.18) 6.42 × 10−2
4–5 21 (4.51) 9 (1.80) 3.38 (1.38–8.27) 7.71 × 10−3

ars73329476-T, rs4320486-C and rs117626015-C alleles were assumed as risk alleles.

bThe reference group was those with ‘0–1'or ‘2’ risk alleles. Adjusted for smoking status.

In combined data sets of all 692 CWP and IWP cases, the protective allele of rs4320486 was also significantly associated with slower disease progression (HR = 0.76, 95% CI = 0.65–0.89, P = 7.68 × 10−4, Supplementary Material, Table S5). The Kaplan–Meier plot shows that patients with the rs4320486 heterozygous CT genotype or the homozygous TT genotype had a median progression time (MPT) of 35.87 years, whereas the MPT was shortened to 32.98 years in those with the wild-type homozygote CC (Supplementary Material, Fig. S4). However, the other two SNPs (rs73329476 and rs117626015) were not significantly associated with disease progression (Supplementary Material, Table S5).

To characterize the functional relevance of the genome-wide significant SNP rs73329476, we further evaluated the relationship of this variant with the expression levels of three surrounding genes (CPM, MDM2 and CPSF6). We examined CPM, MDM2 and CPSF6 mRNA expression levels in total cellular RNA from whole blood of 156 healthy individuals using quantitative RT–PCR, and observed that the relative expression of CPM was significantly lower in subjects with T allele of rs73329476 (n = 31) when compared with those carrying CC genotype (n = 125) (P = 0.0252) (Fig. 2). However, non-significant results were observed for MDM2 (P = 0.5312) and CPSF6 (P = 0.3498).

Figure 2.

CPM relative expression levels according to rs73329476 genotypes in total cellular RNA from whole blood of 156 healthy individuals by quantitative RT–PCR. Lines indicate the median with inter-quartile range (25 and 75%). P-values were calculated by the Wilcoxon rank test.

Figure 2.

CPM relative expression levels according to rs73329476 genotypes in total cellular RNA from whole blood of 156 healthy individuals by quantitative RT–PCR. Lines indicate the median with inter-quartile range (25 and 75%). P-values were calculated by the Wilcoxon rank test.

Using imputation analyses based on the 1000 Genomes Project (Phase I integrated variant set, v3, October 2012), we tested the associations of the SNPs around all three regions and found 183 SNPs showing association (P < 1.0 × 10−3) with pneumoconiosis (Supplementary Material, Table S6). We also observed a series of significant signals surrounding the three marker SNPs in a 500 kb window (Fig. 3).

Figure 3.

Regional plot of the three marker SNPs with pneumoconiosis risk. Results (−log10 P) are shown for SNPs in the 500 kb flanking region on each side of the target SNP. The marker SNP was shown in purple diamond and the LD for each SNP is presented as colors representing r2 values. The genes within the region of interest are annotated below each plot, with arrows indicating the direction of transcription. (A) rs73329476; (B) rs4320486; (C) rs117626015.

Figure 3.

Regional plot of the three marker SNPs with pneumoconiosis risk. Results (−log10 P) are shown for SNPs in the 500 kb flanking region on each side of the target SNP. The marker SNP was shown in purple diamond and the LD for each SNP is presented as colors representing r2 values. The genes within the region of interest are annotated below each plot, with arrows indicating the direction of transcription. (A) rs73329476; (B) rs4320486; (C) rs117626015.

As presented in Supplementary Material, Table S7, three pathways were significant (FDR ≤ 0.1 and FWER ≤ 0.2) associated with pneumoconiosis risk in BioCarta database, including crebPathway (FDR = 0.092 and FWER= 0.080), mPRPathway (FDR = 0.067 and FWER= 0.110) and dreampathway (FDR = 0.070 and FWER= 0.167). However, in KEGG and GO databases, no pathways were identified. Association results of representative SNPs at P < 0.05 for genes in three identified pathways were shown in Supplementary Material, Table S8.

## DISCUSSION

CWP and IWP, which have some similar clinical patterns, result from the inhalation of mixed dust that usually contains free crystalline silica (13). Several factors such as dust concentration, years of exposure may affect the risk of developing pneumoconiosis (14). This study attempted, for the first time, to systematically evaluate genetic variants that were associated with pneumoconiosis susceptibility and identified a genome-wide significant (P < 5.0 × 10−8) association for rs73329476 and two additional replicated associations for rs4320486 (P < 0.05) and rs117626015 (P < 0.05) that were significantly associated with pneumoconiosis in Han Chinese. These findings may provide new insights into the pathogenesis of pneumoconiosis, and may also have some clinical utility for risk prediction for pneumoconiosis and high-risk population screening for workers with occupational silica exposure.

The SNP rs73329476 at 12q15 is between three genes: carboxypeptidase M (CPM, 180kb upstream), murine double minute 2 (MDM2, 265 kb downstream) and cleavage and polyadenylation specific factor 6 (CPSF6, 125 kb upstream) (Fig. 3A). CPM is an extracellular peptidase attached to the outer membrane by a glycosyl-phosphatidylinositol (GPI) anchor (15). CPM is highly expressed in the lung, especially on the plasma membrane of alveolar type I cells (16), and the activity of which was much higher in the lung than that found in other tissues (17). The high CPM activity may stem from its induction and release in acute lung disease (16). Importantly, CPM is associated with monocyte to macrophage differentiation and regarded as a marker of macrophage maturation, indicating it may be important in inflammatory responses (18,19). Correlation analysis results indicate that rs73329476 may be associated with the expression of CPM, subjects carrying the risk CT and TT genotypes at rs73329476 had lower mRNA expression level of CPM as compared with those carrying CC genotype (P = 0.0252) (Fig. 2). As expected, it is biologically plausible that the T risk allele at rs73329476 may result in lower expression of CPM, hence activating the crucial signaling pathways that drive lung fibrosis. However, these results are very preliminary and merit further experimental investigations.

MDM2 is a target gene of the transcription factor tumor protein p53. It has been reported that the attenuation of p53–Mdm2 conjugation may be involved in the epithelial cell apoptosis seen in idiopathic pulmonary fibrosis (20). Recent evidence showed that p53-regulated death pathway activated by oxidative stress play an important role in the pathobiology of asbestosis, which is pulmonary fibrosis resulting from asbestos exposure (21). Furthermore, the mRNA expression levels of p53 and mdm-2 were increased in mice after exposure to MeHgCl and lung tissue sections of MeHgCl-treated mice showed pathological fibrosis as compared with vehicle control (22). For CPSF6, the protein encoded by which is one subunit of a cleavage factor that required for 3’ RNA cleavage and polyadenylation processing (23). CPSF6 was reported to involve in the evasion of innate immune sensors and induction of a cell-autonomous innate immune response in primary human macrophages (24).

We also applied the publicly available data from GTEx (Genotype-Tissue Expression) eQTL (Expression quantitative trait loci) Browser, eQTL.Chicago.edu and Gene Expression Analysis to perform cis-eQTL analysis and evaluated the relationship between rs73329476 and the expression of nearby genes in a variety of cells/tissues. The GTEx project revealed that rs73329476 was a cis-eQTL of MDM2. Meanwhile, the SNP rs61928707, which is in strong LD with rs73329476 (r2 = 0.92), was also significantly associated with the mRNA expression level of MDM2. However, we did not found significant cis-eQTL for CPM and CPSF6 in any of the three open data sets, although our eQTL data showed that the T allele for rs73329476 was associated with lower mRNA expression of CPM. Since the GTEx project as well as the other two was performed mainly in Caucasians, the difference may due to genetic heterogeneity between different ethnics. Nevertheless, further eQTL analyses in lung tissues are required to evaluate the relationship between rs73329476 and CPM, MDM2 or CPSF6 expression levels because eQTL are usually tissue/cell specific (25).

The SNP rs4320486 at 7q31.1 is located 150 bp upstream of laminin, beta 1 (LAMB1) (Fig. 3B). Laminin-1 is critical for basement membrane organization and each of the laminin-1 α, β and γ chains is expressed by epithelial and mesenchymal cells during lung development (26). Furthermore, a significantly increased expression of the laminin and laminin-1 receptor was demonstrated in the lungs of mice with chronic allergen-induced airway remodeling, and the increased airway laminin expression may lead to lung fibrosis and deterioration in respiratory function (2628).

The SNP rs117626015 at 12q12 is located at intron 3 of PDZ domain containing ring finger 4 (PDZRN4) and 160 kb downstream of contactin-1(CNTN1) (Fig. 3C). PDZRN4 belongs to the LNX (PDZRN) family, which is implicated in cell fate determination through the inhibition of the Notch signaling (29). CNTN1, a member of the immunoglobulin superfamily, is a GPI—anchored neuronal membrane protein that functions as a cell adhesion molecule and acts as a ligand for the Notch signaling (30). Coincidently, both of PDZRN4 and CNTN1 play an important role in the Notch signaling, which may regulate lung development as well as pathological conditions of the lung, such as chronic obstructive pulmonary disease, fibrosis and lung cancer (31,32). In addition, the PDZRN4 region was speculated to be associated with multiple sclerosis (33).

Recently, a GWAS on Caucasian population was published using a discovery scan of 1616 patients with fibrotic idiopathic interstitial pneumonias (IIPs) and 4683 controls, with follow-up replication analyses in 876 cases and 1890 controls. They confirmed associations with variants in TERT at 5p15, MUC5B at 11p15 and the 3q26 region near TERC, and identified seven newly associated loci including FAM13A (4q22), DSP (6p24), OBFC1 (10q24), ATP11A (13q34), DPP9 (19p13), chromosomal regions 7q22 and 15q14-15 for IIPs (34). Considering that pneumoconiosis shares common characteristics with IIPs, we evaluated whether these IIPs related SNPs were also associated with pneumoconiosis. Among the 11 reported SNPs for IIPs, 7 SNPs were included in our original GWAS data set, and the marker SNPs for FAM13A (4q22) and ATP11A (13q34) were also associated with CWP at P < 0.05 (rs2609255 for FAM13A: OR = 1.40, P = 0.031; rs1278769 for ATP11A: OR = 0.68, P = 0.021) (Supplementary Material, Table S9). Moreover, we investigated whether the identified SNPs of pneumoconiosis risk also affect lung function in healthy populations by querying the results from a large meta-analysis GWAS of lung function (35). However, among the three SNPs, only rs4320486 was included in the data set, and the P-values were 0.73 and 0.82 for FEV1/FVC and FEV1, respectively.

Pathway-based approach that evaluates the cumulative contribution of the function-related genes may provide new insights into the biology of pneumoconiosis. In the current study, using BioCarta database, we identified three pathways (crebPathway, mPRPathway and dreampathway) that may play an important role in the development of pneumoconiosis in Han Chinese population. For the crebPathway (transcription factor CREB and its extracellular signals), it is reported that the cAMP response element-binding protein (CREB) signaling pathway, in particular the phosphorylation of CREB, is altered in pulmonary fibrosis (36). In addition, recent evidences suggest an important role for the CREB-binding protein signaling in the development of human fibrotic lung disease (3739). In the dreampathway (Repression of Pain Sensation by the Transcriptional Regulator DREAM), DREAM (DRE-antagonist modulator) is a Ca2+-regulated transcriptional repressor (40) and plays an important role in human Calcitonin gene expression (41). Calcitonin gene-related peptide (CGRP) is associated with fibrosis in the lung of male Wistar rats exposed to dusts (42), while degradation of CGRP may promote the transition from lung inflammation to fibrosis (43). For the mPRPathway (How Progesterone Initiates Oocyte Membrane), progesterone binds to its intracellular receptor, which induces activation of the MAPK signaling pathway (44), hence associated with the promotion of pulmonary fibrosis (4547).

In the current study, we observed a significant joint effect of the three loci on the pneumoconiosis risk when counting the number of risk alleles. Considering the effect of winner's curse bias in the discovery stage might overestimate the joint effect, the joint effect of the three SNPs was calculated only based on the replication stages. Besides, since using the most common group as a reference would make the results more meaningful (48), we also calculated joint effect of the three SNPs using ‘2’ as a reference group and found that, as expected, those with ‘0–1’ risk alleles had a decreased risk while those with more than two risk alleles showed higher risk. In addition, a weighted risk score was further calculated to analyze the joint effect of the three SNPs, which may be more accurate for joint effect estimation than that of risk allele counting especially for numerous loci with differential effects.

There are several strengths in the present study. First, this is the first relatively large GWAS of occupational disease, and successfully identified a genome-wide significant (P < 5.0 × 10−8) association and two additional replicated (P < 0.05) associations for pneumoconiosis susceptibility. Second, the pneumoconiosis cases and controls in each study stage were matched with occupational silica exposure levels, which are major determinants in pneumoconiosis pathogenesis. However, several limitations of our study also need to be addressed. First, considering the potential difference of pathogenesis of pneumoconiosis between coal workers and iron ore workers, homogeneity might be limited, since we used iron ore workers for the second-stage replication, although three loci identified from CWP were successfully validated by IWP. Second, the sample size may still limit the statistical power of the study. Further studies are warranted to validate and extend our findings, and to re-sequence the identified regions and to evaluate the potential functional significance of the loci.

## MATERIALS AND METHODS

### Study design

A three-stage case–control study was designed to evaluate the associations between genetic variants across human genome and the risk of pneumoconiosis. For the three stages, all the roentgenologic diagnosed pneumoconiosis cases were included from three independent cohorts of underground miners. The controls were also the underground miners with similar exposure level (e.g. age and years) as the cases.

With an effort to control the exposure level between cases and controls and improve the statistical power and study efficiency, we carried out an m: n matched case–control design for each stage based on CDE and AFDE (11). In principle, a case in each stage was generally matched to an exposed control with the most similar CDE and AFDE. In our study, we tried to retain all cases to evaluate the association between genetic variants and the susceptibility of developing pneumoconiosis. In particular, when there were no more desirable controls left in each cohort for a case to form new-matched risk set, then the cases were matched to a control in the matched risk set that was already generated. According to methods described by Hosmer and Lemeshow (12), there were no more than four cases or controls in each matched case–control set and the subjects were not overlapped between any two case–control sets.

### Study populations

The GWAS stage contains 205 CWP cases and 203 exposed controls, which were recruited from four coal mines of Datun Mining Business Group Co., Ltd from 2006 to 2012 in Jiangsu, China. The first replication (Replication I) includes independent 323 CWP cases and 245 exposed controls, which were recruited from five coal mines of Xuzhou Mining Business Group Co., Ltd from 2006 to 2012 in Jiangsu, China. The second replication (Replication II) samples are consisted of 167 IWP cases and 296 exposed controls, which were recruited from two iron mines from 1960 to 2003 in Hubei, China. All these subjects were underground coal or iron ore miners with occupational silica exposure and were unrelated ethnic Han Chinese who was used in previously published asso­ciation studies (6,7). In fact, of all the coal workers in our study, 82.13% worked at mine face and excavation face. During 2003 and 2008, we monitored the silica exposure level in Xuzhou Mining Business Group Co., Ltd and found that individuals at mine face and excavation face were exposed to silica in all monitoring sites (49). All subjects provided written informed consent, and the institutional review boards of each participating institution approved this collaborative study.

The subjects exposed to silica dust received routine health surveillance including physical examination and chest radiograph every 2 years, even after cessation of dust exposure. The cases were diagnosed as stage I, stage II or stage III pneumoconiosis according to the size, profusion and distribution range of opacities on chest X-ray by three national certified occupational physicians based on the China National Diagnostic Criteria for Pneumoconiosis (GBZ 70-2002), which is similar to the 1980 International Labor Organization on the classification of pneumoconiosis in the judgment of opacity profusion (50). In addition, it is worth mentioning that tuberculosis could contribute to some of the radiographic abnormalities that we have observed. In our study, the diagnoses of pneumoconiosis were made based on continuous chest X-rays for several years. Moreover, the chest X-rays of each case were further assessed by two independent physicians again and the cases with active tuberculosis were excluded. Therefore, the possibility of radiographic changes of tuberculosis can be minimized in this study.

We collected information about smoking habits through interviews. Those who had smoked an average of less than one cigarette per day and <1 year in their lifetime were defined as non-smokers; otherwise, they were considered as smokers.

### Occupational exposure data

The dust concentrations were measured monthly and calculated for each calendar year as recorded in historical industrial monitoring records. Meanwhile, duration of employment and job type for each subject was taken from personal employment records. We calculated the CDE for each underground employee from the starting date of dust-exposed work until the time of pneumoconiosis diagnosis or ending with dust-exposure as described previously (4).

### Quality control in GWAS

The GWAS scan was conducted in 205 cases and 203 exposed controls using IlluminaOmniExpressZhonghua chips with 900 015 SNPs. Before the association tests, we conducted a systematic QC procedure on the raw genotyping data to filter out both unqualified samples and SNPs (Supplementary Material, Fig. S1). One sample with overall genotype completion rates <95% was excluded. Additional four probable relatives were excluded based on pairwise identity by state according to their PI_HAT value in PLINK (all PI_HAT > 0.125). A set of 58 480 common autosomal SNPs [minor allele frequency (MAF) >0.25] with low LD (r2 < 0.05) were used to identify population outliers (n = 3) in the samples that had passed QC, using the founders of the HapMap trios of Yoruba in Ibadan (YRI, N = 90), Utah residents of Northern and Western European ancestry (CEU, n = 90), Han Chinese in Beijing (CHB, n = 45) and Japanese in Tokyo (JPT, N = 44) as the controls. SNPs were excluded if they met those criterion: (i) SNPs were not mapped on autosomal chromosomes; (ii) SNPs had a call rate <95%; (iii) SNPs had a minor MAF <0.05 in all screening samples and (iv) the genotype distributions of SNPs deviated from those expected by the Hardy–Weinberg equilibrium in all samples (P < 1 × 10−5). After the QC procedure, 202 cases and 198 controls with 710 999 SNPs were remained for further analysis.

### SNP selection and genotyping in the replication study

After genome-wide association analyses, we selected SNPs for the first-stage replication based on the following criteria: (i) SNPs had P < 1 × 10−4 for GWAS samples; (ii) they had clear genotyping clusters; (iii) only the SNP with the lowest P-value was selected when multiple SNPs were observed in a strong LD (r2 ≥ 0.8). As results, a total of eight SNPs satisfied the criteria (i) and (ii), and seven SNPs survived according to criterion (iii). Therefore, we genotyped these seven SNPs in the first-replication stage (Table 2 and Supplementary Material, Table S1). Although the SNP rs6660355 also met the significance threshold, it was excluded from further analysis as it was in strong LD with the selected SNP rs9970018 (Supplementary Material, Table S2). The SNPs showed significant associations with pneumoconiosis risk with P < 0.05 in the first-stage replication were selected for the second replication stage.

Genotyping analysis in the first-stage replication was performed using the iPLEX Sequenom MassARRAY platform (Sequenom, Inc). A TaqMan allelic discrimination assay (Applied Biosystems, Inc.) was used for the second-stage replication. A series of methods were used to control the quality of genotyping: (i) case and control samples were mixed on each plate; (ii) genotyping was performed without knowing case or control status; (iii) two water controls were used in each plate as blank controls; (iv) 10% of the samples were randomly selected to repeat the genotyping as blind duplicates, and the reproducibility was 100%.

### Statistical analysis

Genome-wide association analysis was performed using CLRM in an additive genetic model with covariate adjustment based on the matched study design. We used PLINK 1.07 for gen­eral statistical analysis. Population structure was evaluated by the principal component analysis in the software package EIGENSTRAT 4.2 (51). The Manhattan plot of −log10P was generated using R 2.15.1. Ungenotyped SNPs were imputed in the GWAS discovery samples using Shapeit 1.0 (phasing step), IMPUTE2 (imputation step) and the haplotype information from the 1000 Genomes Project (the Phase I integrated variant set across all 1092 individuals, v3, Oct 2012). The regional plot was generated using an online tool, LocusZoom 1.1. Principal component analysis revealed two significant (P < 0.05) eigenvectors when they were further included in the CLRM with other covariate of smoking status. For the GWAS scan, ORs and 95% CIs were estimated for each SNP from CLRM adjusted for smoking status and eigenvectors. For the combined association analyses of the GWAS scan and replication samples, we adjusted for smoking status. The χ2-based Cochran's Q statistic was calculated to test for heterogeneity between groups in a stratified analysis. Overall progression time for each case was defined as the duration from AFDE to the time of pneumoconiosis diagnosis. The hazard ratios (HRs) and 95% CIs were acquired from a multivariate Cox regression model adjusted for AFDE, CDE and smoking status. The Kaplan–Meier survivor function graph was plot using Stata version 9.2 (StataCorp LP, TX, USA). Joint effect was calculated based on a combination of the risk alleles. For any SNP, subjects carrying two risk alleles were coded as 2, while those carrying one or no risk alleles were coded as 1 or 0, respectively. The total risk alleles for each subject were further calculated by combining risk alleles for the three SNPs (T allele of rs73329476, C allele of rs4320486 and C allele of rs117626015). Joint effect of the three loci was estimated using CLRM with individuals carrying ‘0–1’ risk allele or ‘2’ risk alleles as reference groups. The weighted risk score for each individual was calculated using the equation of Y = SNPa × βa + SNPb × βb + SNPc × βc, where Y, a, b and c denote the risk score, rs73329476, rs4320486 and rs117626015, respectively, and βa, βb and βc are the reported ORs of the corresponding SNPs. Differences in expression between the rs73329476 CC, CT and TT genotypes were assessed by a Wilcoxon-signed rank test. The P-value was two sided, and the OR used in the manuscript was based on an additive model using CLRM, if not specified otherwise.

### Pathway analysis

We collected pathways from three public resources: Kyoto Encyclopedia of Genes and Genomes (KEGG), BioCarta and Gene Ontology (GO) database. Pathways containing genes between 10 and 200 were included in this study. SNPs were assigned to a gene if they located within 50 kb downstream or upstream of the gene. Of the total 710 999 genotyped SNPs, 499 250 SNPs were mapped into 23 487 genes within 50 kb upstream or downstream. Among them, 134 535 SNPs were ultimately assigned to 4592 genes within the pre-defined pathways (84 551 SNPs to 2919 genes in KEGG pathways, 299 10 SNPs to 985 genes in BioCarta pathways and 73 619 SNPs to 2491 genes in GO pathways). All genes were assigned to pathways containing 10–200 genes (178 pathways used in KEGG database, 191 pathways used in BioCarta database and 418 pathways used in GO database) in following pathway analysis. We then used the gene set enrichment analysis (GSEA) to perform pathway analysis. Briefly, three steps are used for pathway analysis in GSEA. First, individual–SNP association analysis was conducted to determine the effect for each SNP. Second, the representative SNP with the lowest P-value was mapped to each gene, and all genes were assigned to pre-defined biological pathways. Finally, the association between pneumoconiosis risk and each pathway was evaluated by the GenGen software using the weighted Kolmogorov–Smirnov-like running sum statistic (denoted by enrichment score, ES), which reflected the over-representation of a cluster of genes within this pathway at the top of the entire ranked list of genes in the genome. We randomly shuffled the case–control status for 10 000 times, and repeated these above steps to get the permuted pathway association results. Thus, the normalized ES after adjusted for different sizes of genes could be acquired by the permutation procedure. The significance level of pathways analysis was set to be FDR ≤ 0.1 and FWER ≤ 0.2.

### Quantitative reverse transcription polymerase chain reaction

Quantitative reverse transcription polymerase chain reaction (qRT–PCR) was performed to determine the mRNA expressions of CPM, MDM2 and CPSF6. RNAs from whole blood of 156 healthy individuals were isolated with the QIAamp RNA Blood Mini Kit (QIAGEN). We used TaqMan gene expression probes (Applied Biosystems, Inc.) to perform the qRT-PCR assay. All real-time PCR reactions, including no-template controls and real-time minus controls, were run by using the ABI7900 Real-Time PCR System (Applied Biosystems, Inc.) and performed in triplicate. β-Actin gene was used to normalize the expression levels. A relative expression was calculated using the equation 2−ΔCt (Ct, Cycle Threshold), in which $ΔCt=Ctgene−Ctβ−actin$.

## SUPPLEMENTARY MATERIAL

Supplementary Material is available at HMG online.

## FUNDING

This work was funded by the National Key Basic Research Program Grant (2011CB503805), the National Science and Technology Support Program (2011BAI09B02), the Major Program of the National Natural Science Foundation of China (81390543), the National Natural Science Foundation of China (81273044), the State Key Program of National Natural Science of China (81230067), Jiangsu Province Clinical Science and Technology Projects (BL2012008), the Priority Academic Program for the Development of Jiangsu Higher Education Institutions (Public Health and Preventive Medicine) and Collaborative Innovation Center For Cancer Personalized Medicine in Jiangsu Province.

## ACKNOWLEDGEMENTS

The authors wish to thank all the study participants, research staff and students who participated in this work.

Conflict of Interest statement. The authors had full access to the data and take full responsibility for its integrity. All authors have read and agree to the manuscript as written.

## REFERENCES

1
Ministry of Health of the People's Republic of China
National Report of Occupational Diseases in 2012
2012
Beijing
Ministry of Health of the People's Republic of China
2
Calvert
G.M.
Rice
F.L.
Boiano
J.M.
Sheehy
J.W.
Sanderson
W.T.
Occupational silica exposure and risk of various diseases: an analysis using death certificates from 27 states of the United States
Occup. Environ. Med.

2003
60
122
129
3
Castranova
V.
Vallyathan
V.
Silicosis and coal workers’ pneumoconiosis
Environ. Health. Perspect.

2000
108
Suppl. 4
675
684
4
Chen
W.
Liu
Y.
Wang
H.
Hnizdo
E.
Sun
Y.
Su
L.
Zhang
X.
Weng
S.
Bochmann
F.
Hearl
F.J.
et al.
Long-term exposure to silica dust and risk of total and cause-specific mortality in Chinese workers: a cohort study
PLoS. Med.

2012
9
e1001206
5
Yucesoy
B.
Luster
M.I.
Genetic susceptibility in pneumoconiosis
Toxicol. Lett.

2007
168
249
254
6
Ji
X.
Hou
Z.
Wang
T.
Jin
K.
Fan
J.
Luo
C.
Chen
M.
Han
R.
Ni
C.
Polymorphisms in inflammasome genes and risk of coal workers’ pneumoconiosis in a Chinese population
PLoS One

2012
7
e47949
7
Wang
M.
Wang
S.
Song
Z.
Ji
X.
Zhang
Z.
Zhou
J.
Ni
C.
Associations of IL-4, IL-4R, and IL-13 gene polymorphisms in coal workers’ pneumoconiosis in China: a case-control study
PLoS One

2011
6
e22624
8
R.
Mintz
M.
Rivas-Fuentes
S.
Jedlicka
A.
Lavergne
E.
Rodero
M.
Kauffmann
F.
C.
Kleeberger
S.R.
Polymorphisms in chemokine and chemokine receptor genes and the development of coal workers’ pneumoconiosis
Cytokine

2006
33
171
178
9
Ates
I.
Yucesoy
B.
Yucel
A.
Suzen
S.H.
Karakas
Y.
Karakaya
A.
Possible effect of gene polymorphisms on the release of TNFalpha and IL1 cytokines in coal workers’ pneumoconiosis
Exp. Toxicol. Pathol.

2011
63
175
179
10
Hardy
J.
Singleton
A.
Genomewide association studies and human disease
N. Engl. J. Med.

2009
360
1759
1768
11
Woodward
M.
Epidemiology study design and analysis
1999
Boca Raton
Chapman and Hall
264
266
1999
12
Hosmer
D.W.
Lemeshow
S.
Applied logistic regression
1989
New York
John Wiley & Sons
248
252
1989
13
Borm
P.J.
Tran
L.
From quartz hazard to quartz risk: the coal mines revisited
Ann. Occup. Hyg.

2002
46
25
32
14
Zhang
Q.
Huang
X.
Induction of ferritin and lipid peroxidation by coal samples with different prevalence of coal workers’ pneumoconiosis: role of iron in the coals
Am. J. Ind. Med.

2002
42
171
179
15
Denis
C.J.
Deiteren
K.
Hendriks
D.
Proost
P.
Lambeir
A.M.
Carboxypeptidase M in apoptosis, adipogenesis and cancer
Clin. Chim. Acta.

2013
415
306
316
16
Dragovic
T.
Schraufnagel
D.E.
Becker
R.P.
Sekosan
M.
Votta-Velis
E.G.
Erdos
E.G.
Carboxypeptidase M activity is increased in bronchoalveolar lavage in human lung disease
Am. J. Respir. Crit. Care. Med.

1995
152
760
764
17
Nagae
A.
Abe
M.
Becker
R.P.
Deddish
P.A.
Skidgel
R.A.
Erdos
E.G.
High concentration of carboxypeptidase M in lungs: presence of the enzyme in alveolar type I cells
Am. J. Respir. Cell. Mol. Biol.

1993
9
221
229
18
Fujiwara
N.
Ikeda
M.
Hirabayashi
S.
Mori
H.
Shirasawa
M.
Kansaku
A.
Sunamori
M.
Hata
Y.
Monoclonal antibody 7F9 recognizes rat protein homologous to human carboxypeptidase-M in developing and adult rat lung
Respirology

2007
12
54
62
19
Krause
S.W.
Rehli
M.
Andreesen
R.
Carboxypeptidase M as a marker of macrophage maturation
Immunol. Rev.

1998
161
119
127
20
Nakashima
N.
Kuwano
K.
Maeyama
T.
Hagimoto
N.
Yoshimi
M.
N.
M.
Nakanishi
Y.
The p53-Mdm2 association in epithelial cells in idiopathic pulmonary fibrosis and non-specific interstitial pneumonia
J. Clin. Pathol.

2005
58
583
589
21
Cheresh
P.
Kim
S.J.
Tulasiram
S.
Kamp
D.W.
Oxidative stress and pulmonary fibrosis
Biochim. Biophys. Acta

2013
1832
1028
1040
22
Lu
T.H.
Chen
C.H.
Lee
M.J.
Ho
T.J.
Leung
Y.M.
Hung
D.Z.
Yen
C.C.
He
T.Y.
Chen
Y.W.
Methylmercury chloride induces alveolar type II epithelial cell damage through an oxidative stress-related mitochondrial cell death pathway
Toxicol. Lett.

2010
194
70
78
23
Ruegsegger
U.
Blank
D.
Keller
W.
Human pre-mRNA cleavage factor Im is related to spliceosomal SR proteins and can be reconstituted in vitro from recombinant subunits
Mol. Cell

1998
1
243
253
24
Rasaiyaah
J.
Tan
C.P.
Fletcher
A.J.
Price
A.J.
Blondeau
C.
Hilditch
L.
Jacques
D.A.
Selwood
D.L.
James
L.C.
M.
et al.
HIV-1 evades innate immune recognition through specific cofactor recruitment
Nature

2013
503
402
405
25
GTEx Consortium
The genotype-tissue expression (GTEx) project
Nat. Genet.

2013
45
580
585
26
Elias
J.A.
Zhu
Z.
Chupp
G.
Homer
R.J.
Airway remodeling in asthma
J. Clin. Invest.

1999
104
1001
1006
27
Christie
P.E.
Jonas
M.
Tsai
C.H.
Chi
E.Y.
Henderson
W.R.
Jr.
Increase in laminin expression in allergic airway remodelling and decrease by dexamethasone
Eur. Respir. J.

2004
24
107
115
28
Dekkers
B.G.
Bos
I.S.
Halayko
A.J.
Zaagsma
J.
Meurs
H.
The laminin beta1-competing peptide YIGSR induces a hypercontractile, hypoproliferative airway smooth muscle phenotype in an animal model of allergic asthma
Respir. Res.

2010
11
170
29
Katoh
M.
Identification and characterization of human PDZRN4L gene and mouse Pdzrn4l gene in silico
Int. J. Mol. Med.

2004
13
923
927
30
Hu
Q.D.
Ang
B.T.
Karsak
M.
Hu
W.P.
Cui
X.Y.
Duka
T.
Takeda
Y.
Chia
W.
Sankar
N.
Ng
Y.K.
et al.
F3/contactin acts as a functional ligand for Notch during oligodendrocyte maturation
Cell

2003
115
163
175
31
Xu
K.
Moghal
N.
Egan
S.E.
Notch signaling in lung development and disease

2012
727
89
98
32
Su
J.L.
Yang
C.Y.
Shih
J.Y.
Wei
L.H.
Hsieh
C.Y.
Jeng
Y.M.
Wang
M.Y.
Yang
P.C.
Kuo
M.L.
Knockdown of contactin-1 expression suppresses invasion and metastasis of lung adenocarcinoma
Cancer. Res.

2006
66
2553
2561
33
Baranzini
S.E.
Wang
J.
Gibson
R.A.
Galwey
N.
Naegelin
Y.
Barkhof
F.
E.W.
Lindberg
R.L.
Uitdehaag
B.M.
Johnson
M.R.
et al.
Genome-wide association analysis of susceptibility and clinical phenotype in multiple sclerosis
Hum. Mol. Genet.

2009
18
767
778
34
Fingerlin
T.E.
Murphy
E.
Zhang
W.
Peljto
A.L.
Brown
K.K.
Steele
M.P.
Loyd
J.E.
Cosgrove
G.P.
Lynch
D.
Groshong
S.
et al.
Genome-wide association study identifies multiple susceptibility loci for pulmonary fibrosis
Nat. Genet.

2013
45
613
620
35
Soler Artigas
M.
Loth
D.W.
Wain
L.V.
Gharib
S.A.
Obeidat
M.
Tang
W.
Zhai
G.
Zhao
J.H.
Smith
A.V.
Huffman
J.E.
et al.
Genome-wide association and large-scale follow up identifies 16 new loci influencing lung function
Nat. Genet.

2011
43
1082
1090
36
Liu
X.
Sun
S.Q.
Ostrom
R.S.
Fibrotic lung fibroblasts show blunted inhibition by cAMP due to deficient cAMP response element-binding protein phosphorylation
J. Pharmacol. Exp. Ther.

2005
315
678
687
37
Henderson
W.R.
Jr.
Chi
E.Y.
Ye
X.
Nguyen
C.
Tien
Y.T.
Zhou
B.
Borok
Z.
Knight
D.A.
Kahn
M.
Inhibition of Wnt/beta-catenin/CREB binding protein (CBP) signaling reverses pulmonary fibrosis

2010
107
14309
14314
38
A.G.
Sokolow
A.
Shank
S.
Corey
D.
Myers
R.
Plafker
S.
Kelley
T.J.
Interaction with CREB binding protein modulates the activities of Nrf2 and NF-kappaB in cystic fibrosis airway epithelial cells
Am. J. Physiol. Lung. Cell. Mol. Physiol.

2012
302
L1221
L1231
39
Zhou
B.
Liu
Y.
Kahn
M.
Ann
D.K.
Han
A.
Wang
H.
Nguyen
C.
Flodby
P.
Zhong
Q.
Krishnaveni
M.S.
et al.
Interactions between beta-catenin and transforming growth factor-beta signaling pathways mediate epithelial-mesenchymal transition and are dependent on the transcriptional co-activator cAMP-response element-binding protein (CREB)-binding protein (CBP)
J. Biol. Chem.

2012
287
7026
7038
40
Carrion
A.M.
W.A.
Ledo
F.
Mellstrom
B.
Naranjo
J.R.
DREAM is a Ca2+-regulated transcriptional repressor
Nature

1999
398
80
84
41
Matsuda
M.
Yamamoto
T.A.
Hirata
M.
Ca2+-dependent regulation of calcitonin gene expression by the transcriptional repressor DREAM
Endocrinology

2006
147
4608
4617
42
Morimoto
Y.
Ogami
A.
Nagatomo
H.
Hirohashi
M.
Oyabu
T.
Kuroda
K.
Kawanami
Y.
Murakami
M.
Myojo
T.
Higashi
T.
et al.
Calcitonin gene-related peptide (CGRP) as hazard marker for lung injury induced by dusts
Inhal. Toxicol.

2007
19
283
289
43
Hartopo
A.B.
Emoto
N.
Vignon-Zellweger
N.
Suzuki
Y.
Yagi
K.
Nakayama
K.
Hirata
K.
Endothelin-converting enzyme-1 gene ablation attenuates pulmonary fibrosis via CGRP-cAMP/EPAC1 pathway
Am. J. Respir. Cell. Mol. Biol.

2013
48
465
476
44
Zhu
Y.
Rice
C.D.
Pang
Y.
Pace
M.
Thomas
P.
Cloning, expression, and characterization of a membrane progestin receptor and evidence it is an intermediary in meiotic maturation of fish oocytes

2003
100
2231
2236
45
Meng
Y.
Yu
C.H.
Li
W.
Li
T.
Luo
W.
Huang
S.
Wu
P.S.
Cai
S.X.
Li
X.
Angiotensin-converting enzyme 2/angiotensin-(1–7)/mas axis protects against lung fibrosis by inhibiting the MAPK/NF-κB pathway
Am. J. Respir. Cell. Mol. Biol.

2014
50
723
736
46
S.
Mukherjee
S.
Smith
M.
Das
S.K.
Activation of MAPK/AP-1 signaling pathway in lung injury induced by 2-chloroethyl ethyl sulfide, a mustard gas analog
Toxicol. Lett.

2008
181
112
117
47
S.K.
Schmidt
S.
Davidson
C.
Ikegami
M.
Wert
S.
Hardie
W.D.
MEK-ERK pathway modulation ameliorates pulmonary fibrosis associated with epidermal growth factor receptor activation
Am. J. Respir. Cell. Mol. Biol.

2012
46
380
388
48
Goddard
G.H.
Lewis
C.M.
Risk categorization for complex disorders according to genotype relative risk and precision in parameter estimates
Genet. Epidemiol.

2010
34
624
632
49
Song
Z.F.
Qian
H.Y.
Wang
S.S.
Jia
X.M.
Ye
Y.
Ni
C.H.
Analysis on the incidence of coal workers’ pneumoconiosis from 2003 to 2008 in a coal mining group
Zhonghua Lao Dong Wei Sheng Zhi Ye Bing Za Zhi

2011
29
56
58
50
Ministry of Health of the People's Republic of China
Diagnostic Criteria of Pneumoconiosis GBZ70–2002

2002
Xi'an, China
Law Press
1
10
51
Price
A.L.
Patterson
N.J.
Plenge
R.M.
Weinblatt
M.E.
N.A.
Reich
D.
Principal components analysis corrects for stratification in genome-wide association studies
Nat. Genet.

2006
38
904
909

## Author notes

These authors contributed equally to this work.
These authors contributed equally.