GLI1+ cells are a source of repair-supportive mesenchymal cells (RSMCs) during airway epithelial regeneration

Repair-supportive mesenchymal cells (RSMCs) have been recently reported in the context of naphthalene (NA)-induced airway injury and regeneration. These cells transiently express smooth muscle actin (Acta2) and are enriched with platelet-derived growth factor receptor alpha (Pdgfra) and fibroblast growth factor 10 (Fgf10) expression. Genetic deletion of Ctnnb1 (gene coding for beta catenin) or Fgf10 in these cells using the Acta2-Cre-ERT2 driver line after injury (defined as NA-Tam condition; Tam refers to tamoxifen) led to impaired repair of the airway epithelium. In this study, we demonstrate that RSMCs are mostly captured using the Acta2-Cre-ERT2 driver when labeling occurs after (NA-Tam condition) rather than before injury (Tam-NA condition), and that their expansion occurs mostly between days 3 and 7 following NA treatment. Previous studies have shown that lineage-traced peribronchial GLI1+ cells are transiently amplified after NA injury. Here, we report that Gli1 expression is enriched in RSMCs. Using lineage tracing with Gli1Cre−ERT2 mice combined with genetic inactivation of Fgf10, we show that GLI1+ cells with Fgf10 deletion fail to amplify around the injured airways, thus resulting in impaired airway epithelial repair. Interestingly, Fgf10 expression is not upregulated in GLI1+ cells following NA treatment, suggesting that epithelial repair is mostly due to the increased number of Fgf10-expressing GLI1+ cells. Co-culture of SCGB1A1+ cells with GLI1+ cells isolated from non-injured or injured lungs showed that GLI1+ cells from these two conditions are similarly capable of supporting bronchiolar organoid (or bronchiolosphere) formation. Single-cell RNA sequencing on sorted lineage-labeled cells showed that the RSMC signature resembles that of alveolar fibroblasts. Altogether, our study provides strong evidence for the involvement of mesenchymal progenitors in airway epithelial regeneration and highlights the critical role played by Fgf10-expressing GLI1+ cells in this context. Supplementary Information The online version contains supplementary material available at 10.1007/s00018-022-04599-2.

Introduction
The bronchi of the lung represent important interfaces with the outside environment and are as such exposed to many toxins and pathogens. Therefore, the bronchial epithelium undergoes a process of injury repair. The relevant bronchial epithelial cells for the repair process in the bronchi are the non-ciliated secretory cells or club cells [1,2]. These cells express secretoglobin family 1A member 1 (SCGB1A1, a.k.a. Clara cell 10 kDa secretory protein CC10 or CCSP), and they represent progenitors for both club and ciliated cells [3,4]. A large subset of club cells expresses high levels of cytochrome P450 family 2 subfamily f polypeptide 2 (CYP2F2), which is involved in the transformation of naphthalene (NA) into trans-1R-hydroxy-2R-glutathionyl-1,2-dihydronaphthalene, a toxic intermediate that leads to the depletion of club cells in the lung [5][6][7]. In addition, club cells expressing low levels of CYP2F2, called variant club cells, escape the deleterious effect of NA and represent a pool of progenitors. These cells are present at proximity of neuroepithelial bodies (NEBs) and bronchioalveolar duct junctions (BADJs). These cells are critical for the repair process [8][9][10]. The NA injury model is therefore useful to investigate airway regeneration.
Epithelial-mesenchymal interactions orchestrated by fibroblast growth factors, which have been described in the embryonic lung [11], are also at play after injury in the adult lung. It has been proposed that following NA injury in adult mice, WNT7b expressed by surviving ciliated cells acts on airway smooth muscle cells (ASMCs) leading to the increased expression of Fgf10 [12]. FGF10 binds to and activates FGFR2b expressed in variant club cells to control their amplification [12]. Leucine-rich repeat-containing G-protein coupled receptor 6 (Lgr6)-positive cells represent an ASMC subpopulation that enhances repair of the bronchial epithelium following NA administration and also involves a WNT-FGF10 positive feedback loop [13]. Therefore, there is a consensus in the field that ASMCs are a mesenchymal niche for variant club cells following injury.
In a recently published paper, we showed that during lung repair after NA injury, a mesenchymal population of ACTA2+ PDGFRα+ cells emerges in the non-cartilaginous conducting airway niche, which is normally populated by ASMCs [14]. This cell population, which we termed "repair-supportive mesenchymal cells" (RSMCs), is distinct from conventional ASMCs, which have previously been proposed to contribute to epithelial repair [14]. Gene expression analysis on sorted lineage-labeled cells showed that RSMCs express low levels of ASMC marker genes, but high levels of the pro-regenerative marker Fgf10 [14]. Organoid co-cultures demonstrated an enhanced ability for ACTA2+ cell-derived RSMCs in supporting club cell growth in vitro [14,15].
In this study, we carried out lineage tracing of ACTA2+ cells captured before (defined as Tam-NA condition) or after NA injury (defined as NA-Tam condition). We also analyzed the contribution of FGF10 derived from the GLI1+ lineage to the repair process. Bronchiolosphere cultures using GLI1+ mesenchymal cells isolated from noninjured or injured lungs co-cultured with SCGB1A1+ cells were also carried out. Our study demonstrates the critical role played by Fgf10-expressing GLI1+ cells in the context of repair of the bronchial epithelium.

Immunofluorescence
Lungs were flushed by transcardiac perfusion with PBS and then, either fixed with 4% paraformaldehyde for paraffin embedding and 5 μm sectioning using a microtome or inflated with 1.7% low melting agarose for 200 μm sectioning using a vibratome. Immunofluorescence (IF) was performed using monoclonal anti-ACTA2 (Sigma, 1:200), monoclonal anti-SCGB1A1 (Santa Cruz Biotechnology, 1:200), polyclonal anti-KI67 (Thermo Fisher Scientific, 1:200) and polyclonal anti-RFP antibodies (Thermo Fisher Scientific, 1:200; used only when combined with KI67 staining which requires antigen retrieval). Nuclei were stained using DAPI (Life Technologies). Tile scans of thin sections were acquired using EVOS FL2 Imaging System (Thermo Fischer Scientific). Other fluorescent images were acquired using a Leica DM550 B fluorescence microscope, an SP5 or SP8 confocal microscope, both equipped with a white-light laser and hybrid detectors (Leica Microsystems). Where applicable, three-dimensional (3D) reconstruction of z-stacks was performed using LAS X software (version 3.5; Leica Microsystems). For quantification of stained control and experimental lungs, sections from at least 3 independent lungs were analyzed and multiple images were used (n > 8).

Flow cytometry analysis and cell sorting
Flow cytometry and cell sorting were carried out as we previously described [14]. In brief, lungs were isolated in Hank's Balanced Salt Solution (HBSS, Gibco), minced into small pieces and incubated with 0.5% collagenase type IV in HBSS (Life Technologies) at 37 °C for 45 min. Lung homogenates were then passed through 70 and 40 μm cell strainers (BD Biosciences) to obtain single-cell suspensions. Cells were centrifuged at 4 °C at 1000 rpm for 5 min and then resuspended in MACS buffer and stained with anti-PDGFRα (APC-conjugated, 1:100), anti-EpCAM (APC-Cy7-conjugated, 1:100), anti-CD31 (Pacific Blue-conjugated, 1:100) and anti-CD45 (Pacific Blue-conjugated, 1:100) antibodies (all from Biolegend) for 20 min on ice in the dark. Then, cells were washed with MACS buffer. Flow cytometry and cell sorting were done using FACSAria III cell sorter (BD Biosciences).

Quantitative real-time PCR
For bulk RNA extraction from lungs, accessory lobes were lysed and subjected to RNeasy Mini kit (Qiagen) as we previously described [14]. Extraction of RNA from sorted cells was performed using RNeasy Micro kit (Qiagen). Primers were designed using the Universal Probe Library Assay Design center program (Roche Applied Science), and quantitative real-time PCR (qPCR) was performed using Light Cycler 480 II (Roche Applied Science). Hypoxantine-guanine phosphoribosyltransferase (Hprt) was used as a reference gene.

