To investigate the role of Nf1 in the mammary tumorigenesis, we previously generated three mutant lines representing distinct states of Nf1 deficiency (Figure 1A). (Dischinger et al., 2018) In this model, we targeted exons 21–22 (previously denoted as exons 20–21) within the CSRD (cysteine-serine rich domain) of Nf1. The CSRD is an allosteric activator of the adjacent GRD (GAP-related domain) and near the dimerization region of neurofibromin (Lupton et al., 2021; Young et al., 2023). Recent studies have revealed that hotspot mutations in the CSRD region lead to severe phenotypes in NF1 individuals (Koczkowska et al., 2018; Koczkowska et al., 2020). The IF/+ animals harbor an in-frame indel in exon 21, whereas the PS/+ rats have a premature stop codon in exon 21. The IF; PS-21/+ line was derived from an in-frame indel in exon 21 (previously exon 20) and a distinct mRNA isoform with a premature stop in exon 22 (previously exon 21). By creating indels in the CSRD of one copy of the Nf1 gene, each of these Nf1-deficient immunocompetent rat lines developed highly penetrant, aggressive mammary adenocarcinomas with varying times on tumor onset. The IF/+ line, with the least impact on neurofibromin function, had the slowest tumor onset, the PS-20/+ and IF; PS-21/+ lines with a premature stop indel (thus removing all the domains after the CSRD) resulted in animals with quicker tumor onset (Figure 1B). The difference in tumor initiation between the Nf1 indels offers distinct models in which to study the pretumorigenic changes, with a long window of before tumor onset in IF/+ rats compared to the rapid tumor onset that occurs in the PS-20/+ and IF; PS-21/+ lines.

To understand the mechanism of different tumor onset, it is vital to understand the changes that occur in preneoplastic glands that create an environment of tumor permissiveness. While we have previously shown an increase in mammary epithelial branching and end bud formation during development (Tovar et al., 2022), H&E staining also shows an increased matrix region around the mammary ducts and end buds of Nf1 deficient rats compared to WT. To quantify this difference, we measured the percent area of eosin around the ducts. At 21 days, the odds of having eosin staining are 14 times higher in Nf1-deficient IF; PS-21/+ tissue than WT tissues (p = 0.005, Figure 1C). However, at 33 days, the odds of having eosin staining in Nf1-deficient PS-20/+ and IF; PS-21/+ tissues are 5 times (p = 0.002) and 8 times (p < 0.001) higher respectively compared to WT tissues (Figure 1D). At both 21 and 33 days, the odds of having more eosin staining in Nf1-deficient IF/+ tissues than WT tissues are not statistically significant with the odds of having more eosin staining only 2 times higher odds (p = 0.619 at 21 days; p = 0.959 at 33 days).

To verify that this increase in eosin correlates with an increase in collagen deposition that progresses with age towards tumor formation, we stained tissues from 43-day-old rats with picrosirius red (PSR) that allows for collagen visualization under exposure of polarized light (Coelho et al., 2018). At 43 days, all Nf1-deficient mammary glands had increased collagen content compared to the WT glands. Consistent with the eosin quantification, both Nf1-deficient PS-20/+ and IF; PS-21/+ tissues showed approximately 150% increase in polarized light (POL) signals compared to the WT tissues (p < 0.01; Figure 1E). Interestingly, although still lower than the other two Nf1-deficient genotypes, the IF/+ tissues at 43 days also had an increase in POL signals as well, about 72% more than the WT tissues (p = 0.02; Figure 1E). This increase in collagen indicates that as the IF/+ rats approach an age closer to tumor formation, their collagen content also increases.

This indicates that collagen and matrix changes correspond to Nf1 loss of function and rates of tumor onset. In our model, increases in collagen deposition are present earlier in PS-20/+ and IF; PS-21/+ mammary glands which have a larger loss of neurofibromin and earlier tumor onset (6–12 weeks). In contrast, increase in collagen changes are first observed at 43 days, in IF/+ mammary glands that have a minor alteration in neurofibromin and a much later tumor onset (<10 months). These results show that the timing of increase in collagen mirrors the timing of tumor onset in the Nf1-deficient lines.

