Comparative Transcriptome Analyses Highlight Distinct Pathogenetic Mechanisms for Pleuropulmonary Blastoma and Congenital Pulmonary Airway Malformations

Congenital cystic lung malformations are uncommon. They encompass several types of rare pulmonary disorders, including congenital pulmonary airway malformation (CPAM), one of the most frequent pediatric lung lesions with an incidence of about 1:3000 births, and pleuropulmonary blastoma (PPB), a very rare lung neoplasm often confounded with CPAM [1,2]. In addition to sharing cystic morphological features, PPB and CPAM are thought to originate from developmental defects that impair the lung morphogenesis program. However, the molecular mechanisms underlying the pathogenesis of these clinically significant conditions remain largely elusive.

PPB is the most common pediatric primary lung neoplasm [3]. It manifests in early childhood and progresses through three distinct clinicopathological subtypes. Type I PPB consists of low-grade epithelial cystic lesions that are described to contain rhabdomyoblasts and immature cartilage cells within the cyst walls [4,5]. Overgrowth of these primitive mesenchymal cells may result in a mixed cystic and solid neoplasm (Type II), which can further evolve into a purely solid high-grade sarcoma with poor prognosis (Type III) [4]. Early diagnosis of Type I PPB is associated with improved outcomes, whereas delayed detection or failure to resect cystic PPB worsens prognosis, as Type I lesions may progress into aggressive Type III tumors within a few years [6,7].

PPB is the hallmark tumor of the DICER1 syndrome (also known as DICER1-related tumor predisposition; MONDO: 0100216), an autosomal dominant condition predisposing individuals to multiple cancers [8]. The DICER1 gene encodes a ribonuclease (RNase) III essential for processing pre-microRNAs (miRNAs) into mature, functional miRNAs [9]. DICER1 syndrome typically involves a germline loss-of-function mutation in one allele of the DICER1 gene, coupled with an acquired somatic missense mutation in the second allele [10]. These tumor-specific mutations usually occur in five hotspots, codons E1705, D1709, G1809, D1810, and E1813, affecting critical amino acids within the metal-ion binding catalytic cleft of the RNase IIIb domain of the DICER1 protein [11,12]. The resulting neomorphic DICER1 function disturbs miRNA biogenesis, selectively impairing miRNA-5p processing and favoring miRNA-3p production from the double-stranded precursors [13,14]. These quantitative and qualitative changes in the miRNA expression panel ensuing from DICER1 mutations may thus perturb target gene expression and contribute to PPB pathogenesis.

Diagnosing PPB remains challenging due to its overlap of clinical and histopathological features with other cystic lung malformations. Type I PPB is particularly prone to misidentification as CPAM, which generally has a more favorable outcome [2]. Accurate recognition of PPB at the Type I stage is therefore critical, as it represents the most curable phase of the disease.

A complex set of classification systems and terminologies has emerged over the years to describe CPAM. Historically, CPAM was divided into five subtypes based on cyst location along the tracheobronchial tree, cyst size, and histological characteristics [1,15]. However, cystic Type I PPB and Type IV CPAM are now considered to be the same pathological entity [16]. The widespread use of prenatal ultrasound has also prompted to revisit CPAM and to group them into two broad categories: Type 1, characterized by larger cysts and frequent mucous cell clusters, and Type 2 featuring smaller cysts [17,18]. However, this classification remains debated, as several CPAM cases exhibit overlapping features of both types [2,18,19].

Despite morphologic similarities, notable distinctions between PPB and CPAM have been documented. PPB cases frequently present an alveolar-type epithelial lining, whereas CPAM cysts are predominantly bordered by a ciliated respiratory bronchial-type epithelium [1,15,17]. Moreover, CPAM lesions very rarely contain primitive mesenchymal cells or cartilage nodules within their septa [4]. Finally, proliferation in CPAM samples is comparable to that observed in controls, in contrast to the increased mesenchymal proliferation commonly seen in PPB lung tissues, which likely reflects neoplastic expansion and disease progression [20,21]. While PPB is associated with DICER1 pathogenic variants, the molecular etiology of CPAM remains unresolved. Nevertheless, mutations in genes such as KRAS, FGFR2, and TP53 have been identified in CPAM cases [19,22,23].

Comparative expression analysis of candidate genes in PPB and CPAM lung biopsies has revealed differences in gene expression in PPB that were not observed in CPAM, indicating distinct molecular signatures [21,24]. Among the genes differentially expressed in PPB, several encode transcription factors (TF) and signaling molecules important for lung development [21,25]. Furthermore, misexpression of genes associated with epithelial-mesenchymal transition (EMT) suggests that this process may contribute to PPB pathogenesis. Altered gene expression is expected in PPB due to biallelic compound mutations in the DICER1 gene and the resulting dysregulation of miRNA biogenesis. For instance, miR-125a-3p negatively regulates lung branching as well as the expression of the ubiquitous TF YY1. It is also abnormally overexpressed in PPB lung tissues and in patient sera, a possible explanation for the decreased expression of the YY1 protein in PPB biopsies [21,24,26]. However, despite the identification of candidate genes, the sequence of molecular events driving PPB development remains to be elucidated.

To gain a comprehensive view of the molecular impact occurring in PPB and CPAM, we performed a transcriptome analysis of lesions at the genome-wide level, leveraging valuable patient lung specimens for these two rare lung diseases. RNA-seq data revealed that each pathology exhibits a unique expression profile, indicating distinctive molecular mechanisms. Analysis of the miRNA repertoire in PPB lung samples confirmed the predominance of miRNA-3p over miRNA-5p. It also unveiled the potential role of specific dysregulated miRNAs in the misexpression of components of the Polycomb Repressive Complex 2 (PRC2), with repercussions on developmental TF genes and cancer.

This study was conducted using anonymized specimens from control, CPAM, and PPB patients from the Department of Anatomo-Pathology of Hôpital Necker-Enfants Malades (Paris, France) and from McGill University/Segal Cancer Centre (Montréal, QC, Canada). CCAM and PPB tissue samples were obtained from surgical resections of the respective malformations. Decisions regarding whether to perform surgery and its timing were made at the discretion of the participating centers. All procedures were conducted according to the standard protocols of these centers, each of which is a certified or well-established pediatric surgery facility. Tissue samples for the malformation site were preserved for research purposes, either as frozen or formalin-fixed paraffin-embedded (FFPE) specimens. Informed consent was obtained from all subjects, and all clinical investigations were conducted in accordance with the principles expressed in the Declaration of Helsinki. The procedures were approved by the respective institutional ethics committee and by the CHU de Québec-Université Laval research ethics committee (project MP-20-2020-4809). When possible, the sex of specimens was determined by RNA-seq or RT-qPCR analyses of the Y-specific lung-expressed gene, UTY [27]. Characteristics of patients are summarized in Table S1.

