Collaborative Cross-MexTAg breeding (RA/3/300/106; RA/3/300/107) and experimental (RA/3/100/1,408; RA3/100/1730) protocols were approved by the University of Western Australia animal ethics committee (UWA AEC) in accordance with the Australian code for the use of animals in medical research (Council, 2013). Male CC mice were generously provided by Geniad Pty Ltd. from its colonies at the Animal Resource Centre (Perth, Western Australia). Additional male CC strains (denoted by CCXXX format, Supplementary Table S1) were obtained from the University of North Carolina CC colony. Female parental homozygous 266-MexTAg mice (C57Bl/6, H2-K^b) were bred and housed at the University of Western Australia (UWA) Biomedical Research Facility (Perth, Western Australia).

All CCMT progeny and 266-MexTAg heterozygous control mice received a total of 6 mg asbestos administered as two intraperitoneal (i.p.) injections of 3 mg sterile asbestos (IUCC reference sample; Wittenoom Gorge crocidolite) suspended in 0.5 mL PBS at weeks 0 and 4 and survival calculated from the day of first injection as per published protocols (Robinson et al., 2006; Robinson et al., 2011). All mice were humanely euthanased in accordance with the UWA AEC approved animal care services standard operating procedure “201–8 Euthanasia of the Mouse”, primarily via methoxyflurane inhalation overdose (Medical Developments International, Australia. Penthrox^® methoxyflurane, 99.9% v/v) followed by cervical dislocation.