Additionally, to delineate whether the collagen deposition increased solely as a result of epithelial ductal development, we measured fibrosis in mammary glands of older IF/+ rats that did not develop palpable tumors at 7.5 months old. Interestingly, we discovered a noticeable increase in fibrotic region around epithelial ducts that appear to have ductal carcinoma in situ (DCIS; Supplementary Figure S1). We are unable to measure collagen content at later timepoints in PS-20/+ and IF; PS-21/+ lines because they develop mammary adenocarcinomas with complete penetrance by 14 weeks of age. However, in the slower progressing IF/+ line, increased collagen content does persist past the epithelial gland developmental stage into adulthood. This finding underscores the impact of Nf1 loss of function in remodeling the rat mammary stroma during both pubertal development and adulthood.

To identify which cell population and what type of collagen was predominantly involved in the elevated collagen levels, we performed RNAseq of adipocytes, fibroblasts and epithelial cells isolated from mammary glands before tumor formation. Unsupervised clustering of collagen gene expression revealed that each sample expresses different types of collagens with very little overlap between each other (Supplementary Figure S2). As expected, fibroblasts have the highest expression of fibrillar collagens, such as collagen I, III, and V, yet there were no significant differences between the Nf1-deficient and WT fibroblasts. Epithelial cells and mature adipocytes primarily express distinct subunits of network forming collagens. Surprisingly, we observed that the Nf1-deficient IF/+ mature adipocytes have increased gene expression of collagens (Figure 2A). Because adipocytes are more well known for their endocrine and metabolic functions, we were intrigued by these unexpected findings.

Recent studies have revealed that adipocytes surrounding a tumor can become activated to form cancer-associated adipocytes (CAAs) and have increased collagen secretion (Dirat et al., 2011; Zhang et al., 2012; Bochet et al., 2013; Wei et al., 2019; Zhu et al., 2022). Since the adipocytes we analyzed are from mammary glands before tumor formation, these results indicate that Nf1 plays a role in normal mammary adipocyte differentiation. Next, we examined whether the IF/+ mature adipocytes also exhibit other reported gene expression features of CAAs. Along with increased collagen secretion, these studies also reported that CAAs can remodel collagen architecture and alignment through increased expression of collagen modification genes (Wei et al., 2019). An example of a collagen modification gene is prolyl 4-hydroxylase subunit alpha 1 (P4ha1), which is the most important catalytic subunit of the P4H enzyme and necessary for collagen polypeptide chains folding into stable triple-helical molecules (Myllyharju and Kivirikko, 2004). From our RNAseq data, P4ha1 and other collagen modification genes have increased expression in IF/+ mature adipocytes compared to the WT and two other Nf1-deficient mature adipocytes (Figure 2B). Additionally, gene ontology (GO) analysis of Nf1-deficient IF/+ mature adipocytes gene expression compared to WT showed enrichment in pathways involved in ECM processes and RAS-MAPK signaling (Figure 2C). Pathways involved in MAPK signaling were expected due to neurofibromin’s role as a negative regulator of RAS. Enrichment of pathways in ECM processes substantiates a role for Nf1-deficient IF/+ adipocytes in collagen deposition and a CAA-like phenotypes (Figures 2A, B). Our analysis indicated other ECM components besides collagen may contribute to tumorigenesis in Nf1-deficient IF/+ glands. Several genes typically known for their role in tissue fibrosis (i.e., tenascin, fibrinogen, and fibronectin genes) were increased in Nf1-deficient IF/+ adipocytes (Figure 2D). Overall, these results demonstrated that Nf1-deficient IF/+ adipocytes have increased collagen and fibrosis-related ECM gene expression preceding tumor formation.