Single-cell RNA sequencing
Single-cell RNA-seq of tdTomato+ (hereafter referred to as tdTom+) cells isolated from Acta2-Cre-ERT2; tdTomato flox lungs was performed using the ICELL8 platform (Takara Bio) using the V2 protocol with in-chip preamplification. Samples were presorted into HBSS (Gibco) at a concentration of 20 cells/µL and stained with Hoechst dye (Invitrogen) for 5 min at RT prior to loading onto the MSND dispenser. For each sample, 3 dispenses were made on one 5,184-well chip resulting in 348 and 459 wells containing single cells after selection in CellSelect software (Takara Bio). Selected wells were processed for RT and in-chip preamplification, and cDNA was pooled for final library preparation using standard Nextera protocol (Illumina). The library was checked in Labchip Gx touch 24 (Perkin Elmer) and sequenced on Nextseq500 (Illumina) using V2 chemistry and paired-end setup (11 bp Read1 for cell barcode and 80 bp Read2 for cDNA template sequencing). A total of 537 M reads were obtained.
Single-cell RNA-seq of tdTom+ cells isolated from Gli1 Cre−ERT2 ; tdTomato flox lungs was carried out using the 10x Genomics Single-Cell 3′ v3.1 RNA-seq kit. Gene expression libraries were prepared according to the manufacturer's protocol. MULTI-seq barcode libraries were retrieved from the samples, and libraries were prepared independently.

Single-cell data analysis
ICELL8 platform: cDNA reads were trimmed for poly(A) tails using BBMAP (unpublished, https:// jgi. doe. gov/ dataand-tools/ bbtoo ls/). Reads of length under 50 bases were discarded. The 456 M remaining reads were then aligned using STAR (release v2.4.0a) against Mus musculus build mm10 following ENCODE RNA-seq mapping recommendations. STAR indices were generated using Ensembl GTF file (release 83). For counting based on read counts, we used the Dropseq Core Computational Protocol version 1.12 (dropseq.jar) [19]. Gene-to-cell count table matrix was then used for statistical analysis performed within the Seurat package (http:// satij alab. org/ seurat/) [20]. The pipeline goes through cell and gene filtering, data normalization, then finding the most variable genes to perform principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) dimension reduction for visualization purpose. We first created a Seurat object including all cells and all genes present in at least 5 different cells. We then filtered out cells with more than 95% dropouts or less than 50,000 read counts. We also discarded cells with more than 20% of counts from mitochondrial genes or 10% counts from ribosomal genes. We then normalized the count values for the 305 remaining cells to a scaling factor corresponding to the median read counts per cell. We next identified most variable genes by plotting genes into bins based on X-axis (average expression) and Y-axis (dispersion) using cutoffs X = 3 and Y = 0.1, identifying 1,212 genes. We then performed a PCA, after scaling and centering of the data across these genes. This ensured robust identification of the primary structures in the data. We identified the first 8 principal components as statistically significant and used it for t-SNE calculation (perplexity = 30). The t-SNE procedure returns a two-dimensional (2D) embedding of single cells utilized for dataset visualization. Cells with similar expression signatures of genes within the variable set, and therefore, similar PC loadings, are most likely to localize near each other in the embedding, and hence, distinct cell types form 2D point clouds across the t-SNE map. Identification of the clusters of cells was determined using the FindClusters Seurat function (resolution = 1). The Seurat FindAllMarkers function was used to identify differentially expressed genes across each cluster (min.pct = 0.25, thresh.use = 0.25).
Chromium 10x platform: Single-cell suspensions were processed using the Chromium Next GEM Single-Cell 3ʹ Reagent Kits v3.1 (10x Genomics). Each sample was run separately on a lane in Chromium controller, and gene expression libraries were prepared according to the manufacturer's protocol. Library sample index sequences were incorporated as the sample index read and sequenced together on a Nextseq2000 sequencer. Raw reads were aligned against the mouse genome (mm10, ensemble assembly 104) and mapped and counted by StarSolo [21]. Secondary analyses were carried out following the Seurat workflow. We removed 285 cells that did not express more than 1,000 genes, with less than 20,000 detected UMIs or had a mitochondrial content above 15%. Raw counts were normalized and scaled using the SCTransform function. The dimension of the count matrix was reduced by PCA. We used the first 40 principal components to compute the UMAP embedding and construct the nearest neighbor graph. Clusters were identified using a resolution of 0.3 and annotated according to the expression of markers described in the literature [22]. Gene expression exploration and differential analyses were carried out on log-normalized data.

Statistical analysis and figure assembly
Quantitative data were assembled and analyzed using Graph-Pad Prism software (GraphPad Software). These data are presented as average values ± standard error of mean (sem). Student's t-test was used to compare 2 groups, one-way ANOVA was used to compare 3 or more groups in the presence of 1 variable (e.g., treatment, time or genotype), and two-way ANOVA was used to compare multiple groups in the presence of 2 variables (e.g., treatment and genotype, time and genotype). Differences were considered significant if P < 0.05. The number of biological replicates is indicated in the corresponding figure legends. Figures were assembled using Adobe Illustrator (Adobe).

Mesenchymal cells transiently expressing Acta2 emerge during epithelial injury and repair
In order to analyze the spatial distribution of the ACTA2+ lineage around the injured airways, including our recently reported RSMC population [14], the Acta2-Cre-ERT2; tdTomato flox lineage-tracing model was employed. Animals were randomly allocated into 3 groups receiving either tamoxifen food prior to corn oil injection (Tam-Oil), tamoxifen food prior to NA injection (Tam-NA) or NA injection followed by tamoxifen food (NA-Tam) (Fig. 1A). Animals were euthanized at day 14 (d14) after NA or Oil injections. Such approach allowed labeling (and visualizing) steady-state ACTA2+ cells, including those around the airways (ASMCs) (Tam-Oil), pre-existing ACTA2+ cells (including ASMCs) during repair after injury (Tam-NA) and total ACTA2+ cells during repair after injury (bona fide SMCs and newly recruited ACTA2+ cells; RSMCs) (NA-Tam). While the lineage label, tdTom, colocalized with the ACTA2 stain in the Tam-Oil (Fig. 1B, E, H, K) and Tam-NA groups (Fig. 1C, F, I, L), imaging of NA-Tam lungs revealed, in addition to ACTA2+ tdTom+ cells in the ASMC layer, the presence of ACTA2− tdTom+ cells between the ASMC layer and the airway epithelium, corresponding to RSMCs (Fig. 1D, G, J, M). These cells were clearly distinct from ASMCs and were mostly observed around the proximal lobular airways. They were more abundant around the mainstem bronchus and at the branch points. Quantification of the distance between the ASMC layer and airway epithelium showed significant thickening of the peribronchial region in Tam-NA, similarly to NA-Tam, compared to the Tam-Oil condition (Fig. 1N), confirming that the additional layer of RSMCs also forms in the Tam-NA condition but cannot be visualized due to the timing of tamoxifen administration (before injury).
The number of tdTom+ cells in these three experimental conditions was quantified by flow cytometry (Fig. 1O). Compared to the number of tdTom+ cells in Tam-Oil, our results demonstrate an increase in NA-Tam compared to Tam-Oil and Tam-NA, confirming that additional cells acquire ACTA2 expression after NA injury. We also investigated Acta2, Pdgfra, Fgf10 and Gli1 expression in sorted tdTom+ cells isolated from Tam-Oil, Tam-NA and NA-Tam lungs. Our results indicate that Pdgfra and Fgf10 expression are indeed upregulated in NA-Tam tdTom+ cells (RSMCenriched) compared to Tam-Oil or Tam-NA conditions (both SMC-enriched), while Acta2 expression showed the opposite trend (Fig. 1P). Gli1 expression showed a trend for upregulation in NA-Tam compared to the other groups (Fig. 1P).
The emergence of tdTom+ RSMCs in NA-Tam, but not Tam-NA, lungs indicated that those cells were descendants of a mesenchymal population that transiently acquired Acta2 expression after NA injury. To study the dynamics of tdTom+ cells that are labeled after injury, NA-Tam lungs were analyzed by flow cytometry at d3, d7 and d14 (Fig.  S1A, B). The results revealed a significant increase in the abundance of tdTom+ cells from d3 to d7 after NA injury (Fig. S1C). Interestingly, such difference in tdTom+ cells was no longer seen between d7 and d14 (Fig. S1C), indicating that RSMCs peaked in terms of abundance at d7.

Single-cell transcriptomic analysis reveals a mesenchymal niche that is associated with epithelial regeneration
We next performed scRNA-seq on tdTom+ cells derived from Acta2-Cre-ERT2; tdTomato flox lungs to delineate the complexity and heterogeneity of these cells, and to identify the RSMC cluster. FACS-sorted tdTom+ cells from Tam-NA (n = 3 animals) and NA-Tam (n = 3 animals) were used to dispense single cells onto an ICELL8 chip (Fig. 2). Following library preparation and RNA-seq, quality control yielded 122 cells from the Tam-NA group and 183 cells from the NA-Tam group ( Fig. 2A-C). Bioinformatic analysis revealed the presence of five cellular clusters corresponding to SMCs, endothelial cells, macrophages and 2 additional clusters that showed matrix fibroblast signatures, one expressing Col13a1 and the other one Col14a1 (Figs. 2A-D,  S2A). Similar matrix fibroblast clusters have been previously described [22,24]. Importantly, the analysis showed a significant emergence of the COL13A1+ matrix fibroblast-like population in NA-Tam lungs compared to Tam-NA lungs (Fig. 2B-D). We have previously reported that deletion of Ctnnb1 and Fgf10 in NA-Tam mice led to repair impairment and that Pdgfra marks the RSMC population in the NA-Tam condition [14]. In the current study, Pdgfra indeed marked the emerging COL13A1+ population corresponding to RSMCs (Figs. 2G and S2A). This cluster expressed typical alveolar fibroblast/lipofibroblast markers such as Tcf21, Limch1 and Pdgfra (Figs. 2G and S2). Comparison between NA-Tam RSMCs and Tam-NA SMC confirmed that Pdgfra was among the top genes discriminating between RSMCs and pre-existing SMCs in the context of NA injury (Fig. S2B  and

The GLI1+ lineage is a source of RSMCs in the lung
The peribronchial and perivascular regions of the lung are populated by hedgehog responsive GLI1+ mesenchymal cells (Fig. S1F). These cells have been shown to contribute to bleomycin-associated myofibroblasts in lung fibrosis [25] and were also shown to be amplified around the airways following NA injury [26]. A proportion of GLI1+ cells is intermingled with ASMCs in the peribronchial space (Fig.  S1F). We previously used the Acta2-Cre-ERT2; tdTomatoflox mice to FACS-isolate the RSMC-enriched population (PDGFRα+ tdTom+) and the other mesenchymal populations (SMC-enriched PDGFRα− tdTom+) from NA-Tam animals at d3, d7 and d14 [14] (Fig, S1D). Gene expression analysis showed that Gli1 is highly expressed in the RSMCenriched fraction, and its expression peaks at d7 (Fig. S1E). Therefore, we suspected that GLI1+ cells might serve as a source of RSMCs following NA injury.
To test the hypothesis that pre-existing GLI1+ cells are a source of RSMCs, Gli1 Cre−ERT2 ; tdTomato flox mice were fed tamoxifen-containing pellets for 2 weeks followed by a 2-week washout period and then treated with NA or Oil (Fig. 3A). This approach allowed labeling of pre-existing GLI1+ cells and studying their fate during epithelial repair after injury. A similar experimental setup was carried out using Gli1 Cre−ERT2 ; tdTomato flox ; Fgf10 flox/flox mice to investigate whether FGF10 specifically derived from these cells is required for the repair process (Fig. 3A). Genetic deletion of both Fgf10 alleles in pre-existing GLI1+ cells significantly impaired club cell replenishment (Fig. 3B, C). Immunofluorescence for KI67 showed a significant decrease in club cell proliferation in experimental (Gli1 Cre−ERT2 ; tdTomato flox ; Fgf10 flox/flox ) vs. control (Gli1 Cre−ERT2 ; tdTomato flox ) lungs following NA treatment (Fig. S3). Quantitative real-time PCR using lung homogenates showed significant downregulation of Scgb1a1 at d14 after NA treatment in experimental mice compared to controls (Fig. 3D). Interestingly, we did not observe a significant drop in Fgf10 expression in non-injured (Oil condition) experimental vs. control lungs, suggesting that at the resting state, GLI1+ cells are not the main source of Fgf10 expression in the lung (Fig. 3D). However, following injury, a clear drop in Fgf10 expression is observed in experimental vs. control lungs fitting with the proposed increase in Fgf10 expression in the lung during the repair process (Fig. 3D). Of note, we could not detect significant proliferation in the surrounding mesenchyme in the NA vs. Oil-treated control lungs at d14 (Fig. S3), which is in line with the results shown in Fig. S1C, where no change in tdTom+ cells was seen between d7 and d14, thus supporting the notion that RSMCs were no longer expanding after d7.
To verify if the increase in Fgf10 expression (Fig. 3D) is linked to Fgf10 upregulation in GLI1+ cells, we analyzed Fgf10 expression in sorted GLI1+ cells. Our results indicate no significant difference in Fgf10 expression between the d14 group and the Oil group (Fig. 4E). Therefore, the global increase in Fgf10 expression is likely linked to a quantitative increase in GLI1+ cells. Furthermore, qPCR on tdTom+ cells sorted from Oil-and NA-treated Gli1 Cre−ERT2 ; tdTomato flox lungs revealed a transient downregulation of Gli1 at d3 as previously described [26], in parallel to a significant, transient upregulation of Acta2 at d3 (Fig. 4E). Flow cytometry-based quantification confirmed the acquisition of Acta2 expression in the GLI1+ lineage upon NA injury (Fig. 4F). Altogether, these data suggest that pre-existing GLI1+ cells 1) transiently acquire Acta2 expression and give rise to RSMCs during epithelial regeneration (and are therefore captured when ACTA2+ cells are labeled after injury (Fig. 1)), 2) are a source of FGF10 during both steady state and after injury and 3) are massively amplified and recruited to the airways in response to NA injury.