Alcian Blue staining, for the detection of acid mucus-producing cells, and immunostaining experiments were performed as described [28]. IHC staining was revealed with DAB⁺ Substrate Chromogen kit (3,3′-diaminobenzidine; Dako K3468, Agilent Technologies, Santa Clara, CA, USA). IF labelling was performed with Vectastain ABC kit (Vector Laboratories, Newark, NJ, USA) and HRP activity was detected with the Tyramide Signal Amplification (TSA) plus Cyanine 5 (1/150; NEL745001KT Akoya Biosciences, Marlborough, MA, USA). Slides were counterstained with hematoxylin or nuclear Fast Red. For IF studies, nuclei were visualized with DAPI staining. Primary and secondary antibodies are listed in Table S2. IHC sections were imaged with a Leica DMR microscope (Leica Camera, Wetzlar, Germany) and a HP-PL-Fluostar 20× 0.50, while IF sections were imaged with a Zeiss Axio Observer 7 microscope (Jena, Germany) inverted with a Zeiss Axiocam 705 color camera (Jena, Germany) and a Plan-Apochromat 20× 0.80NA or with a Zeiss microscope LSM 900 AiryScan II (Jena, Germany) inverted with a Plan-Apochromat 20× 0.80NA.

Quantification of SOX2 IF staining was performed using the QuPath open-source software for digital pathology image analysis (v0.6.0; https://qupath.github.io, accessed on 27 June 2025) [29]. Epithelial nuclei were identified by applying the Cell Detection command to the Blue Channel (DAPI) and using QuPath’s integrated segmentation algorithms to ensure accurate nuclear detection. For each detected nucleus, the total SOX2 signal intensity from the RED channel was calculated. To normalize the nuclear size, a ratio of SOX2 signal per nuclear area (mm²) was computed. The basal epithelial length was measured using the Freehand Polyline annotation tool within QuPath. For each image, values were obtained from 55 to 567 epithelial cells, and 3 to 7 images were analyzed per specimen.

Total RNA from four PPB (PPB1-4), four CPAM (CPAM1-4), and four control (CTL1-4) lung biopsies was extracted as previously described (Table S1) [21]. cDNA libraries were constructed using SMART-Seq v4 Ultra-low input RNA Kit (Takara Bio USA, San Jose, CA, USA), and paired-end sequenced at 100 bp reads on Illumina Novaseq 6000 (San Diego, CA, USA) at the Next-Generation Sequencing (NGS) Platform, Genomics Center, Centre de recherche du CHU de Québec-Université Laval (Québec, QC, Canada). The quality of the sequencing reads was assessed using FastQC software (v0.12.1; http://www.bioinformatics.babraham.ac.uk/projects/fastqc, accessed on 1 March 2023). Reads were trimmed with Trim Galore (v0.6.5; https://github.com/FelixKrueger/TrimGalore, accessed on 19 November 2019) and mapped with STAR (v2.7.9a) [30] to the human reference genome from ENSEMBL GRCh38.p14 (release 3 February 2022) [31]. Between 77% and 87% of the reads corresponded successfully to the human genome (Table S3A). The mapped sequencing data were then processed with the software featureCounts to obtain the number of reads associated with genes (v2.0.3) [32]. Among the aligned reads, an average of 44% mapped to exon sequences that were further considered for analysis (Table S3A). Gene expression levels were represented as transcripts per million (TPM). Genes with a mean of TPM < 1 in both control and pathological conditions were filtered out before performing DESeq2 package (Bioconductor) differential expression analysis with an alpha parameter of 0.05 (v1.42.0) [33]. Genes were considered differentially expressed (DE) when they exhibited at least a 2-fold difference (fold change; FC) in expression between control and pathological conditions, with a false discovery rate (FDR) adjusted p value (padj) < 0.05. The value of a |Log2FC| ≥ 1 was chosen considering the small number of replicates and the important biological variation between specimens. Lists of the differentially expressed genes (DEG) are provided in Table S4.

Gene Set Enrichment Analysis (GSEA) was performed on TPM data from CPAM samples versus control samples and from PPB samples versus control samples with the GSEA software (v4.2.3) using the MSigDB database [34]. FDR-adjusted Q value < 0.05 and |Normalized Enrichment Score (NES)| > 1 were used as thresholds to define statistic gene set enrichment.

Total RNA from three PPB (PPB1-3) and four control (CTL1-4) lung biopsies was used for miRNA-seq assay (Table S1). After a purification step with GenElute PCR clean-up kit (Sigma-Aldrich, St. Louis, MI, USA), an appropriate range of cDNA fragments ~150 bp was extracted using a Pippin Prep instrument (Sage Science, Beverly, MA, USA). miRNA libraries were constructed and paired-end sequenced for 100 bp reads on the Illumina Novaseq 6000 of the NGS platform mentioned above. Data quality was evaluated using FastQC software, and reads were trimmed with Trim Galore. Only R1 reads were mapped with the software Bowtie (v1.3.1) [35] to the human reference genome from ENSEMB GRCh38.p14. More than 99% of the reads mapped successfully to the human genome (Table S3B). The mapped sequencing data were processed with the software featureCounts to obtain the number of reads associated with genes (v2.0.3) [32]. Among the aligned reads, an average of 84% mapped to miRNA sequences that were further analyzed (Table S3B). miRNA expression levels were represented as reads per million (RPMs). miRNAs with a mean of RPM < 1 in both control and PPB conditions were filtered out before performing DESeq2 package (Bioconductor) differential expression analysis with an alpha parameter of 0.05 (v1.42.0) [33]. miRNAs were considered to be differentially enriched when they exhibited at least a 2-fold change in expression between PPB and control conditions with a false discovery rate (FDR) padj < 0.05. Differentially expressed miRNAs are listed in Table S5.

Genomic DNA from the four PPB lung specimens (PPB1-4) was prepared using the Agilent SureSelect^XT Low Input Human All Exon V7 (Agilent Technologies). 200 ng of DNA was fragmented on a Covaris instrument (Covaris, Woburn, MA, USA), and fragments of 150–200 bp were purified and used for library preparation. Following hybridization to the target-specific capture library and purification using streptavidin-coated beads (Dynabeads MyOne Streptavidin T1, Thermo Fisher Scientific, Waltham, MA, USA), libraries were amplified for 10 cycles with indexing primers. Exome libraries were paired-end sequenced for 150 bp reads on the Illumina Novaseq 6000 of the NGS platform previously mentioned. Data quality was evaluated using FastQC software, and reads were trimmed with Trim Galore before being mapped with the software Bowtie (v2.5.2) [35] to the human reference genome from ENSEMB GRCh38.p14. More than 99% of the reads mapped successfully to the human genome, and an average of 78% aligned to exon sequences (Table S3C). Variant detection was done using the variant calling software bcftools (v1.19) [36]. We found 17,202 variants, Single Nucleotide Polymorphism (SNP) and Insertion-Deletion (INDEL), for PPB1, 15,265 for PPB2, 16,996 for PPB3, and 17,469 for PPB4. Using this variant detection approach, we were able to identify DICER1 mutations in lung biopsy samples from PPB patients (Figure 1).

Genomic DNA from CPAM specimens was recovered from paraffin sections using QIAamp DNA FFPE Advanced Kit (Qiagen, Venlo, The Nethernands). DNA fragments encompassing the KRAS G12 protein position were amplified by PCR with Taq DNA Polymerase (Bio Basic, Markham, ON, Canada) with three different primer sets: #1: Forward: 5′-TGTTGGATCATATTCGTCCACAA-3′, Reverse: 5′-GTACTGGTGGAGTATTTGATAGTGT-3′; #2: Forward: 5′-CGATACACGTCTGCAGTCAAC-3′, Reverse: 5′-TCATGAAAATGGTCAGAGAAACC-3′; and #3: Forward: 5′-GTGTGACATGTTCTAATATAGTCAC-3′, Reverse: 5′-AAAACAAGATTTACCTCTATTGTTGG-3′. DNA fragments were sequenced by Sanger technology, and the KRAS G12 status was established by analysis of chromatograms.

To estimate the chromatin state of the 92 genes coding for TF and presenting a gain in expression in PPB samples, we used ChIP-seq peak called data for human lung from 16-, 30- and 40-years old individuals obtained from ENCODE for the three following histone post-translational modifications: H3K27ac (ENCODE:ENCSR510RPC/GEO:GSE175318; ENCODE:ENCSR540ADS/GEO:GSM906395; ENCODE:ENCSR494WCX/GEO:GSE167619), H3K27me3 (ENCODE: ENCSR800MYV/GEO:GSE175189; ENCODE:ENCSR204NFO/GEO:GSM1220283; ENCODE:ENCSR975GDL/GEO: GSE167807), and H3K36me3 (ENCODE:ENCSR793LDG/GEO:GSE175165; ENCODE:ENCSR671NXL/GEO:GSM956014; ENCODE:ENCSR906YQO/GEO:GSE167769). Only peaks present in the three data sets with an overlap of at least 50% were kept using the software BEDTools option intersect (v2.31.0) [37]. The genome-wide distribution of the three histone marks was analyzed. The combination of histone marks at the gene level enabled to define a 3-state model. The sole presence of H3K27me3 corresponds to heterochromatin, while H3K27ac and H3K36me3 combined or alone are associated with permissive chromatin. Finally, the absence of all three marks, or the presence of the repressive H3K27me3 mark with the activating H3K36me3 and H3K27ac marks, was considered an undetermined state.

Statistical analyses were performed in R (v4.3.1). Random sampling was used to study the variations in proportions of chromatin states in seven different groups of genes: (i) ‘all genes’, (ii) ‘all TF’, (iii) ‘DEG down in PPB’, (iv) ‘DE TF down in PPB’, (v) ‘DEG up in PPB’, (vi) ‘DE TF up in PPB’, and (vii) ‘DE TF with a gain in expression in PPB’. Gain in expression was defined as expression < 1 TPM in controls, but ≥1 TPM in PPB specimens. To determine if the gain in expression of the 92 DE TF genes in PPB was random or significant, a proportion test based on the chromatin states was performed. To do so, 92 genes from the ‘all genes’, ‘all TF genes’, and ‘DEG up in PPB’ categories were randomly sub-sampled from the genome 10,000 times [38]. Observed and simulated Chi-square statistics were then compared to assess the observed differences in chromatin state proportions of the 92 DE TF genes with a gain in expression in PPB. The difference was statistically significant with a p value < 2.2 × 10⁻¹⁶ for all three comparisons, indicating that the variation in proportions was not due to chance. Significant differences in gene expression (Log2(mean of TPM + 1)) between PPB and control specimens for the seven different groups of genes mentioned above were assessed using the pairwise Wilcoxon rank sum test.

To characterize transcriptome-wide changes occurring in PPB and CPAM, we performed a comparative RNA-seq analysis of lung lesions from four PPB patients and four CPAM patients, as well as of normal lung tissues from four controls (Table S1). We first examined the DICER1 sequences to determine whether mutations were present in PPB patients. All four PPB specimens carried a missense mutation in the RNase IIIb catalytic domain, corresponding to one of the previously reported mutational hotspots (Figure 1) [12,14]. For patients PPB1, PPB3, and PPB4, the variant frequency ranged from 53% to 61% indicating a heterozygous mutation. In contrast, PPB2 exhibited an allelic frequency of 99%, suggesting that the mutation was either homozygous or hemizygous with a loss of heterozygosity. The availability of DNA from the four PPB lung tissues also enabled whole exome sequencing, which confirmed the presence of these missense mutations. Again, the variant frequency of the mutation was higher in PPB2 (89%) than in the other three specimens (24% to 42%). Moreover, whole exome sequencing data from PPB1, PPB3, and PPB4 revealed a loss-of-function mutation caused by nucleotide deletion, resulting in a frameshift and premature termination of the DICER1 protein. No loss-of-function mutation was found in the PPB2 specimen. Since we did not have access to DNA from paired normal samples, we could not determine whether the DICER1 mutations identified in the lung lesions were germline predisposition variants or tumor-specific alterations.

RNA-seq data did not reveal DICER1 mutations in CPAM specimens or in controls, with the exception of CPAM3 lung tissue, which carried the DICER1 missense mutation Y1835C. This mutation lies in a region with no known functional annotation, located between the RNase IIIb and the double-stranded RNA binding domain (Figure 1). Previous studies have shown that carriers of non-hotspot DICER1 mutations do not present biases in miRNA biogenesis, suggesting that the functional impact of the Y1835C mutation on miRNA processing is probably limited [39].

Further analysis of the whole exome sequence for cancer gene mutations using the Catalogue of Somatic Mutations in Cancer (COSMIC) and the IntOGen databases revealed that the PPB4 lesion carried a missense mutation in the TP53 gene (Figure 1) [40,41]. This C275F mutation, located in the DNA-binding domain of the TP53 protein, is considered a hotspot in several carcinomas, including lung adenocarcinoma and lung squamous cell carcinoma (cBioPortal.org, accessed on April 2018). The high allelic frequency (95%) of this TP53 mutation suggested it was either homozygous or hemizygous with a loss of heterozygosity. Biallelic TP53 mutations are a frequent event in PPB [14].

To determine whether PPB and CPAM have similar transcriptomic impacts in the lung, RNA-seq assays were carried out on PPB, CPAM, and control specimens. Principal Component Analysis (PCA) was conducted on normalized gene counts to assess the variability within the RNA-seq dataset. This analysis revealed a clear clustering of specimens according to their pathological status (Figure 2A). Despite the limited sample size, PC1 captured a well-defined separation of PPB from CPAM and controls, while PC2 showed some overlap among CPAM1, CPAM2, and control specimens. Clustering was further supported by the comparison of the differentially expressed genes (DEG) exhibiting at least a 2-fold difference in expression with a padj < 0.05. A Venn diagram showed 2416 DEG between PPB and control, 2301 DEG between PPB and CPAM, and only 250 DEG between CPAM and control specimens. Notably, the DEG profiles of PPB versus control and PPB versus CPAM superimposed by more than 70% (1725 DEG; Figure 2B). This was reflected in the comparable numbers of downregulated and upregulated DEG in PPB: 890 and 876 DEG were downregulated in PPB versus control and PPB versus CPAM comparisons, respectively, while 1526 and 1425 DEG were upregulated (Figure 2C,D). In contrast, only seven DEG were downregulated in CPAM versus control, while 243 DEG were upregulated, a pattern similar to that reported by Zhang et al., 2024 (Figure 2E) [42].

Heatmap analysis of DEG expression across all 12 samples confirmed the hierarchical clustering by pathology. For the 2416 PPB DEG, PPB specimens were grouped together while CPAM and control samples clustered jointly (Figure 2F). Similarly, for the 250 CPAM DEG, CPAM specimens were sorted together (Figure 2G). However, as shown by the PCA analysis, heterogeneity between CPAM cases was observed. CPAM1 and CPAM2 formed a subgroup that clustered with the control samples, which may explain the moderate number of DEG between CPAM and controls. In contrast, CPAM3 and CPAM4 were part of another subgroup with markedly increased DEG expression (Figure 2G). This subgroup also included PPB3, which had been identified as a Type I cystic form. This may justify why its transcriptome shares molecular characteristics with CPAM, while retaining a distinct PPB profile (Figure 2F,G; Table S1). The dual behavior of PPB3 may offer molecular clues on the intermediate phenotypes that arise during PPB tumor progression. In summary, transcriptome-wide data indicate that CPAM and PPB are globally different pathological entities with a unique molecular signature, supporting our previous findings [21]. Furthermore, transcriptional heterogeneity among CPAM specimens highlights the complexity of this disease.

Gene Set Enrichment Analysis (GSEA) was applied to CPAM and PPB datasets. For CPAM, upregulated transcripts were significantly enriched for genes specific to basal, ciliated, secretory, goblet, and club cell types from the bronchial epithelium, known to line CPAM cysts (Figure 3A–C) [1]. GSEA also revealed several terms related to cilia, such as cilia development and movement, axoneme assembly, microtubule motor activity, and dynein complex (Figure 3A). This analysis agreed with previous studies that have suggested the importance of cilia in CPAM disease [42,43,44]. Heatmap representation of upper airway cell type-specific gene expression data confirmed the higher expression of proximal cell markers for basal, ciliated, goblet, club, and submucosal gland (SMG) cells (Figure 3D). However, the four CPAM specimens behaved differently. Increased expression of bronchial cell markers was more important in CPAM3 and CPAM4 than in CPAM1 and CPAM2, which presented similarities with controls (Figure 3D). Elevated expression of proximal cell markers was also detected in PPB3, consistent with the clustering of PPB3 with CPAM3 and CPAM4. Moreover, the high expression of specific bronchial cell genes in CPAM3, CPAM4, and PPB3 correlated with the increased expression of the TF sex-determining region Y (SRY)-box 2 (SOX2) gene, a lineage commitment marker of proximal endodermal lung progenitor cells (Figure 3D) [45]. Increased SOX2 expression in CPAM was previously reported [20,43,46]. In addition, gain of Sox2 expression in mice hampers lung branching and causes cystic formation, supportive of a role in CPAM pathogenesis [46,47]. Interestingly, expression of FOXA1 and SIX1 genes, which encode the winged-helix TF FOXA1 and the sine oculis homeobox transcriptional regulator SIX1, respectively, was also augmented in CPAM3, CPAM4, and PPB3. FOXA1 and SIX1 proteins were both shown to act positively on SOX2 expression by direct binding to SOX2 specific enhancers [48,49]. Moreover, the loss of function of Six1 or Foxa1 in mice affects lung morphogenesis [50,51]. Finally, SOX2 was shown to directly activate the expression of the basal cell specific Trp63 gene [46]. Together, these correlations suggest that the concomitant increased expression of FOXA1 and SIX1 can trigger SOX2 upregulation and consequently perturb bronchial epithelial cell differentiation in CPAM lung tissues.

To validate the changes in expression of proximal cell genes revealed by RNA-seq, we performed immunostaining assays. For these analyses, we considered the KRAS status of the CPAM specimens. Somatic gain-of-function mutations in the KRAS (Kirsten rat sarcoma viral oncogene homolog) gene are established driver mutations for the development of several cancers, including lung mucinous adenocarcinoma, a condition that can arise from CPAM [52,53,54]. KRAS mutations at codon 12 were also reported in CPAM cases containing mucinous cell clusters, raising questions on the cancer risk associated with the cystic disease [18,19,22,55,56,57]. Using Sanger sequencing of PCR-amplified KRAS specific DNA fragments, we determined the KRAS status of CPAM lung biopsies. Out of 25 CPAM lesions, 12 (48%) carried either a p.G12D (83%) or a p.G12V (17%) constitutive active mutation (KRAS^mut), while no KRAS mutation (KRAS^wt) was detected in 13 specimens (52%) (Table S1). The prevalence of KRAS p.G12D in CPAM cases was previously reported [57]. Foci of goblet cell hyperplasia and regions of pseudostratified epithelium and papillary architecture were frequent in KRAS^mut specimens but rarely seen in KRAS^wt ones (Figure 4) [22].

Variations between specimens were detected, but an overall trend of augmented SOX2 protein expression in KRAS^mut specimens was observed, as shown by IF staining and SOX2 intensity measurement using QuPath (Figure 4A–C and Figure S1). SOX2 staining in CPAM KRAS^wt tissues was comparable to that of controls (Figure 4A,B). As well, expression of SIX1, FOXA1, and TP63 was also increased in KRAS^mut specimens versus CPAM KRAS^wt and controls (Figure 4D–I,M–O). A similar result was obtained for FOXJ1, a marker of ciliated cells and CC10, which labels club secretory cells (Figure 4P–U). Mucus-producing cells detected by Alcian Blue staining were also more frequent in CPAM KRAS^mut specimens (Figure 4V–X). Ki-67 staining was similar in the three conditions analyzed, indicating that the increased protein staining detected in CPAM KRAS^mut tissues was not due to enhanced proliferation (Figure 4J–L).

Thus, gene expression changes revealed by RNA-seq assay were confirmed by differential protein expression detection in a larger CPAM cohort. It also showed a correlation between the presence of a KRAS constitutive active mutation, altered epithelial architecture, and variations in expression of bronchial cell specific markers.

While CPAM DEG were almost exclusively upregulated, PPB DEG included both downregulated and upregulated genes (Figure 2). GSEA for PPB downregulated transcripts showed significant enrichment for genes associated with cell types from the alveolar epithelium, such as alveolar epithelial Type 1 cells that permit gas exchanges with the underlying vascular endothelium, and Type 2 cells that secrete surfactant proteins (Figure 5A,B and Figure S2) [58]. Moreover, a significant enrichment for downregulated genes specific to capillary endothelial cells, alveolar fibroblasts, aerocytes, pericytes, and the Gene Ontology (GO) term “lung alveolus development” was observed, suggesting a loss of cell identity in the distal lung. Heatmap analysis clearly demonstrated the decreased expression of markers for alveolar epithelial and mesenchymal cells (Figure 5C,D). No significant variation was detected in SOX9 and ID2 expression, both of which encode lineage markers of distal epithelial lung progenitors, suggesting that cellular commitment was not affected (Figure 5C) [58]. In contrast, expression of GATA6, NKX2-1, and ICAM1, known to be important for distal airway cell differentiation, was diminished in PPB [59,60,61]. Thus, impaired alveolar cell differentiation likely contributes to PPB pathogenesis.

GSEA also revealed decreased expression of genes associated with epithelial cell characteristics, including cell adhesion, cell junctions, epithelial cell development, and tube branching. Interestingly, epithelial genes known to be suppressed in EMT were found to be downregulated in PPB (Figure 5A,E). In parallel, upregulated transcripts were enriched in genes involved in mesenchyme-derived cell types and in mesenchyme development, as well as in genes highly expressed in EMT, supporting the notion that EMT contributes to the progression of PPB (Figure 5A,E) [21]. Finally, upregulated transcripts were enriched in genes involved in pediatric cancer markers, cell cycle, translation, ribosome assembly, and chromatin remodeling, suggestive of an enhanced proliferative status (Figure S2).

GSEA analysis of PPB data also unveiled enrichment in upregulated genes associated with embryonic organ morphogenesis, proximal-distal pattern formation, and limb development, categories that were not observed in CPAM upregulated DEG. Among these PPB upregulated genes, a significant number corresponded to TF included in the GO term “DNA-binding transcription factor binding activity” (Figure 5A).

We questioned what proportion of TF was abnormally expressed in PPB. Among the PPB DEG, 7.5% corresponded to upregulated TF and 2% to downregulated TF (Figure 6A). Interestingly, in the top 10, 50, and 100 upregulated DEG, TF represented respectively 80%, 56% and 39%. Thus, despite TF constituting only a small proportion of upregulated DEG, they were among the most strongly upregulated ones. This pattern was not observed for downregulated DEG, where TF corresponded to only 2% of the top 50 DEG (Figure 6A,B). Most upregulated TF genes (51%) were associated with a gain of expression, i.e., their expression was <1 TPM in controls, but ≥1 TPM in PPB specimens (Figure 6C). Thus, ectopic expression of developmental TF not normally found in the lung occurs in PPB.

As listed in Figure 6B, most of the upregulated TF were members of gene families having an established role in embryo development. One of the main representatives was the HOX gene family, which is involved in embryo patterning, axis specification, and organogenesis [62]. In PPB, 22 out of the 39 HOX genes presented changes in their expression, 15 showing a gain in expression (HOXA1, HOXA9, HOXA10, HOXA11, HOXA13, HOXC5, HOXC8, HOXC9, HOXC10, HOXC11, HOXC12, HOXC13, HOXD10, HOXD11 and HOXD13). These 15 HOX genes are not normally expressed in the lung [63]. Repression of developmental TF is known to be associated with heterochromatin, and silencing is mediated by the histone H3K27me3, an epigenetic negative mark (Figure 6D) [64]. Therefore, ectopic gain in TF expression could be explained by defective chromatin regulation due to the loss of H3K27me3. This hypothesis was further supported by the GSEA data showing an enrichment in upregulated genes for the “Polycomb Repressive Complex 2 (PRC2) targets” category (Figure 5A). PRC2 is the chromatin-modifying multi-protein complex that adds methyl groups to lysine 27 on histone 3 causing gene silencing [65].

To further determine if PPB TF DEG are mainly associated with a repressive chromatin state compared to other genes in the lung, we looked at the histone mark distribution over all genes of the entire human genome. To do so, we exploited ChIP-seq data for lungs of healthy human adults obtained from ENCODE for three histone post-translational modifications: H3K27ac, H3K36me3, and H3K27me3. H3K27ac is enriched in promoters and enhancers of transcribed genes, H3K36me3 is present in coding sequences of actively transcribed genes, while, as mentioned, H3K27me3 is a repressive mark [64,66,67]. Combinations of these three marks define three broad chromatin states: permissive, repressive, and undetermined (Figure 6E). A similar distribution of chromatin states was observed for the ‘all genes’, ‘all TF’, ‘DEG down in PPB’, ‘DE TF down in PPB’ and ‘DEG up in PPB’ groups, which were mainly associated with a permissive chromatin state (from 53% to 71%) with a low % corresponding to the repressive state (from 5% to 21%; Figure 6E). In contrast, ‘DE TF up in PPB’ and ‘DE TF with a gain in expression in PPB’ were found predominantly in a repressive chromatin state (50% and 78.5%, respectively) with a small % in the permissive state (28% and 4%; Figure 6E). To confirm that the pattern of chromatin states for the 92 DE TF showing a gain in expression in PPB was not due to arbitrary events, we randomly sub-sampled 92 genes from the ‘all genes’, ‘all TF genes’, and ‘DEG up in PPB’ categories from the genome 10,000 times as described in the Materials & Methods section. This analysis demonstrated that the enrichment in H3K27me3 observed for the ‘DE TF with a gain in expression in PPB’ group was significant, indicating that the chromatin state distribution was not due to chance (Figure 6E).

We also looked at the changes in transcript levels between PPB and control specimens for the seven different groups. No significant variation was detected for the ‘all genes’ and ‘all TF’ categories, whereas a significant difference was measured for the other five categories (Wilcoxon test, p value < 0.001; Figure 6F). Moreover, the expression levels observed in control specimens corroborated with the chromatin state deduced from the histone mark distribution. For the ‘DE TF with a gain in expression in PPB’ group, transcript levels were close to 0 in controls, supporting the repressive chromatin state, while significant expression levels were detected in PPB specimens (Figure 6F). Together, these analyses highlighted a gain of expression for a high proportion of TF associated with embryonic development but normally repressed in the developing lung. Correlation with the presence of the H3K27me3 repressive histone post-translational modification on these TF genes in normal lung tissue suggests a connection between their upregulation in PPB and impaired PRC2 function. This is further supported by the loss of H3K27me3 signal reported in PPB specimens [68].

DICER1 mutations in PPB favor the production of miRNA-3p over miRNA-5p [13,14]. miRNA-seq assay was therefore carried out on PPB and control specimens to identify differentially expressed miRNA-5p and -3p. Data obtained agreed with previous reports, as 98% of upregulated miRNAs were miRNA-3p, while miRNA-5p constituted 83% of the downregulated DE miRNAs in PPB specimens (Figure 7A,B). The majority (52%) of the remaining 17% downregulated DE miRNAs were miRNA-3p, for which miRNA-5p expression was also decreased, suggesting impaired synthesis of the pre-miRNA. Heatmap analysis revealed heterogeneity between the three PPB specimens studied. PPB1 and PPB3 presented similar miRNA expression profiles, while PPB2 showed a more pronounced deregulation of miRNA-3p expression. The latter might be explained by the sole presence of the missense mutation (Figure 1 and Figure 7A).

Despite differences in their miRNA patterns, the impact on gene expression was largely comparable between PPB specimens (Figure 2F). We thus hypothesized that miRNAs showing similar expression changes in the three PPB specimens should be involved in the disease. We separately examined upregulated miRNA-3p and downregulated miRNA-5p (Figure 7 and Figure S3). To identify functional categories statistically over-represented among the misregulated miRNA, a microRNA Enrichment Analysis and Annotation (miEAA; v2.1) with the Over-Representation Analysis (ORA) option was performed for each population [69]. For the upregulated miRNA-3p, miEAA revealed significant enrichment in lung development, lung disease, and lung cancer, including EMT, TP53, and cell cycle. Categories associated with gene regulation, transcription factors, and chromatin were also significantly enriched (Figure 7D). For the downregulated miRNA-5p, significant enrichment was obtained for categories associated with embryonic development, FGF pathway, cancer, miRNA processing, chromatin and translation (Figure S3B).

Attention was given to upregulated miRNA-3p. Among the DE upregulated miRNA-3p common to all specimens, some were known to target protein subunits of the PRC2 complex. Furthermore, loss of expression of SUZ12 or EED protein subunits was shown to cause an increased expression of transcription factors associated with embryonic development [70,71]. We thus decided to look closely at components of the PRC2 complex (Figure 7B,C). Interestingly, 68% of the DE TF showing a gain in expression in PPB, including the 15 HOX genes listed above, were found upregulated in non-small-cell lung murine cancer cells carrying an Eed deletion [72]. In PPB, we observed that miR199a-3p, miR-199b-3p, and miR-148b-3p, as well as miR-128-3p and miR148a-3p, all known to target SUZ12, were either upregulated or presented a trend for overexpression [73,74,75]. Moreover, miR-323a-3p, which was shown to target Eed in mice directly, was upregulated in PPB specimens [76]. To determine if miR-323a-3p overexpression correlates with changes in EED protein expression, we performed immunofluorescence staining on sections from control lung tissues and Type I, Type II, and Type III PPB lung lesions. In controls, EED protein was detected throughout the tissue as previously reported (https://www.proteinatlas.org/ENSG00000074266-EED/tissue/bronchus#, release 19 June 2023) [77]. In contrast, the EED signal appeared mosaic with regions devoid of EED expression in the PPB specimens tested (Figure 7E). Thus, the increased expression of specific miRNA-3p due to the DICER1 mutations in PPB may negatively affect PRC2 activity with repercussions on chromatin epigenetic landscape and expression of developmentally involved transcription factors not normally expressed in the lung.

The morphologic overlap between PPB and CPAM poses a significant challenge to accurate diagnosis. Improving the ability to distinguish PPB from CPAM would represent a significant advance in diagnostic precision and clinical management. To address this, we investigated the molecular differences between the two diseases by profiling their transcriptional landscapes at a genome-wide scale using comprehensive bulk RNA-seq and miRNA-seq analyses. The transcriptome-wide data support our previous gene expression analysis, indicating that PPB and CPAM are different pathologies possessing distinct molecular signatures [21]. Expression of proximal lung epithelium markers was specifically increased in CPAM (Figure 3). In contrast, distal lung epithelium and mesenchyme markers were less expressed in PPB, with the exceptions of mesenchyme markers associated with EMT and selected developmental transcription factors that showed increased and even gain of expression (Figure 5 and Figure 6). Apart from specimen PPB3, diagnosed as a Type I cystic form and presenting an RNA-seq profile sharing both PPB and CPAM characteristics, a clear difference was thus seen between PPB and CPAM specimens (Figure 2F,G).

We questioned if DEG common to PPB and CPAM transcriptomes, as compared to CTL, may participate in cyst formation for both pathologies. Among the 25 common DEG, EYA1 appears as a candidate that could play such a role (Table S4). EYA1 is a transcriptional co-activator of SIX1 [78]. In the lung, EYA1 and SIX1 were shown to interact synergistically in the control of epithelial branching and mesenchymal development by inhibiting SHH signaling activity, coherent with a role in cyst formation [79]. EYA1 and SIX1 expression levels were increased in PPB and CPAM specimens, with a more pronounced effect in PPB. This may contribute to the reduced SHH expression observed in PPB (Figure 3D; Table S4). However, in CPAM, there is no evidence for altered expression of SHH pathway members, suggesting that either the modest upregulation of EYA1 and SIX1 expression has minimal impact on SHH signaling, or that EYA1/SIX1 acts through a distinct cascade that remains to be identified (Table S4) [21,80].

A point to consider is that, given the limited number of CPAM samples analyzed, we disregarded their CPAM classification in the RNA-seq data analysis. This could explain the heterogeneity between transcriptomic data and the low number of DEGs retrieved compared to other studies with larger cohorts or that considered CPAM classification based on cyst size [42,43,44]. This raises the possibility that other genes not detected by our transcriptomic analysis may be shared between PPB and CPAM during cyst formation.

Analysis of CPAM RNA-seq data allowed us to distinguish two subgroups with distinct molecular signatures. CPAM3 and CPAM4 specimens showed strong expression of proximal pulmonary epithelial cell markers, a feature associated with Type 1 CPAM. In contrast, gene expression changes were less pronounced in CPAM1 and CPAM2 (Figure 2G) [17]. Type 1 CPAM is generally characterized by large cyst lesions lined with a pseudostratified epithelium and papillary structures, mucus cell hyperplasia, and the presence of a KRAS oncogenic mutation at codon 12 [17,18,57]. Although CPAM3 and CPAM4 displayed transcriptomic profiles consistent with Type 1 CPAM, they were clinically diagnosed as Type 2, and they lacked a KRAS mutation. Conversely, CPAM2, which had a transcriptomic profile resembling that of controls, was diagnosed as a Type 1 and carried a KRAS^G12V mutant allele (Table S1). These inconsistencies highlight the challenges of CPAM classification. Discrepancies between KRAS mutation status and clinical diagnosis may also reflect intralesional heterogeneity of the mutation, as previously reported [57]. Moreover, limited dissection of cystic tissue may hinder KRAS mutation detection by the Sanger sequencing method. Finally, the possibility that Type 1 CPAM cases lack KRAS mutations cannot be excluded.

For the validation by immunostaining of DEG encoding proximal lung epithelium markers, we considered the KRAS status of CPAM samples from a larger cohort. A clear correlation was established between the KRAS^mut status of CPAM lesions, increased expression of bronchial epithelial cell markers, presence of pseudostratified epithelium with papillary projections and abrupt transitions, and the occurrence of mucus islands. This agrees with previously established links between KRAS activating mutations and increased differentiation of secretory and ciliated cells, epithelial cell hyperplasia and mucinous cell clusters [57,81,82]. However, we did not observe that the thickening of the cyst epithelium in CPAM KRAS^mut specimens was associated with changes in cell proliferation, as observed in human airway basal cells overexpressing a constitutively active KRAS protein (Figure 4) [82]. Our findings concurred with the lack of change in lung epithelial proliferation seen in mice carrying a Kras^G12D mutation and presenting a defective lung branching phenotype that causes the formation of large, fluid-filled sacs reminiscent of lung cysts [83]. While these data further support the developmental origin of CPAM, they also suggest that the presence of KRAS oncogenic mutation may cause an intrinsic epithelial anomaly that leads to a disordered epithelial growth and differentiation of proximal airway structures and participates in cyst development.

The increased expression of SOX2 observed in CPAM KRAS^mut specimens may constitute an additional and accurate molecular marker to classify CPAM cystic lesions. Evidence from genetically modified mice has established a role for SOX2 in cyst formation [46]. Moreover, transcription of SOX2 is under the control of SIX1 and FOXA1, two TF with augmented expression in CPAM KRAS^mut specimens. It is therefore possible that the KRAS activating mutation may affect the expression of several molecules, including SOX2, SIX1, and FOXA1, and consequently impact bronchial epithelial cell differentiation in CPAM lesions.

On the other hand, SOX2, SIX1, and FOXA1 were also highly expressed in PPB3, a specimen positive for both DICER1 loss-of-function and missense mutations but negative for the KRAS mutation. This observation raises the possibility that misregulation of the SOX2, SIX1, and FOXA1 transcription factors occurs independently of the KRAS mutation status.

In PPB, genomic heterogeneity was observed between specimens. Notably, PPB2 only exhibited a DICER1 missense mutation, rather than the typical combination of a loss-of-DICER1 function mutation on one allele and a missense hotspot mutation on the other. Due to limited material, we were unable to determine whether the mutation present in PPB2 was homozygous or hemizygous with a loss of heterozygosity. Cases of tumor-specific, biallelic DICER1 pathogenic variants were reported to give rise to sporadic, non-syndromic PPB [8,10]. Furthermore, we cannot exclude the possibility that the DICER1 missense mutation in PPB2 is heterozygous and exerts a dominant-negative effect over the wild-type allele, as recently reported [84].

Despite its specific genomic status, PPB2 presented an RNA transcriptome profile similar to that of the three other specimens analyzed, each with a significant decrease in alveolar epithelial cell markers that distinguishes this pathology from CPAM (Figure 5) [43]. The reduced expression of alveolar epithelial cell markers did not result from fewer progenitor cells, since the alveolar epithelial progenitor cell markers SOX9 and ID2 did not present significant changes in their expression levels. Moreover, diminished expression of markers of specialized alveolar mesenchymal cells, such as capillaries, alveolar fibroblasts, aerocytes, and pericytes, suggested a widespread de-differentiation process occurring in the distal lung that might be indicative of tumoral progression. Tumor development was also supported by changes in the expression of genes promoting EMT as observed in the PPB specimens except for PPB3, coherent with the Type I classification of this sample [85].

One striking feature observed in all four PPB specimens was the abnormal upregulation of genes encoding developmental transcription factors (Figure 6). These TF are not normally involved in lung development and are primarily associated with heterochromatin-enriched regions, where their expression is repressed in healthy lung tissue. Genomic regions encompassing these genes are characterized by the presence of the histone H3K27me3 mark, product of the PRC2 complex catalysis, and essential for the transcriptional repression of PRC2 target genes [64]. In ES cells, it was shown that the PRC2 complex occupies numerous genes, most coding for transcription factors playing key roles in development, such as the HOX genes that are canonical PRC targets [65].

PRC2 complex comprises numerous protein subunits, including EED (for Embryonic Ectoderm Development). EED deficiency was shown to cause de-repression of TF genes due to the loss of H3K27me3 repressive mark [70,86]. Furthermore, inactivation of the PRC2 complex through EED or SUZ12 protein loss has been identified in human cancers, leading to abnormal expression of transcription factors linked to embryonic development and morphogenesis, including HOX TF [71,72]. Defective deposition of the repressive chromatin mark H3K27me3 has previously been reported in PPB, consistent with the aberrant gain of expression of developmental TF and the reduced expression of EED detected by immunostaining in PPB specimens (Figure 6 and Figure 7) [68]. The decreased EED expression in PPB may be explained by the increased expression of miRNA-323-3p, which is known to target EED. This suggests a mechanistic link between altered miRNA production, resulting from DICER1 mutations, and the abnormal expression of developmental transcription factors [76].

A close relationship exists between development and cancer. It is well established that misregulation of developmental TF genes can perturb downstream effector expression, leading to inappropriate activation of embryonic developmental cascades. This, in turn, can interfere with normal growth and differentiation programs, contributing to neoplastic transformation [87]. Aberrant activation of HOX genes has been reported in numerous cancers [88]. For instance, HOXC10 expression is abnormally upregulated in non-small cell lung cancer due to functional suppression of the PRC2 complex [89].

In conclusion, for the first time, and despite the limited sample size, integration of multi-omics analyses has enabled the identification of a potential pathogenetic mechanism in PPB that connects DICER1 mutations and altered miRNA profiles to impaired PRC2 activity, misexpression of developmental transcription factor genes, and tumorigenesis. To validate the significance of this pathway in PPB pathogenesis, additional RNA- and miRNA-seq experiments on larger cohorts will be required, along with functional assays to assess the impact of DICER1 mutations on other PRC2 complex components and to identify the transcription factors driving the oncogenic process. scRNA-seq assays should also be performed to dissect the dynamic cellular states involved. Similarly, further studies involving more patients are needed to clarify the relationship between KRAS status and SOX2 expression in CPAM pathology, and to disentangle the underlying molecular mechanisms.

The following supporting information can be found at: https://www.sciepublish.com/article/pii/822, Table S1: General characteristics of patients; Table S2: List of primary and secondary antibodies; Table S3: RNA-seq, miRNA-seq and whole exome statistics; Table S4: RNA-seq data—Gene expression PPB-CPAM-CTL; Table S5: miRNA-seq data—miRNA expression PPB-CTL; Figure S1: Box-and-whisker plot representation of SOX2 fluorescence signal along lung epithelium in control, CPAM KRAS^wt and CPAM KRAS^mut specimens; Figure S2: Statistically significant enrichment plots of terms associated to GSEA analysis of PPB transcriptome with an FDR-adjusted Q < 0.05 and a NES > 1; Figure S3: (A) Heatmap and hierarchical clustering of downregulated miRNA-5p common to the three PPB specimens analyzed. (B) A functional enrichment analysis of downregulated miRNA-5p was performed using microRNA Enrichment Analysis and Annotation (miEAA) with the Over-Representation Analysis (ORA) option surveying across GO (violet) and pathway (yellow) databases.

We thank Olivier Boucherat, Jean Charron and Tom Moss for comments on the manuscript and sharing reagents. We also thank Venkata Manem for guidance in bioinformatic analyses, Stéphane Dubois from the Next-Generation Sequencing Platform, Genomics Center, CRCHUQ-ULaval for help in RNA-seq assays, and Carl St-Pierre from the Cellular Imaging Unit CRCHUQ-ULaval for advice.

Conceptualization: J.G., N.H., K.L.-T. and L.J.; Methodology: J.G., N.H., K.L.-T. and B.F.; Validation: J.G., N.H. and L.J.; Formal analysis: J.G., N.H. and L.J.; Resources: C.D., W.D.F. and L.J.; Data curation: J.G., N.H. and L.J.; Writing—original draft preparation: J.G. and L.J.; Writing—review and editing: J.G., N.H., K.L.-T., B.F., C.D., W.D.F. and L.J.; Supervision: L.J.; Project administration: L.J.; Funding acquisition: L.J. All the authors read and approved the manuscript before submission.

The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the respective institutional ethics committees and by the CHU de Québec-Université Laval research ethics committee (project MP-20-2020-4809, approved: 5 August 2019).

Informed consent was obtained from all subjects involved in the study.

RNA-seq and miRNA-seq data are available through NCBI GEO (accession numbers GSE311205 and GSE311209) and exome data are available through SRA (accession number PRJNA1356925).

This work was supported by the Canadian Institutes of Health Research (FRN169158 to L.J.) and the Cancer Research Society (24217 to L.J.). B.F. holds studentships from NSERC (BESC M) and Fonds de Recherche du Québec en Santé (BF1-325367).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.