Another reported phenotype of CAAs is their ability to dedifferentiate to a preadipocyte- and fibroblast-like cell through a process called adipocyte mesenchymal transition (AMT) (Dirat et al., 2011; Zhang et al., 2012; Bochet et al., 2013; Zhu et al., 2022). To determine if AMT was present in the Nf1-deficient IF/+ mature adipocytes, we evaluated adipocyte differentiation expression signatures. Principal component analysis (PCA) of adipocyte and fibroblast samples show that the Nf1-deficient IF/+ mature adipocytes cluster further away from the other Nf1-deficient genotypes, indicating that their gene expression profile has deviated from a typical mature adipocyte (Figure 3B). Looking at the different genes expressed at adipocyte, preadipocyte, and fibroblast stages, our RNAseq data showed the Nf1-deficient IF/+ mature adipocytes had lower expression of mature adipocyte markers, but higher expression of fibroblast and preadipocyte markers. WT and Nf1-deficient PS-20/+ and IF; PS-21/+ mature adipocytes do not have altered mesenchymal lineage gene expression. For example, adiponectin (Adipoq) and perilipin-1 (Plin1), genes which are primarily expressed in adipose tissues and mature adipocytes but are both downregulated in the IF/+ mature adipocytes. In contrast, fibroblast activation protein alpha (Flap) and Cd34 are genes primarily expressed in activated fibroblast and preadipocytes, respectively, but have increased expression in the IF/+ mature adipocytes compared to the other genotypes (Figure 3A). Gene set enrichment analysis (GSEA) of the mature adipocytes revealed an enrichment in adipocyte differentiation gene signatures (obtained from the MSigDB database). To highlight one example, the Nf1-deficient mature adipocytes showed an enrichment in a fat cell differentiation gene set (GO0045444) which supports the difference in mesenchymal lineage markers we observed (NES = 1.29, p = 0.08; FDR = 0.024; Figure 3C). Overall, these results show that the IF/+ indel in one Nf1 allele causes the gene expression profile of mature adipocytes before tumor formation to mimic a preadipocyte- and fibroblast-like cell state.

To verify the altered differentiation gene expression observed in our RNAseq data (Figure 3), we isolated the adipose-derived stem cell (ASC) population from 33-day-old rat mammary glands and differentiated them in vitro as illustrated in Figure 4A. After an overnight collagenase digest, we cultured the stromal vascular fraction (SVF) for 2-3 passages to enrich for the ASC population before culturing them in commercially available differentiation media for 2 weeks. After 2 weeks, fully differentiated cells containing lipid droplets were counted as a measure of differentiation ability. After 14 days in differentiation media, we observed that Nf1-deficient IF/+ ASCs have 60% fewer cells with lipid droplets compared to WT, thus indicating a decreased differentiation ability (WT vs. IF/+ odds ratio = 0.4, p = 0.08; Figure 4C).

Mature adipocytes and preadipocytes differentiate through a multistep process originating from mesenchymal stem cells (MSCs), and in adipose tissues specifically, these tissue resident MSCs are known as ASCs (Berry and Rodeheffer, 2013; Berry et al., 2015; Brown et al., 2022). To investigate whether the reduced differentiation ability of the Nf1-deficient IF/+ ASC that we observed (Figure 4) was due to a decreased ASC population, we used flow cytometry to quantify the amount of ASC present in the SVF of our rat mammary glands (Figure 5A). After red blood cell removal, we derived a lineage-negative population by eliminating CD45⁺ hematopoietic cells, CD31⁺ endothelial cells, and E-Cad + epithelial cells. Using three stem cell markers, CD34⁺ CD29⁺ CD90⁺, we further stratified the ASC population, before excluding CD54⁺ fibroblasts (Rodeheffer et al., 2008; Berry and Rodeheffer, 2013; Berry et al., 2015; Jeffery et al., 2015; Jeffery et al., 2016). The abundance of the resulting ASC population (Lin- CD34⁺ CD29⁺ CD90⁺ CD54⁻) was measured compared to the whole SVF population. Even though substantial heterogeneity of ASC was observed within the genotypes, we did not detect a decrease of the Lin- CD34⁺ CD29⁺ CD90⁺ CD54⁻ ASC population in Nf1-deficient IF/+ rats compared to WT (Figures 5B, C). Together, our differentiation assays and flow cytometry analysis indicate that Nf1-deficient IF/+ adipocytes have impaired differentiation that is not due to a decrease in the adipocyte progenitor population. In addition to demonstrating that Nf1 does not impact the early stages of differentiation, these data suggest a potential of Nf1 regulating the later stages of adipocyte differentiation instead, specifically through PPARγ signaling.