GLI1+ cells support club cell growth in an in vitro co-culture organoid system
We then opted to test whether GLI1+ cells derived from d14 Oil-and NA-treated lungs possess an intrinsic capacity to support the growth of club cells in a co-culture organoid system. Club cells sorted using uninjured Scgb1a1 CreERT ; tm14 tdTomato flox mice or epithelial progenitors sorted from wild-type mice were co-cultured with tdTom+ cells derived from Tam-Oil or Tam-NA Gli1 Cre−ERT2 ; tdTomato flox mice (Fig. 4G). Representative images of large bronchiolospheres are shown in Fig. 4H, I. Further analysis of these bronchiolospheres indicated a slight trend for increased organoid number and size in the NA vs. Oil condition (Fig. 4J, K), suggesting that the beneficial effect of GLI1+ cells appearing after injury is likely linked to their quantity/recruitment rather than a qualitative change. Additionally, loss of function of Fgf10 in GLI1+ cells led to decreased organoid size compared to controls (Fig. S4).

Single-cell RNA-seq of GLI1+ cells suggests that alveolar fibroblasts, co-expressing Gli1 and Fgf10, converge on the RSMC lineage
The organoid assay suggested that GLI1+ cells have an intrinsic club cell-supportive potential even in the absence of injury (Fig. 4G-K). We therefore explored the heterogeneity of GLI1+ cells in non-injured lungs using scRNAseq (Fig. 5). UMAP analysis showed the presence of seven clusters (Fig. 5A). Cluster 1 corresponds to alveolar fibroblasts, while cluster 6 represents a subcluster of alveolar fibroblasts expressing high levels of Mmp3 (Fig. 5A, B). Clusters 0 and 2 contain two adventitial fibroblast subsets. Clusters 3 and 4 include peribronchial fibroblasts, while cluster 5 corresponds to SMCs. Examination of the RSMC transcriptomic signature defined in Fig. 2 shows that most of the signature genes are contained mostly in alveolar fibroblasts (clusters 1 and 6), and to a less extent in adventitial fibroblasts (Fig. 5B lower panel). Interestingly, cluster 5 (SMCs) did not display significant RSMC signature, thus supporting the notion that RSMCs are distinct from SMCs. Given the important role played by Fgf10 expressed by GLI1+ cells in the repair process and the fact that Fgf10 is not significantly upregulated in these cells after NA injury, we reasoned that GLI1+ FGF10+ cells exist in the non-injured lung. This is rather counter-intuitive as Gli1, a downstream target of sonic hedgehog signaling, and Fgf10 expression are thought to be mutually exclusive especially during lung development [27]. Confirming our hypothesis, these cells were found in alveolar fibroblasts and adventitial fibroblasts but completely absent in the SMCs (Fig. 5C). Quantification of the expression of Gli1, Fgf10 and Acta2 confirmed that alveolar fibroblasts represent the   (Fig. 5D). Using our novel Fgf10 Cre−ERT2 mice recently reported [28], we also examined the localization of FGF10+ cells in non-injured lungs. In addition to being located in the alveolar region [28], FGF10+ cells were also observed around the bronchi and distinct from ASMCs (Fig. 5E), similarly to what was observed for GLI1+ cells (Fig. S1F).

Discussion
In this study, we demonstrated that GLI1+ cells are a source of RSMCs during airway epithelial regeneration. Lineage tracing of ACTA2+ cells (using the Acta2-Cre-ERT2 driver line) labeled either before (Tam-NA) or after (NA-Tam) injury indicates that these cells appear mostly following injury. Nevertheless, lineage-tracing experiments using the Gli1 Cre−ERT2 driver line indicate that RSMCs also arise from pre-existing GLI1+ cells that acquire Acta2 expression and accumulate in the peribronchial region. Deletion of Fgf10 in GLI1+ cells leads to impaired repair. Using bronchiolosphere assays, we show that GLI1+ cells isolated from either non-injured or NA-treated lungs similarly support club cell growth in vitro. Overall, our data indicate that GLI1+ cellderived RSMCs produce FGF10 that is required for the repair process after NA injury. As such, RSMCs represent a privileged mesenchymal niche for club cell progenitors.