Tissue samples including macroscopic tumour, spleen, kidney, liver, and diaphragm were collected for histological and genetic analysis (RNAlater, QIAEN). Cell lines were generated from peritoneal ascitic fluid and macroscopic tumour when possible (see below). Tissues were fixed in 4% formaldehyde (Amber Scientific Pty Ltd., Perth, Western Australia) for 24–48 h, preserved in 70% ethanol prior to embedding in paraffin blocks (Surgipath Paraplast paraffin Leica Biosystems, Australia). Five micrometre (5 µm) sections were cut and stained with haematoxylin and eosin (H&E). Histopathological analysis was performed via bright field microscopy (Nikon Eclipse E200 microscope (Minato City, Tokyo, Japan), with selected sections scanned using Leica (Aperio) Scanscope Digital Slide Scanner. Ten (10) histological features were reviewed with each feature defined as representative of either benign or malignant disease. Additionally, the histological subtype of each mouse was also noted (epithelioid, sarcomatoid or biphasic).

Ascites was collected under aseptic conditions immediately post euthanasia and placed in tissue culture flasks with at least 2X volume of DMEM (Sigma) supplemented with 15% NCS (HyClone Cytiva), under 5% CO₂ and 95% humidity. Likewise, when possible, solid tumours were dissociated using scalpels and small fragments cultured as described above. Cells were passaged as required and frozen stocks stored in liquid nitrogen.

Kaplan–Meier survival curves were analysed by log rank test (Mantel-Cox) with >95% confidence intervals (CI). Correlations were analysed by Pearson’s test for correlation with 95% CI. Statistical analyses were performed using Graph Pad Prism V8.4.2 (Graph Pad Software Inc., United States of America). One-way ANOVA test for variance was used to analyse data from three or more groups. The non-parametric, unpaired, two-tailed t-test was used to compare data from two test groups. Frequencies of histological features between groups was compared using Chi-Square frequency table analysis and logistic regression. Endpoints for survival analyses included disease latency; disease progression; overall survival (OS) and cancer-specific survival (CSS), defined as time from first injection until death from mesothelioma. Disease latency and disease progression were analysed using linear regression. OS and CSS were analysed using Cox proportional hazard modelling. Multivariate analysis was performed using stepwise backward selection with all significant variables in univariate analysis included in the starting model with a p < 0.05 required to remain in the model. Variables were considered significant at α < 0.05. Features with less than 5% of a feature either present or not present across the cohort were excluded to reduce type 1 error. For histological analyses all statistical analyses were performed in SAS v9.4 (IBM; Cary, NC, United States).

Candidate modifier genes were identified using the GeneMiner bioinformatic suite (https://www.sysgen.org/Geniad2/) that combines the HAPPY (Mott et al., 2000; Durrant et al., 2011) and DO-QTL programs to define the founder haplotypes associated with respective biological traits as previously described (Ram et al., 2014). Genome wide scans were performed to define chromosomal locations of peak SNPs associated with each respective phenotype, such as overall survival, disease progression, latency, and ascites volume. Multinomial logistic regression models were fitted for each trait at each locus and ANOVA chi-square tests used to estimate the p-value of association. A false discovery rate (FDR) of p < 0.001 was used to define significant genome-wide linkage (Benjamini et al., 2001). The founder strain(s) contributing to each trait were determined by deriving coefficients (log odds ratio) of the fit from the multinominal regression model, as implemented in DO-QTL.

Chromosomal regions containing peak eQTL were interrogated using either the Mouse Genome Informatics database (https://www.informatics.jax.org/) or Sanger Sequencing database (via the GeneMiner informatics portal; https://www.sysgen.org/Geniad2/) to identify candidate modifier genes.

RNAseq datasets from the Bueno (Bueno et al., 2016) and The Cancer Genome Atlas (TCGA-MESO) human mesothelioma cohorts were used to investigate the association between the expression of candidate modifier genes in tumours and age of mesothelioma patients at the time of diagnosis/surgery. To utilize RNAseq datasets from the Bueno cohort, we imported RNAseq datasets in raw FASQ file format from the European genome-phenome archive (EGA, accession code: EGAS00001001563) and aligned them with Kallisto (v0.46.1) against the human reference genome (GRCh38). The “TCGAbiolinks” R package was used to retrieve RNAseq data from TCGA datasets in the STAR aligned raw gene count format from the Genomic Data Commons (GDC) portal. We converted CCMT candidate modifier gene symbols into homologous human gene symbols and performed univariate cox regression analysis using the “coxph” function from the “survival” R package to identify genes associated with age at the time of surgery, and age at the time of diagnosis in the Bueno and TCGA cohorts respectively. Genes with p-value < 0.05 were selected for Kaplan Meier survival analysis using the “survminer” R package. All analyses were conducted using Rstudio (version 4.1.0).

The generation of MexTAg mice has been described previously (Robinson et al., 2006; Robinson et al., 2011). To generate CCMT mice, female homozygous 266-MexTAg mice, that contain two copies of the SV40 TAg transgene, were crossed with male mice from 71 distinct CC strains in a sequential, staggered, batch-breeding protocol (Figure 1A). All CC-MexTAg progeny carried a single copy of the SV40 TAg transgene but differ genetically based on the genotype of the parental CC strain. All attempts were made to match CCMT experimental groups for age and gender balance prior to asbestos exposure. Characteristics of CCMT groups are described in Supplementary Table S1.

Seventy-one (71) groups of CCMT mice and a single group of 266-MexTAg heterozygous controls (266Het, containing a single copy of TAg; cohort median n = 37/group, range 18–45/group, 2,565 total mice) received intraperitoneal asbestos injections as previously described (Robinson et al., 2006; Robinson et al., 2011) and were monitored for ARD development (Figure 1B). Asbestos-exposed mice were assessed for traits including overall survival (time from first asbestos exposure to cull); ARD latency (time from first asbestos exposure to first signs of disease; FSD); ARD progression (time from FSD to euthanasia) and mean ascites volume (Figure 1C, Supplementary Table S2). Mice were culled at predefined welfare endpoints (namely, presence of ascites or loss of condition), or at the 18 months (548 days) experimental endpoint, whichever occurred first. Tissue samples were collected upon euthanasia and analyses of phenotypic parameters and histological sections performed (Figure 1D).

To assess the impact of host genetics on ARD development, overall survival (time from first asbestos exposure to cull) was analysed (Figures 2A, B). Significant variation in overall survival was observed between different asbestos exposed groups (p > 0.001), with a 3.75-fold difference observed between groups with the shortest, compared to those with the longest median survival (Figure 2C; red and green bars respectively). ARD incidence was consistent across most CCMT groups (cohort median 87%, range 50%–100%, with only progeny of four CC strains (SEH 50%; LIV 62.5%; ZIE2 66.7% and LUF 67.5%) having an ARD incidence more than two-standard deviations lower than the cohort mean (87% ± 18.4% (2SD), Supplementary Table S1). Interestingly, heterozygous parental 266-MexTAg control mice, which contain a single copy of the SV40 TAg transgene, displayed similar median survival (536.5 days) and incidence (93.8%) to the longest surviving CCMT groups (groups CC059 to DAVIS green bars Figure 2C; 548 days; 87.5% ± 5% respectively, Supplementary Table S1); suggesting that a) the B6 strain (one of the CC founders and host of the MexTAg transgene) confers more protection than the other founders; and b) such modifiers are not linked to the single copy of TAg as linked genes do not mask the influence of host genetics on ARD development.

We next assessed variation in ARD latency (time from first asbestos exposure to FSD) and progression (time from FSD to Cull; Figures 3A, B). The variation in ARD latency was identical to that observed for overall survival (Pearson’s correlation coefficient r = 0.9988), while an inverse correlation was observed between overall survival and disease progression (r = −0.1291, Figures 3D, E). The implications of the strong positive correlation between overall survival and disease latency are two-fold: firstly, it suggests that overall survival and disease latency are likely the same trait and secondly, it indicates that the variation observed in overall survival occurred during the latency period and not once disease was established (i.e., during ARD progression).

We further assessed the relationship between each of the different ARD phenotypic traits (Figures 3C–I). Consistent with the strong positive correlation between latency and survival, both survival and latency demonstrated a positive correlation with mean ascites volume (Figures 3F, G). However, little correlation was observed between ARD progression and all other phenotypes (Figures 3H–J). Taken together, these data demonstrate that asbestos-exposed mice developed ARD during a long latency period, consistent with the observed long period for human mesothelioma development.

We next sought to characterise the histological features associated with asbestos exposed CCMT mice. Spleen, kidney, liver, and diaphragm were harvested from all asbestos-exposed mice at predefined disease-associated or experimental endpoints and were assessed for signs of histological disease. Bulk tumour tissue was also assessed when available. Histological review for benign and malignant features was performed on sections from a subset of 403 individual mice from 12 distinct CCMT groups representing either short (XAC2, ZIE2, HIP, ZIF2, LUF and LEM), medium (LOT, TOP, LUV, 266-Het) or long (NUK and LIL) survival groups, based on each respective groups’ median survival time (Figure 2C: “short” and “long” groups represent either the lowest 30% or highest 30% survival percentile respectively). Histological features representative of benign changes included mesothelial thickening, plaque development, paucicellularity within mesothelial layers, presence of giant cells and regular nuclei without atypia. Representative malignant histological changes included overt tumour (100% sarcomatoid), hypercellularity, nuclear atypia (multiple nucleoli, coarse chromatin), the presence of mitotic figures and invasion of surrounding tissue (Figure 4).

A range of benign/fibrotic and neoplastic histological features was observed in all asbestos-exposed mice. Frequency table analysis and logistic regression used to compare histologically features between short and long survival groups indicated that CCMT mice whose samples contained plaques (p = 0.004), regular nuclei (p = 0.019), mitotic figures (p < 0.001) and areas of invasion (p = 0.022; Table 1) were significantly associated with short survival. No significant differences between short and long survival groups were found for thickening, presence of giant cells, overt tumour, or regions of hypercellularity (Table 1).

Consistent with our previous observations of phenotypic traits, univariate analysis revealed no histological variables that were significantly associated with disease progression (Supplementary Table S3). Conversely, for disease latency, the presence of plaques (p = 0.033) and regular nuclei (p = 0.015) were associated with groups with shorter latency, while regions of paucicellularity in the mesothelial layer (p = 0.002), hypercellularity (p = 0.038) and invasion (p = 0.012) were associated with longer disease latency (Table 2). Multivariate analysis revealed that presence of plaques was independently associated with shorter disease latency time (adjusted p = 0.035), whereas the presence of paucicellular regions (adjusted p = 0.003) and evidence of invasion (adjusted p = 0.016) were independently associated with longer disease latency (Table 2).

Cox proportional hazards regression was used to analyse the overall (OS) and cancer-specific survival (CSS) data from individual asbestos exposed animals. Mice that survived until the end of the study (18 months) were censored for both OS and CSS.

Features found to be good prognostic factors for OS included paucicellular regions (HR 0.66; 95%CI 0.52–0.84; p < 0.001) and regions of hypercellularity (HR 0.31; 95%CI 0.14–0.70; p = 0.005). Regions of regular nuclei (HR 2.19; 95%CI 1.17–4.11; p = 0.015) and mitotic figures (HR 1.51; 95%CI 1.13–2.02; p = 0.005) were poor prognostic factors (Table 3). In multivariate analysis only paucicellularity was independently associated with overall survival with a 44% reduction in the hazard of death from any cause (Table 3). Similarly for cancer specific survival, paucicellularity (HR 0.65; 95% CI 0.51–0.84; p = 0.001) and hypercellularity (HR 0.25; 95% CI 0.11–0.57; p = 0.001) were considered good prognostic factors, whereas regular nuclei (HR 2.47; 95% CI 1.22–4.99; p = 0.012) and mitotic figures (HR 1.87; 95% CI 1.35–2.59; p < 0.001) were significant poor prognostic factors in univariate analysis (Table 4). Paucicellularity (adj. HR 0.74; 95% CI 0.57–0.97; adj. p = 0.027) and the presence of mitotic figures (adj. HR 1.68; 95%CI 1.94–2.36; adj. p = 0.003) were independently associated with improved cancer specific survival in multivariate analysis (Table 4).

Taken together the histological data from a subset of asbestos exposed CCMT mice indicates that mice develop many of the ARD related histological features that are characteristic of human mesothelioma.

To identify candidate modifier genes associated with asbestos related traits, we used the GeneMiner bioinformatic portal. For each trait, multinomial logistic regression models were fitted at each locus and ANOVA chi-square tests used to estimate the p-value of association. A summary of candidate modifier genes associated with each phenotypic trait is provided in Figure 5; Table 5. With respect to overall survival/latency as traits, we identified major effect QTL (LOD ≥5.8, founder coefficient range −100 to +150) on chromosomes 6, 12, and X (Figure 5A), with numerous protein-coding genes located around each peak QTL (Table 5; Supplementary Table S4). Similarly, eight additional QTL located on chromosomes 2, 3, 4 (two QTL), 5, 9, 16 and X were identified when using ascites volume as the trait; although the founder coefficient range was considerably smaller (−1 to +1.5), suggesting each of these genes alone would only have a minimal influence on ascites volume (Table 5; Supplementary Table S5; Supplementary Figure S1).

Of the 97 CCMT candidate modifier genes identified across the three major effect QTL for ARD survival/latency, four occurred at the peak QTL on chromosome 6 (Kcnd2, Tspan12, Ing3, Cped1); six on chromosome 12 (Matn3, Wdr35, Ttc32, Osr1, Nt5c1b, Rdh14); and 17 on the X chromosome (Il13ra1, Pgrmc1, Septin6, Ndufa1, Nkrf, Steep1, Akap17b, Zcchc12, Lonrf3, Slc25a5, Ct47, Gm14819, Rhox1, Gm14569, Gm10486, Akap14, Dock11; Figure 5B, Table 5; Supplementary Table S4). While no individual candidate modifier gene was observed to have a significant influence over any other gene at each respective QTL, many are known to be associated with advanced or metastatic cancers (Supplementary Table S4), suggesting a complex polygenic interplay between host genes affects ARD survival/latency.

We next sought to assess the potential for CCMT candidate modifier genes to influence disease outcome in human mesothelioma. All CCMT candidate modifier genes spanning the major effect QTL were converted into homologous human gene symbols and univariate Cox regression analysis performed on human mesothelioma RNAseq datasets to identify genes associated with age at the time of surgery and age at the time of diagnosis in the Bueno (Bueno et al., 2016) and TCGA-MESO cohorts, respectively. Expression of two candidate genes, odd-skipped related 1 (OSR1, HR 1.4, p = 0.026) and cadherin-like and PC-esterase domain containing 1 (CPED1; HR 1.5, p = 0.016), were significantly associated with a poorer outcome in the Bueno cohort (Figure 6A; Supplementary Table S6). Likewise, expression of five of the six genes identified in the TCGA cohort (Figure 6B; Supplementary Table S7), namely, HCLS1 binding protein 3 (HS1BP3, HR 2, p = 0.016); interleukin 13 receptor alpha 1 (IL13RA1; HR 2.1, p = 0.0057); LSM8 homolog, U6 small nuclear RNA associated (LSM8, HR 1.9, p = 0.021); NADH: ubiquinone oxidoreductase subunit A1 (NDUFA1, HR 1.8, p = 0.35) and testin LIM domain protein (TES, HR 2.1, p = 0.0075) were significantly associated with poorer outcome. Only expression of tetraspanin 12 (TSPAN12, HR 0.55, p = 0.03) was associated with improved outcome in the TCGA-MESO dataset (Figure 6B; Supplementary Table S7). No human homologues of CCMT candidate modifier genes were common to both cohorts.