Female and male immunocompetent Nf1-deficient Sprague Dawley rats were bred in-house at the Van Andel Research Institute. The four genotypes include WT, and three Nf1-deficient lines: IF/+, PS-20/+, and IF; PS-21/+. Male Nf1-deficient rats were bred to wild-type female rats and female pups were used at the indicated developmental time-point. Animals were housed two to a cage, had ad libitum access to food and water, and were kept under a 12 h light - 12 h dark cycle. All animal protocols were approved by the Van Andel Research Institute Animal Care and Use Committee.

Mammary tissues were formalin fixed and paraffin embedded utilizing routine histological procedures. Sections for H&E staining were cut at 5 μm thickness and was performed by Van Andel Institute Pathology and Biorepository Core using a Leica Rotary microtome. Staining was performed on the Tissue-Tek Prisma Plus autostainer, using the Sakura H&E Kit #1. Slides were scanned at ×20 magnification on either a Leica Aperio XT or a Leica Aperio AT2 scanning system by the Van Andel Institute Pathology and Biorepository Core.

The H&E collagen area analysis was conducted by Van Andel Institute Optical Imaging Core. Full resolution svs images of H&E-stained mammary tissue sections were imported into QuPath (version 0.4.2) for analysis (Bankhead et al., 2017). QuPath’s Estimate Stain Vector tool was used to automatically deconvolve the H&E stain color vectors for each section. A region of interest (ROI) contouring the mammary gland and excluding large blood vessels and lymph nodes was manually generated for each tissue section. A pixel classifier was then created using QuPath’s Pixel Classifier Tool to measure the collagen within the mammary gland ROI. The pixel classifier used an eosin stain-based threshold of 0.05, pixel resolution of 1.0062 um (equivalent to a two factor down sample), and Gaussian blur of 1.5; positive pixels were converted into an area measurement within QuPath. The areas of the mammary gland contours and pixel classifier measurements were exported into a csv file for further statistical analysis.

Picrosirius red (PSR) staining was conducted by Van Andel Institute Pathology and Biorepository Core using a standard protocol (Dolber and Spach, 1987). Briefly, tissue sections 5-μm-thick tissue sections were deparaffinized, treated with phosphomolybdic acid before staining with Sirius red. Slides were washed in hydrochloric acid and 70% EtOH before cover slipping for polarized light imaging.

Polarized light images were acquired by Van Andel Institute Optical Imaging Core on a Zeiss AxioScan 7, using an EC Plan-Neofluar 20x/0.5 N.A. Pol M27 air objective, and ZEN Blue software (version 3.7). Tissues were illuminated with a Colibri 7 TL LED light source, imaged in brightfield at 5% intensity and 200 μs exposure, and by circularly polarized light at 10% intensity and 3.78 m exposure, sequentially. Samples were detected by an AxioCam705 color CMOS camera with 5-megapixel resolution. Resulting images were 24 bits and at 0.173 × 0.173 um scaling per pixel.