The lineage map of the lung mesenchyme is still poorly understood
Despite extensive research, the lineage map of the lung mesenchyme is still poorly understood, with the best-characterized cell types being ASMCs/VSMCs and pericytes, and, to a less extent, alveolar fibroblasts (mostly during neonatal stages) and myofibroblasts. The perivascular region of the lung is populated by GLI1+ mesenchymal cells that have been proposed to constitute a subpopulation of PDGFRβ+ pericytes [25]. These cells have been shown to possess mesenchymal stem cell (MSC)-like characteristics such as trilineage differentiation potential in vitro (giving rise to adipocytes, osteoblasts and chondrocytes) and also represent a repertoire of myofibroblast progenitors in fibrotic disease [25,29]. On the other hand, the peribronchial domain of the postnatal mouse lung consists of cartilaginous rings (around main-stem bronchi), ASMCs (around conducting airways) and other poorly characterized mesenchymal cells, some of which also express Gli1.

RSMCs are a component of the niche required for activating club cell progenitors
In this study, we demonstrate the complexity of the mesenchymal niche in the mature lung and identify GLI1+ cellderived RSMCs as a component of the niche required for activating club cell progenitors. GLI1+ cells are amplified, transiently acquire Acta2 expression and contribute to the repair process by producing FGF10 (Fig. 6). We previously published side-by-side comparison of PDGFRα+ tdTom+ cells (RSMC-enriched population) and PDGFRα− tdTom+ cells (SMC-enriched population), indicating that RSMCs display higher WNT signaling activation (using Axin2 mRNA levels as readout), higher Fgf10 expression and greater potential in supporting club cell growth in an in vitro organoid system [14]. However, one limitation was that the PDGFRα− tdTom+ population contains a mixture of bona fide ASMCs/VSMCs and other subsets such as matrix fibroblasts rather than a pure population of ASMCs, and this might reduce the colonyforming efficiency and consequently bronchiolosphere number and size. Moreover, we could not exclude the possibility that parenchymal lineage-labeled cells (such as alveolar fibroblasts) also contribute to repair in the conducting airway epithelium. Our data strongly indicate that GLI1+ cells, which are normally found in the peribronchial, perivascular and alveolar regions, serve as a source of FGF10+ RSMCs in the lung.