Full resolution czi images of polarized light collected from picrosirius red-stained mammary tissue sections were imported into QuPath for analysis by Van Andel Institute Optical Imaging Core. A ROI contouring the full tissue section was manually drawn for each slide. A pixel classifier was created using QuPath’s Pixel Classifier Tool to measure the area covered by the signal from polarized light within the ROI. A training image was created from 12 equal sized rectangular regions from three different sections using QuPath’s create training image tool. This resulting training image was annotated for examples of signals from the polarized light (n = 53) and background (n = 33) and used for training the pixel classifier using an artificial neural network at 1.38 μm/pixel resolution. All other settings in the pixel classifier were left at their default options. The trained pixel classifier was applied to the ROIs of all images to segment the total area covered by the polarized light signal. The segmented areas were exported from QuPath into a csv file for further statistical analysis.

Freshly isolated fourth mammary fat pads (MFP) were digested into a single cell suspension as described in Dischinger et al. (2018) and Tovar et al. (2020). Briefly, the fourth mammary fat pads from female rats are removed and digested overnight. Mature adipocytes are present in the top lipid layer and are carefully collected for RNA extraction. Red blood cells are depleted by centrifugation and lysis, respectively. The remaining cells are negatively selected using Miltenyi CD31 and CD45 biotin conjugated antibodies. Next, addition of Miltenyi anti-biotin microbeads allows capture of the antibody tagged lineage positive populations on an LS column with a Miltenyi QuadroMACS magnet. The eluted mammary epithelial cells (MECs) and fibroblasts are cultured in ultra-low attachment (ULA) plates for 48 h to deplete fibroblasts and isolate only MECs. Alternatively, to deplete MECs and isolate only fibroblast, cell pellets are cultured using differential attachment method or with fibroblast media containing gentamicin.

For RNA extraction of whole fat pads and tumors, fresh harvested samples were cut up and added to MPbio lysing matrix E tube. Cells were homogenized for 20 s in an MPbio homogenizer, then RNA was isolated from the supernatant using the Trizol-chloroform-isopropanol RNA extraction method. For RNA extraction of mature adipocytes, RNA was isolated using the Trizol-chloroform-isopropanol RNA extraction method. For RNA extraction from cultured epithelial cells and fibroblast, 1.2×10⁵ cells were plated into one well of a six-well plate. At 90%–100% confluency, RNA was isolated using QIAGEN QIAshredder (QIAGEN #79654) and RNeasy Kit (QIAGEN #74104) manufacturer’s instructions.

Libraries were prepared by the Van Andel Institute Genomics Core from 500 ng of total RNA using the KAPA mRNA Hyperprep kit (v4.17) (Kapa Biosystems, Wilmington, MA United States). RNA was sheared to 300–400 bp. Prior to PCR amplification, cDNA fragments were ligated to IDT for Illumina TruSeq UD Indexed adapters (Illumina Inc., San Diego CA, United States). Quality and quantity of the finished libraries were assessed using a combination of Agilent DNA High Sensitivity chip (Agilent Technologies, Inc.), QuantiFluor^® dsDNA System (Promega Corp., Madison, WI, United States), and Kapa Illumina Library Quantification qPCR assays (Kapa Biosystems). Individually indexed libraries were pooled and 50 bp, paired-end sequencing was performed on an Illumina NovaSeq6000 sequencer using an S2 sequencing kit (Illumina Inc., San Diego, CA, United States) to an average depth of 45M reads per sample. Base calling was done by Illumina RTA3 and output of NCS was demultiplexed and converted to FastQ format with Illumina Bcl2fastq v1.9.0.

Raw reads were first trimmed of adapters with Trim_Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and then mapped to Rnor 6.0 assembly with STAR v2.5.2b using options “--twopassMode Basic” and “--quantMode GeneCounts” to directly output counts for all features from the ensemble annotation 6.0.90 (Dobin et al., 2013). A negative binomial generalized log-linear model was then fit to the filtered count data with edgeR using the weighted trimmed mean of M-values to normalize for library size and composition biases (Robinson et al., 2010). Different genotype or tissue groups were contrasted, and p-values were generated using empirical Bayes quasi-likelihood F-tests, and then adjusted using the BH method; adjusted p-values less than 0.05 were considered significant. Heatmaps were generated from library-size normalized counts centered across genes (z-scores) using the R package ComplexHeatmap (Gu et al., 2016). Using the normalized counts per million, the fold change against matched genotype normal mammary gland expression was calculated for tumor samples and shown as log2 (fold change) capped at −3 and 3.

The gene ontology (GO) functional enrichment analysis was performed using g:Profiler (version e106_e.g.,53_p16_65fcd97) with g:SCS multiple testing correction method applying significance threshold of q < 0.05 (Raudvere et al., 2019). Results were plotted using ggplot2 (v3.3.6) (Wickham, 2016). For GSEA, protein-coding genes were first pre-ranked by fold change of mutant vs. WT and loaded onto GSEA v4.3.2 (Subramanian et al., 2005). Analysis was performed on adipocyte differentiation gene signatures obtained from the MSigDB database (Liberzon et al., 2011; Liberzon et al., 2015).

SVF isolation was adapted from Picon-Ruiz et al. (2020). After mammary fat pad digestion and red blood cell removal as described above, the SVF fraction was washed with PBS, filtered through a 70-μm cell strainer, resuspended, and plated at 3×10⁵ cells/cm². After 2-3 passages in DMEM supplemented with 10% FBS and 1% penicillin/streptomycin, the cells remaining in culture were taken as the ASC population.

Adipogenic differentiation was induced using Human MesenCult™ Adipogenic Differentiation Kit (STEMCELL Technologies #05412). According to manufacturer’s protocol, DMEM media was replaced with differentiation media once the ASC culture reached 90%–100% confluency. ASC culture was differentiated for 14 days with media changes every 3 days before they were used for downstream applications.

After 14 days in adipogenic differentiation media, cells were stained with Oil Red O following a standard histological protocol (Suvarna et al., 2019). Briefly, cells were fixed with 4% paraformaldehyde for 1 h at room temperature, then incubated with 60% isopropanol for 5 min before staining with Oil Red O (Sigma-Aldrich #O0625) for 30 min. After staining, cells were washed with distilled water, stained with Hematoxylin (Vector Laboratories #H-3401) for 1 min and finally washed with tap water for bluing before visualization/imaging.

For flow cytometry analysis, 1 × 10⁶ cells in 100 μL were stained for 30 min at room temperature using the following antibodies: CD90-BUV395 (BD, #740260, 0.1 mg/μL), ECAD-AF405 (Novus Bio, #NBP2-98433AF405, 1.74 mg/μL), CD45-BV605 (BD, #740371, 0.4 mg/μL), CD29-SB702 (Thermo, #67–0291–82, 0.1 mg/μL), CD34 (Thermo, #PA5-47849, 1 mg/μL) with secondary AF488 (Thermo, #A-11015, 1 mg/μL), CD24-PE (Thermo, #562104, 0.2 mg/μL), CD54-R718 (BD, #751649, 0.1 mg/μL), CD140a-APC (BD, #568234, 0.2 mg/μL), and CD31-PECy7 (Thermo, #25–0310–82, 0.025 mg/μL). All fluorescence minus one (FMO) and single stain controls were stained in HF buffer (HBBS 1X with calcium chloride and magnesium chloride, supplemented with 2% FBS and 1% Pen/Step). Full panel staining was done with HF buffer and 10 μL Brilliant Stain Buffer Plus (BD Horizon #566385) per sample. Finally, cells were washed, stained with 0.005 μg/mL DAPI, filtered, and acquired on the Cytek Biosciences Aurora spectral flow cytometer by the Van Andel Institute Flow Cytometry Core. Data was analyzed using FlowJo v10.8.1. For in vitro phenotypic validation, cells were similarly stained and then sorted on the BD FACSymphony S6. Sorted cells were cultured in DMEM with 10% FBS and 1% Pen/Strep.