How different are RMSCs from the resident mesenchymal cell (rMC) niche?
Resident mesenchymal cells supporting the growth of alveolar stem cells were further characterized [30]. Using flow cytometry to isolate cells as well as the alveolosphere model, Fig. 6 Model summarizing the role of RSMCs as a niche that facilitates epithelial regeneration after naphthalene injury. GLI1+ cells are distinct from ASMCs. After naphthalene injury, GLI1+ cells transiently acquire Acta2 expression and get amplified, giving rise to RSMCs and contributing to epithelial regeneration by producing FGF10. Pre-existing ASMCs have also been shown to mediate the repair process via a WNT-FGF10-mediated mechanism [12,13]. CGRP calcitonin gene-related peptide, CYP2F2 cytochrome P450 2F2, NE neuroendocrine, TUB tubulin CD31− CD45− EpCAM− cells (defined as rMC) could be subdivided into SCA1+ and SCA1−. Only SCA1+ rMCs displayed niche activity for AT2 stem cells. Using reporter lines for Fgf10 expression, these cells were further subdivided into SCA1+ FGF10+ rMCs and SCA1+ FGF10− rMCs, with the FGF10+ subset displaying the niche activity. These cells are also positive for PDGFRα.
We also recently used scRNA-seq to further characterize the different subpopulations present in SCA1+ rMCs [31]. We identified eighteen different clusters. Among them, five different clusters matched different subsets of lipofibroblasts and specifically expressed Fgf10, Fgf7 and Pdgfra. These clusters could also be differentiated via the level of Col13a1, Col1a1 and Col14a1 expression. One of these five clusters expresses high levels of Col13a1, Fgf10, Limch1, inmt, Npnt and Pdgfra, and could therefore be very similar to RSMCs. Interestingly, this cluster expresses lower levels of Sca1 compared to the other clusters. This puzzling similarity between RSMCs around the bronchi and rMC niche cells in the alveolar region raises the possibility that rMC niche cells could be recruited from the adjacent alveolar region surrounding the bronchi. Further investigation using high-resolution imaging and more specific lineage tracing using the Dre/Cre technology should allow answering this important question.
Author contributions SB and EEA contributed to the study conception and design. Material preparation, data collection and analysis were performed by XC, AL, AM, VK, AIVA, JK, AK, GK, MSF, MZ, IS, IA, SG, KL, MT, SB and EEA. AG, TB, BM, CS, XL, WS, SH, JSZ, SB and EEA contributed to platforms and methodology. The first draft of the manuscript was written by SB and EEA, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Data availability Primary data can be made available upon reasonable request. Single-cell data have been deposited on gene expression omnibus (GEO) GSE215094.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Ethical approval Animal experiments were conducted at JLU Giessen after approval by the Regierungspräsidium Giessen (protocol # G 2/2016) and at Wenzhou Medical University (protocol # 2021-0017).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.


s2fieldsofstudy: Biology Chemistry

text:
Antiparasitic Ovalicin Derivatives from Pseudallescheria boydii, a Mutualistic Fungus of French Guiana Termites

Social insects are in mutualism with microorganisms, contributing to their resistance against infectious diseases. The fungus Pseudallescheria boydii SNB-CN85 isolated from termites produces ovalicin derivatives resulting from the esterification of the less hindered site of the ovalicin epoxide by long-chain fatty acids. Their structures were elucidated using spectroscopic analysis and semisynthesis from ovalicin. For ovalicin, these compounds displayed antiprotozoal activities against Plasmodium falciparum and Trypanosoma brucei, with IC50 values of 19.8 and 1.1 µM, respectively, for the most active compound, i.e., ovalicin linoleate. In parallel, metabolomic profiling of a collection of P. boydii strains associated with termites made it possible to highlight this class of compounds together with tyroscherin derivatives in all strains. Finally, the complete genome of P. boydii strains was obtained by sequencing, and the cluster of potential ovalicin and ovalicin biosynthesis genes was annotated. Through these metabolomic and genomic analyses, a new ovalicin derivative named boyden C, in which the 6-membered ring of ovalicin was opened by oxidative cleavage, was isolated and structurally characterized.

Introduction
Insects might be the most represented and diverse living creatures on earth [1]. In nature, insects live in symbiosis with bacteria and fungi to protect themselves from pathogens [2]. These microorganisms play a key role in the ecological success of their host [3], providing them antipathogenic metabolites that could eventually find applications in therapeutics [4]. Insect/microorganism associations are well described for social insects, such as ants, bees, and some wasp species [5,6]. However, little work has been conducted thus far on termite/microorganism interactions (outside the context of trophobioses) [7][8][9][10].
Termites (Blattodea: Isoptera) are eusocial insects comprising more than 3000 species, mostly distributed in tropical and subtropical regions. All termites live in nests containing 1 to 1.5 million individuals divided into castes (reproductive individuals, workers, and soldiers) that cooperatively perform colony tasks. These insects play an important role in ecosystems through the decomposition of plant lignocellulose, having a positive impact in nutrient cycling and soil-water dynamics through their nest construction and tunnel networks in the soil [11]. Termites are vital members of old-growth tropical forests, being perhaps the main decomposers of dead plant material at all stages of humification (decay). Therefore, termite colonies viability is considered a marker of human influence on the environment [12,13]. Nests protect termites from aggressions of biotic and abiotic origin In parallel, new methodologies for dereplication of complex mixtures have been developed to accelerate and expand the discovery of new bioactive natural molecules. In order to avoid focusing research efforts on already known molecules, dereplication strategies have been developed [29,30]. Dereplication allows for the annotation of known molecules in a mixture but also allows for targeting and focusing on potentially novel molecules to be isolated and characterized. Mass spectrometry-based molecular networking approaches can be used to overcome these difficulties, allowing the targeting of specific metabolites through the tandem mass spectrometry data set organization they confer [29][30][31][32][33][34][35]. In this study, a molecular network approach allowed us to navigate the chemical space of P. boydii metabolomes and was used to annotate the chemical diversity of these fungi. A new compound was identified by this method and was isolated and fully characterized. Then, it was chosen to compare the metabolomic and genomic data of all P. boydii in our collection. For this purpose, the whole genome sequence of the four P. boydii strains was established, and the cluster of potential genes for biosynthesis of ovalicin and its fatty acid derivatives was annotated.

Results
The study of the P. boydii SNB-CN85 ethyl acetate extract led to the isolation of ovalicin (1) and the five ovalicin analogs 2-6 ( Figure 1).
Preliminary inspection of the 1 H and 13 C NMR spectra confirmed that Compounds 2-5 were ovalicin analogs (Figures S1-S29 in the Supplementary Materials, Table 1, Tables S1-S5 in the Supplementary Materials). In fact, Compounds 2-5 exhibited the same 1 H and 13 C NMR shifts as ovalicin, and the only difference was observed at position 14 with protons and carbons at δH 4.11 and 4.17/δC 69.4-69.5 ppm instead of δH 2.63 and 2.97/δC 51.0 ppm for ovalicin. For Compound 6, the same chemical shift but also the loss of the Omethyl group at position 7 was observed (Table S5 in    Preliminary inspection of the 1 H and 13 C NMR spectra confirmed that Compounds 2-5 were ovalicin analogs (Figures S1-S29 in the Supplementary Materials, Table 1, Tables S1-S5 in the Supplementary Materials). In fact, Compounds 2-5 exhibited the same 1 H and 13 C NMR shifts as ovalicin, and the only difference was observed at position 14 with protons and carbons at δ H 4.11 and 4.17/δ C 69.4-69.5 ppm instead of δ H 2.63 and 2.97/δ C 51.0 ppm for ovalicin. For Compound 6, the same chemical shift but also the loss of the O-methyl group at position 7 was observed (  13 C, HSQC spectra and 1 H-1 H COSY and 1 H-13 C HMBC correlations indicated the presence of an octadecadienoate chain, as shown by the presence of four olefinic methines (4H at δ H 5.33-5.36), 12 methylenes (δ H 1.28-2.78), and one triplet methyl (δ H 0.91). The chain was also determined following the sequence of 1 [36,37]. Due to too small quantities, this procedure was not possible for Compounds 5 and 6. To confirm the structures and configurations of Compounds 2, 3, and 4, ovalicin esters were synthesized from Compound 1 by esterification with linoleic acid, oleic acid, and stearic acid, respectively ( Figure S34 in the Supplementary Materials). The comparison of the NMR spectra obtained between the isolated and synthetic compounds allowed us to confirm the identification and configurations of fungal metabolites 2-4 ( Figures S35-S43 in the Supplementary Materials). The absolute configurations Compounds 2-4 were also confirmed first by comparison with circular dichroism data of 1 ( Figure S33 in the Supplementary Materials) for the ovalicin moiety. Furthermore, all these compounds showed negative optical rotation. Secondly, for the fatty acid moiety, our 13 C NMR analyses were in accordance with literature data concerning the specific chemical shifts of methylene at C-11 position between two double bonds ( Table S6 in the Supplementary Materials) [38,39]. lowing the sequence of 1 H-1 H COSY cross-peaks H-8′/H-9′/H-10′/H-11′/H-12′/H-13′. The fatty acid chain was linked to C-14 thanks to the H-14/C-1′ correlation ( Figure 2). Compounds 3-4 and 6 presented similar spectral data. The only difference was observed in the fatty acid chains. The exact mass of each molecule, the presence of olefinic protons between δH 5.33 and 5.36 ppm, and the 1 H-1 H COSY correlations confirmed the presence of two double bonds for 2 and 6, one double bond for 3 and 5, and none for 4 ( Figure 2).   [36,37]. Due to too small quantities, this procedure was not possible for Compounds 5 and 6. To confirm the structures and configurations of Compounds 2, 3, and 4, ovalicin esters were synthesized from Compound 1 by esterification with linoleic acid, oleic acid, and stearic acid, respectively ( Figure S34 in the Supplementary Materials). The comparison of the NMR spectra obtained between the isolated and synthetic compounds allowed us to confirm the identification and configurations of fungal metabolites 2-4 ( Figures S35-S43 in the Supplementary Materials). The absolute configurations Compounds 2-4 were also confirmed first by comparison with circular dichroism data of 1 ( Figure S33 in the Supplementary Materials) for the ovalicin moiety. Furthermore, all these compounds showed negative optical rotation. Secondly, for the fatty acid moiety, our 13 C NMR analyses were in accordance with literature data concerning the specific chemical shifts of methylene at C-11′ position between two double bonds ( Table S6 in the Supplementary Materials) [38,39].
Compounds 1-4 antiplasmodial activity against Plasmodium falciparum and antitrypanosomal activity against Trypanosoma brucei brucei are summarized in Table 2. Ovalicin 1 and fatty acid ester of ovalicin 2 and 4 showed high inhibition percentages at IC50 values Compounds 1-4 antiplasmodial activity against Plasmodium falciparum and antitrypanosomal activity against Trypanosoma brucei brucei are summarized in Table 2. Ovalicin 1 and fatty acid ester of ovalicin 2 and 4 showed high inhibition percentages at IC 50 values below 20.6 µM against Plasmodium falciparum. Against Trypanosoma brucei brucei, all compounds exhibited IC 50 values below 5 µM, but 1 to 2 were the only compounds to show activity (IC 50 < 1.1 µM). Compounds 5 and 6 were isolated in too small amounts and could not be tested. To the best of our knowledge, only ovalicin is known to show antiangiogenic properties [40] and potent antiplasmodial activity [41]. Consequently, the isolated compounds were also tested for their cytotoxicity on normal primary human umbilical vein endothelial cells (HUVECs). All compounds demonstrated low cytotoxicity on the tested cell lines (IC 50 > 20 µM), contrary to reported literature for ovalicin with an IC 50 of 13 µM [42,43]. Finally, isolated compounds were also tested for their potential antimicrobial activities against C. albicans. The minimal inhibitory concentrations (MICs) observed for molecules 1-3 were greater than 128 µg/mL, showing no significant inhibition. As it has been previously shown that strains of P. boydii were able to biosynthesize derivatives of tyrocherin and ovalicin [20][21][22], an evaluation was conducted of whether our different strains were also capable of biosynthesizing the fatty acid derivatives of ovalicin isolated from SNB-CN85. For this purpose, EtOAc extracts of four P. boydii strains, SNB-CN71, -CN73, -CN81, and -CN85, were analyzed by reversed-phase liquid chromatography-high-resolution mass spectrometry (RPPLC-HRMS/MS) using the datadependent acquisition mode. The resulting spectral data were further investigated by the molecular network approach (MN) [29][30][31]. Data were processed by MZmine 2.51 [32,33], and an MN was generated by Metgem [34]. The resulting network grouped 3021 nodes into 83 clusters ( Figure 3 and Figure S44 in the Supplementary Materials). Strain mapping was applied by attributing a given color code to each strain in the collection, allowing us to visually navigate through the MN and check for the presence and distribution of specific compound classes within the collection. Relative quantification of the ions was represented using pie charts, with their proportions based on the respective areas of the peaks detected on extracted ion chromatography (XIC).
Molecules 2022, 27, x FOR PEER REVIEW  8 of 21 tions, although they were grown under the same conditions. Compound 6 was not detected in the cluster, likely because this compound is a demethylated derivative of ovalicin but mostly because it seems to be present in very small quantities in the crude extract. Following these results, the biosynthetic pathway of ovalicin was investigated. For this purpose, whole genome sequencing of the four P. boydii strains SNB-CN71, -CN73, -CN81, and -CN85 was performed by combining short and long sequences to achieve de novo hybrid assemblies. The quality of the four 42-43 Mb genomes was assessed by the BUSCO score (95.9-96.1% completeness) ( Table S6 in the Supplementary Materials) [44]. It was considered that ovalicin is a known fumagillol analog, an intermediate in the fumagillin biosynthetic pathway [45,46]. While the fumagillin biosynthetic pathway has been described [47], this is not the case for ovalicin and its fatty ester derivatives. Furthermore, β-trans-bergamotene is a biosynthetic intermediate that ovalicin shares with fumagillin and fumagillol [45,46]. When our four genomes were analyzed using the antiSMASH pipeline [48], a cluster of biosynthetic genes common to all four strains, comprising enzymes responsible for producing ovalicin and its fatty acid esters, was identified. This cluster, composed of 18 genes conserved in the four P. boydii strains, was intertwined with the pseurotin A biosynthetic gene cluster (6 genes out of 18). Pseurotin A is another specialized metabolite produced by P. boydii strains whose protonated species were identified in the molecular networks ( Figure 3, Figure S44 in the Supplementary Materials). These two gene clusters appear to be frequently intertwined, as the pseurotin/fumagillin gene , were present in a large relative proportion in SNB-CN71, -CN73, and -CN81. N-Methyltyrocherin and pseurotin A were not detected in SNB-CN85 and SNB-CN71, but they were previously isolated from these strains. Therefore, pseudallicins A-D isolated previously from SNB-CN85 extract strains were not found during this analysis. This may be because the concentration of these compounds in the crude extracts appeared to be very low, whereas they had been concentrated in the previously studied fractions. Sodium-cationized ovalicin and Compounds 1-5 were detected in a specific cluster present in all extracts. By studying the fragmentation spectra of the compounds grouped in the ovalicin fatty acid cluster, it was possible to determine the presence of ovalicin analogs bearing C18 fatty acid chains with three double bonds (Compound 9), one to two epoxides (Compounds 10, 11 and 12), or with a C16 fatty acid chain (Compounds 5 and 8). Compound 11 exhibited both an epoxide and a double bond ( Figure S45 in the Supplementary Materials). It is interesting to note that this collection of fatty acid esters of ovalicin is quite large, and some of these specialized metabolites are present in the different species of the collection. Finally, not all strains produced the different derivatives in the same proportions, although they were grown under the same conditions. Compound 6 was not detected in the cluster, likely because this compound is a demethylated derivative of ovalicin but mostly because it seems to be present in very small quantities in the crude extract.
Following these results, the biosynthetic pathway of ovalicin was investigated. For this purpose, whole genome sequencing of the four P. boydii strains SNB-CN71, -CN73, -CN81, and -CN85 was performed by combining short and long sequences to achieve de novo hybrid assemblies. The quality of the four 42-43 Mb genomes was assessed by the BUSCO score (95.9-96.1% completeness) ( Table S6 in the Supplementary Materials) [44]. It was considered that ovalicin is a known fumagillol analog, an intermediate in the fumagillin biosynthetic pathway [45,46]. While the fumagillin biosynthetic pathway has been described [47], this is not the case for ovalicin and its fatty ester derivatives. Furthermore, β-trans-bergamotene is a biosynthetic intermediate that ovalicin shares with fumagillin and fumagillol [45,46]. When our four genomes were analyzed using the antiSMASH pipeline [48], a cluster of biosynthetic genes common to all four strains, comprising enzymes responsible for producing ovalicin and its fatty acid esters, was identified. This cluster, composed of 18 genes conserved in the four P. boydii strains, was intertwined with the pseurotin A biosynthetic gene cluster (6 genes out of 18). Pseurotin A is another specialized metabolite produced by P. boydii strains whose protonated species were identified in the molecular networks ( Figure 3, Figure S44 in the Supplementary Materials). These two gene clusters appear to be frequently intertwined, as the pseurotin/fumagillin gene cluster has also been reported in Aspergillus fumigatus [49]. The eight genes coding for enzymes involved in the ovalicin biosynthetic pathway were named Pbo_Ova_A-H ( Figure 4, Table S7 in the Supplementary Materials), and the six genes coding for pseurotin A were named Pbo_Pse_A-F (Table S8 in the Supplementary Materials). Their functions were predicted using the Pfam and UniProt databases [50,51]. Like fumagillin, the ovalicin biosynthetic pathway is initiated by the action of the terpene synthase Pbo_Ova_A, which converts farnesyl pyrophosphate to β-trans-bergamoten (Btb). A cytochrome P450 monooxygenase Pbo_Ova_B whose af510 is the closest homolog (85% similarity) was predicted to catalyze the formation of intermediate (Btb1) Figure S46 in the Supplementary Materials). Finally, fatty acids were likely added to dihydroovalicin via an acyltransferase, the most likely candidate of which was Pbo_Ova_G, whose gene was present at the edge of the Pbo_Ova/Pse biosynthetic gene cluster.  This biosynthetic scheme was fully consistent with experimental data proposed previously, concerning configurations of compounds 2-6 and sheds light on subunits' ovalicin absolute configurations.
To more deeply explore the chemical diversity of the P. boydii metabolome, we looked for nodes with intense signals but could not annotate them by searching standards and analogs in MetGem databases. It was chosen to isolate the specialized metabolites showing an ion signal at m/z 323.1433 (Figure 3 and Figure S44 in the Supplementary Materials), which could be a potential original derivative of ovalicin in view of the observed mass. The molecular formula of Compound 7 was determined to be C 15 H 24 O 6 by ESI-HRMS (m/z 323.1473 [M + Na] + , calcd. for C 15 H 24 O 6 Na + , 323.1471, err 1.3 ppm), indicating the presence of two additional oxygen and the lack of one carbon compared with the molecular formula of ovalicin (1). The 1 H NMR spectra of Compounds 1 and 7 displayed some similarities, suggesting that 7 was an analog of 1. The remarkable difference between Compounds 1 and 7 in 1 H NMR was the absence, in 7, of the singlet methoxy signal integrating for 3H at δ H 3.49 (Table 3).
However, Compound 7 exhibited characteristic dimethylallyl signals at δ H 5.24 (m, 1H), δ H 1.73 (m, 3H), and 1.67 (m, 3H). The examination of the COSY correlations allowed us to link this dimethylallyl moiety to H-8 at 3.26, (t, J = 6.5 Hz, 1H) [8]. The complete proton assignment of Compound 7 was deduced from the careful examination of COSY, HSQC, and HMBC correlations ( Figure 5). HMBC analysis indicated that methylene protons at position 2 at δ H 2.22 ppm correlated with a carbonyl carbon of a carboxylic acid at δ C 178.0 ppm. The analysis of the HSQC correlations showed that Compound 7 had additional oxygenated methylene proton signals at δ H 3.64 and 3.60 ppm. Consequently, we deduced that Compound 7 was an analog of demethylated ovalicin in which the C-1-C-7 bond was cleaved by oxidation or hydrolysis. The relative and absolute stereochemistry at C-7 was defined as drawn in Figure 5 based on biosynthetic considerations (Figure 4). Compound 7 was named boyden C. Unfortunately, as this compound was isolated in small quantities, it could not be biologically tested.

General Experimental Procedures
Optical rotations ([α]D) were measured on an Anton Paar MCP 300 polarimeter in a 100 mm long 350 µL cell (Anton Paar, Graz, Austria). IR spectra were recorded using a Perkin Elmer Spectrometer BX FT-IR (Thermo Fisher scientific, Les Ulis, France). Nuclear magnetic resonance (NMR) spectra were recorded on a Bruker 600 MHz spectrometer equipped with a 3 mm inverse detection cryoprobe or a Bruker 500 and 800 MHz spectrometer equipped with a 5 mm inverse detection probe (Bruker, Rheinstetten, Germany). Chemical shifts (δ) are reported as ppm based on the TMS signal. HR-ESI-MS measurements were performed using a Waters Acquity UPLC system with column bypass coupled to a Waters LCT Premier time-of-flight mass spectrometer equipped with an electrospray interface (ESI) (Waters, Manchester, England). Flash chromatography was performed on a Grace Reveleris system with dual UV and ELSD detection equipped with a 40 g C18 column. Analytical and preparative HPLCs were conducted with a Gilson system equipped with a 322 pump device, a GX-271 fraction collector, a 171 diode array detector, and a prepELSII electrospray nebulizer detector (Gilson, Middelton, USA). Columns used for these experiments included a Phenomenex Luna C18 5 µm 4.6 × 250 mm analytical column and a Phenomenex Luna C18 5 µm 21.2 × 250 mm preparative column. The flow rate was set to 1 or 21 mL/min using a linear gradient of H2O mixed with an increasing proportion of CH3CN. Both solvents were HPLC grade and modified with 0.1% formic acid. Potato dextrose agar (PDA) was purchased from Fluka Analytical (Janssen Pharma-

General Experimental Procedures
Optical rotations ([α]D) were measured on an Anton Paar MCP 300 polarimeter in a 100 mm long 350 µL cell (Anton Paar, Graz, Austria). IR spectra were recorded using a Perkin Elmer Spectrometer BX FT-IR (Thermo Fisher scientific, Les Ulis, France). Nuclear magnetic resonance (NMR) spectra were recorded on a Bruker 600 MHz spectrometer equipped with a 3 mm inverse detection cryoprobe or a Bruker 500 and 800 MHz spectrometer equipped with a 5 mm inverse detection probe (Bruker, Rheinstetten, Germany). Chemical shifts (δ) are reported as ppm based on the TMS signal. HR-ESI-MS measurements were performed using a Waters Acquity UPLC system with column bypass coupled to a Waters LCT Premier time-of-flight mass spectrometer equipped with an electrospray interface (ESI) (Waters, Manchester, England). Flash chromatography was performed on a Grace Reveleris system with dual UV and ELSD detection equipped with a 40 g C18 column. Analytical and preparative HPLCs were conducted with a Gilson system equipped with a 322 pump device, a GX-271 fraction collector, a 171 diode array detector, and a prepELSII electrospray nebulizer detector (Gilson, Middelton, WI, USA). Columns used for these experiments included a Phenomenex Luna C18 5 µm 4.6 × 250 mm analytical column and a Phenomenex Luna C18 5 µm 21.2 × 250 mm preparative column. The flow rate was set to 1 or 21 mL/min using a linear gradient of H 2 O mixed with an increasing proportion of CH 3 CN. Both solvents were HPLC grade and modified with 0.1% formic acid. Potato dextrose agar (PDA) was purchased from Fluka Analytical (Janssen Pharmaceuticalaan, Geel, Belgium). Molecular analyses were performed externally by BACTUP, Saint-Priest, France. MS/MS experiments were performed by direct injection analysis (DIA, 10 µL/min) on an Agilent Q-ToF 6540 (APCI interface) mass spectrometer (Agilent, Santa Clara, CA, USA). Compounds were dissolved in CH 2 Cl 2 at 2 mg/mL and then diluted in CH 3 CN at 200 µg/mL. The corona current was fixed to 2 µA, while the desolvation gas flow was set to 8 L/min (heated at 300 • C). The capillary, fragmentor, and skimmer voltages were set to 3000 V, 100 V, and 45 V, respectively. A collision energy of 15 eV was used.

Collection and Identification of Pseudallescheria Boydii
Termite workers were collected from an aerial termite nest located in Rémire-Montjoly, French Guiana, in July 2011. The termites were identified as Termes cf. hispaniolae and Nasutitermes corniger by Prof. Reginaldo Constantino (University of Brasilia, Brasilia, Brazil). The workers were surface-sterilized by successive soakings in 70% EtOH (2 min), 5% NaOCl (2 min), and sterile water. The termite was subsequently placed in a Petri dish containing solid PDA medium. After 1 week at 25 • C, the first fungal hyphae to emerge from the insect were sampled and transferred into other Petri dishes. One of the microbial colonies consisted of a pure fungus, which was saved in triplicate at −80 • C in H 2 O−glycerol (50/50). A sample was submitted for amplification of the nuclear ribosomal internal transcribed spacer region ITS. Sequencing allowed for strain identification by NCBI sequence comparison (Blastn ® , Rockville, MD, USA). The sequences have been registered in the NCBI GenBank database (http://www.ncbi.nlm.nih.gov, accessed on 20 December 2021) under the registry numbers KJ023743 (SNB-CN85), KC684883 (SNB-CN73), KJ023735 (SNB-CN71), and KJ023720 (SNB-CN81).

Culture, Extraction, and Isolation
The P. boydii strains were cultivated on PDA at 26 • C for 15 days, initially on a small scale. The fungus and culture medium were then transferred into a large container and macerated with EtOAc for 24 h. The organic solvent was then collected by filtration, washed with H 2 O in a separatory funnel and evaporated, yielding the four crude extracts.
The P. boydii strain SNB-CN85 was cultivated on a large scale in 150 × 14 cm Petri dishes, yielding 3.57 g extract. A portion of the extract (1.13 g) was subjected to reversedphase flash chromatography using a gradient of H 2 O mixed with an increasing proportion of CH 3 CN to afford 6 fractions (A to E). The column was then eluted with CH 3 CN/CH 2 Cl 2 1:1 (Fraction F). Fraction D (300 mg, eluted with H 2 O/CH 3 CN, 2:8) was purified by reversedphase flash chromatography using the same solvent system as mentioned above with a modified gradient. The major compound was collected in pure form and was identified as ovalicin (1) Figure S34).

Esterification
The free fatty acids were dissolved in MeOH (500 µL). Each solution was cooled down to 0 • C. SOCl 2 (20 µL, 275.4 µmol) was added, and the solution was stirred for 1 h at 50 • C. The solution was concentrated under reduced pressure, and H 2 O (500 µL) was added. The aqueous solution was extracted with a mixture of Et 2 O/heptane, 5:95 (2 mL). The organic layers were filtered on a Na 2 SO 4 /SiO 2 column and concentrated to afford fatty acid methyl esters.  HCl (0.1 mL), the aqueous solution was extracted with Et 2 O (2 × 2.0 mL). The combined organic layers were washed with water (2 mL), filtered, and evaporated to afford a crude mixture (32.2 mg), which was subjected to preparative HPLC (Luna C18, isocratic 5:95 elution over 20 min) to afford synthetic ovalicin linoleate (2)

Identification of the Double-Bond Position in Fatty Acid Methyl Esters by Liquid Chromatography/Atmospheric Pressure Chemical Ionization Mass Spectrometry
The double-bond positions in the fatty acid methyl esters synthetized from ovalicin linoleate (2) and ovalicin oleate (3)  Compounds were eluted at a flow rate of 0.4 mL·min −1 with a gradient from 5% B to 100% B in 20 min and then 100% B for 3 min. The injection volume was fixed at 5 µL for all analyses. Mass spectra were recorded in electrospray ionization positive ion mode with the following parameters: gas temperature 325 • C, drying gas flow rate 10 L·min −1 , nebulizer pressure 30 psi, sheath gas temperature 350 • C, sheath gas flow rate 10 L·min −1 , capillary voltage 3500 V, nozzle voltage 500 V, fragmentor voltage 130 V, skimmer voltage 45 V, Octopole 1 RF Voltage 750 V. Internal calibration was achieved with two calibrants, purine and hexakis (1H,1H,3Htetrafluoropropoxy) phosphazene (m/z 121.0509 and m/z 922.0098), providing a high mass accuracy better than 10 ppm. The data-dependent MS/MS events were acquired for the five most intense ions detected by full-scan MS, from the 200-1000 m/z range, above an absolute threshold of 1000 counts. Selected precursor ions were fragmented at a fixed collision energy of 30 eV and with an isolation window of 1.3 amu. The mass range of the precursor and fragment ions was set as m/z 200-1000.

Data Processing and Analysis
The data files were converted from the .d standard data format (Agilent Technologies, Waldbronn, Germany) to .mzXML format using MSConvert software, part of the ProteoWizard package 3.0 [53]. All mzxml values were processed using MZmine2v51, as previously described [33,34]. Mass detection was realized with an MS1 noise level of 1000 and an MS/MS noise level of 0. The ADAP chromatogram builder was employed with a minimum group size of scans of 3, a group intensity threshold of 1000, a minimum highest intensity of 1000, and m/z tolerance of 0.008 (or 20 ppm). Deconvolution was performed with the ADAP wavelet algorithm according to the following settings: S/N threshold = 10, minimum feature height = 1000, coefficient/area threshold = 10, peak duration range 0.01-1.5 min, and t R wavelet range 0.00-0.04 min. MS/MS scans were paired using a m/z tolerance range of 0.05 Da and t R tolerance range of 0.5 min. Isotopologues were grouped using the isotopic peak grouper algorithm with a m/z tolerance of 0.008 (or 20 ppm) and a t R tolerance of 0.2 min. Peaks were filtered using a feature list row filter, keeping only peaks with MS/MS scans (GNPS). Adduct identification, i.e., sodium-or potassium-cationized species, was performed on the peak list with a retention time tolerance of 0.1 min, a m/z tolerance of 0.008 or 20 ppm, and a maximum relative peak height of 150%. A complex search, such as dimers, was performed with a retention time tolerance of 0.1 min, a m/z tolerance of 0.008 or 20 ppm, and a maximum relative peak height of 150%. Peak alignment was performed using the join aligner with a m/z tolerance of 0.008 (or 20 ppm), a weight for m/z at 20, a retention time tolerance of 0.2 min, and weight for t R at 50. The MGF file and the metadata were generated using the export/submit to GNPS option (https://doi.org/10.5281/zenodo.5733173, accessed on 10 January 2019).
Molecular networks were calculated and visualized using MetGem 1.34 software (https://metgem.github.io/, accessed on 1 March 2019) [34], and MS/MS spectra were window-filtered by choosing only the top 6 peaks in the ±50 Da window throughout the spectrum. The data were filtered by removing all peaks in the ±17 Da range around the precursor m/z. The m/z tolerance windows used to find the matching peaks were set to 0.02 Da, and cosine scores were kept in consideration for spectra sharing at least 2 matching peaks. The number of iterations, perplexity, learning rate, and early exaggeration parameters were set to 5000, 25, 200, and 12 for the t-SNE view.
3.8. Genomic Analyses 3.8.1. gDNA Extraction of P. boydii Strains SNB-CN71, -CN73, -CN81, and -CN85 For high molecular size gDNA extraction, the strains were cultivated on PDA medium on permeable cellophane membranes that were prepared as described by Fauchery et al. [54]. The cellophane membranes were trimmed to the size of the Petri dish, put into boiling distilled water containing EDTA (1 g·L −1 ) 20 min in order to permeabilize the membrane, rinsed four times in a big container with distilled water, and autoclaved. After 7 days at 26 • C microorganisms were removed from cellophane membrane, snap-frozen, and milled using mortar and pestle. gDNA was extracted from microorganisms using NucleoBond Buffer Set III and AXG 100 (Macherey-Nagel, Hoerdt, France) with slight modifications: 250 mg of milled mycelium was resuspended in 5 mL of buffer G3, 37 • C for 14 h; 1.2 mL of buffer G4 was added; and the solution was gently mixed then incubated at 50 • C for 4 h. Sample was clarified through centrifugation 5000× g, 5 min. After ethanol precipitation, samples were re-dissolved in 100 µL of mQ water. Purity control was carried out on a Nanodrop 2000 spectrophotometer, DNA quantity using a Qubit dsDNA BR assay kit according to the manufacture's recommendation and a Qubit fluorometer, and quality control on 0.7% agarose gel.
3.9. Biological Tests 3.9.1. Antimicrobial Assays Pure isolated compounds were tested on human pathogenic microorganisms such as Candida albicans (ATCC 10213) to evaluate their antifungal activities. The test was performed in conformance with a reference protocol from the European Committee on Antimicrobial Susceptibility Testing (EUCAST). The MIC value was obtained after 48 h for C. albicans, and fluconazole (for fungi) was used as a positive control [58].

Antiparasitic Activities Assay for In Vitro Inhibition of T. brucei Brucei Growth
Bloodstream forms of Trypanosoma brucei brucei strain 93 were cultured in HMI9 medium supplemented with 10% FCS at 37 • C under an atmosphere of 5% CO 2 . In all experiments, log-phage cell cultures were harvested by centrifugation at 3000× g and immediately used. Drug assays were based on the conversion of a redox-sensitive dye (resazurin) to a fluorescent product by viable cells and were performed according to the manufacturer recommendations (AlamarBlue ® Assay, Invitrogen Corporation, Carlsbad, CA, USA) [59,60]. Drug stock solutions were prepared in pure DMSO. T. brucei brucei bloodstream forms (3 × 10 4 cells/mL) were cultured as described above in 96-well plates (200 µL per well) either in the absence or in the presence of different concentrations of inhibitors, with a final DMSO concentration that did not exceed 1%. After a 72 h incubation, resazurin solution was added to each well at with a final concentration of 45 µM. Fluorescence was measured at 530 nm excitation and 590 nm emission wavelengths after a further 4 h incubation. The percentage of inhibition of the parasite growth rate was calculated by comparing the fluorescence of parasites maintained in the presence of drug with that in the absence of drug. DMSO was used as control. Pentamidine was used as an anti-trypanosomal drug control (IC 50 = 0.011 ± 0.0017 µM). IC 50 s were determined from the dose-response curves with drug concentrations ranging from 100 µM to 50 nM. The IC 50 value is the mean ± the standard deviation of three independent experiments. Pentamidine was used as an anti-trypanosomal drug control (IC 50 = 0.011 ± 0.0017 µM).

In Vitro Growth Inhibition of Chloroquine-Resistant P. falciparum
The FcB1/Colombian strain of Plasmodium falciparum was maintained in vitro in human erythrocytes in an RPMI 1640 reaction mixture supplemented with 8% (v/v) heatinactivated human serum at 37 • C under an atmosphere of 3% CO 2 , 6% O 2 , and 91% N 2 . The chloroquine-resistant strain FcB1/Colombia of Plasmodium falciparum was obtained from the National Museum Natural History collection, Paris, France (n • MNHN-CEU-224-PfFCB1). Parasites were maintained in vitro on human A+ erythrocytes in RPMI 1640 medium supplemented by 8% (v/v) heat-inactivated human serum at 37 • C under an atmosphere of 3% CO 2 , 6% O 2 , and 91% N 2 . Human red blood cells and serum were provided by the Etablissement Français du Sang, Ile de France (EFS, Paris, France) under the C-CPSL-UNT approval n • 13/EFS/126. In vitro drug susceptibility assays were measured by [ 3 H]-hypoxanthine incorporation as described previously [61]. Drug solutions were prepared in DMSO at a 10 mM concentration. Compounds were serially diluted 2-fold with 100 µL culture medium in 96-well plates. Asynchronous parasite cultures (100 µL, 1% parasitemia and 1% final hematocrit) were then added to each well and incubated for 24 h at 37 • C prior to the addition of 0.5 µCi of [ 3 H]-hypoxanthine (GE Healthcare, Tremblayen-France, France, 1 to 5 Ci·mmol/mL) per well. After a further incubation of 24 h, plates were frozen and thawed. Cell lysates were then collected onto fiberglass filters and counted in a liquid scintillation spectrometer. The growth inhibition for each drug concentration was determined by comparison of the radioactivity incorporated in the treated culture with that in the control culture maintained on the same plate. The concentration causing 50% growth inhibition (IC 50 ) was obtained from the drug concentration-response curve, and the results are expressed as the mean values ± standard deviations determined from several independent experiments. Chloroquine was used as an antimalarial drug control (IC 50 = 0.072 ± 0.0074 µM).

Cytotoxic Assay
Human umbilical vein endothelial cells (HUVECs) were obtained from Clonetics (Lonza; Walkersville, MD, USA) and cultured according to the supplier's instructions. Briefly, HUVECs from three to six passages were subcultured to confluence onto 0.2% gelatin-coated tissue culture flasks in endothelial cell growth medium (EGM2) containing growth factors and 2% FCS. All cell lines were maintained at 37 • C in a humidified atmosphere containing 5% CO 2 . Cell viability was assessed using Promega (Promega France, Charbonnière-les-Bains, France) CellTiter-BlueTM reagent according to the manufacturer's instructions. Cells were seeded in 96-well plates (5 × 103 cells per well) containing 50 mL growth medium. After 24 h of culture, the cells were supplemented with 50 mL of test compound dissolved in DMSO (<0.1% in each preparation). After incubation for 72 h, 20 mL resazurin was added for 2 h before recording fluorescence (λex = 560 nm, λem = 590 nm) using a Victor microtiter plate fluorimeter (PerkinElmer France, Villebon-sur-Yvette, France). IC 50 values correspond to the concentration of test compound that elicited a 50% decrease in fluorescence for drug-treated cells relative to untreated cells. Experiments were performed in triplicate. Fumagillin was used as a positive control and showed an IC 50 of 0.006 ± 0.002 µM.

Conclusions
This study combined genomic and metabolomic approaches to study metabolome produced by mutualistic fungal strains of termites from French Guiana of Pseudallescheria genus. It allowed the annotation of tyroscherin analogs previously described, as well as pseurotin A, but also new ovalicin analogs.
Our P. boydii strain collection is a source of diverse ovalicin derivatives, including new fatty acid esters of ovalicin and boyden C that have been isolated and structurally characterized. Compounds 2-4, similar to ovalicin, exhibited antiplasmodial and antitrypanozomal properties. Our work also confirms that P. boydii strains are associated with termite biosynthesized compounds with biological activities. We were able to demonstrate that antiparasitic compounds such as Compounds 1 to 4 are produced in different amounts by the different strains of P. boydii. These same compounds reduce cytotoxicities, with IC 50 s greater than 10 µM on HUVECs.
The molecular networking analysis allowed the annotation of all isolated compounds and many other ovalicin analogues, together with the isolation and identification of the new ovalicin analog boyden C. Finally, it was possible by genome mining to propose a common biosynthetic pathway for the ovalicin analogs compounds produced by Pseudallescheria sp. strains of our collection. The biosynthesis also permitted the deduction of absolute configurations of ovalicin analogs.
Fungal natural products possess biological activities that are of great value for medicine, agriculture, and manufacturing. Recent "omics" studies accentuated the immense taxonomic diversity of fungi, and the accompanying specialized metabolic diversity provides a significant and still largely untapped resource for natural product discovery. Particular attention will be given in the future to integrative pan-genomic approaches that use a combination of genomic and metabolomic data for parallelized natural product discovery across multiple strains. Such novel technologies will not only expedite the natural product discovery process but will also allow the assembly of a high-quality toolbox for the re-design or even de novo design of biosynthetic pathways using synthetic biology approaches, while continuing to focus on specific environments such as symbiotic microorganisms.