Elevated miR-29c-5p Expression in Nipple Aspirate Fluid Is Associated with Extremely High Mammographic Breast Density

Simple Summary High mammographic density is a known risk factor for breast cancer. However, the underlying mechanisms of high mammographic density development and breast cancer are unknown. MicroRNAs are potential biomarkers indicative of carcinogenesis and can be assessed in nipple aspirate fluid. We used nipple aspirate fluid from women with very low and extremely high mammographic density to examine differences in expression of multiple miRNAs between both extremes in the spectrum of mammographic density. We found that hsa-miR-29c-5p was upregulated in an extremely high mammographic density context and potential targets were identified that might provide clues of the relationship between high mammographic density and breast cancer risk. Understanding the relationship between high mammographic density and breast cancer is of great value for early breast cancer diagnosis and treatment. With our research we provide new insight into this relationship and further research could determine the effects of dysregulated hsa-miR-29c-5p on the identified candidate targets. Abstract High mammographic density (MD) is associated with an increased risk of breast cancer, however the underlying mechanisms are largely unknown. This research aimed to identify microRNAs (miRNAs) that play a role in the development of extremely dense breast tissue. In the discovery phase, 754 human mature miRNAs were profiled in 21 extremely high MD- and 20 very low MD-derived nipple aspirate fluid (NAF) samples from healthy women. In the validation phase, candidate miRNAs were assessed in a cohort of 89 extremely high MD and 81 very low MD NAF samples from healthy women. Independent predictors of either extremely high MD or miRNA expression were identified by logistic regression and linear regression analysis, respectively. mRNA targets and pathways were identified through miRTarBase, TargetScan, and PANTHER pathway analysis. Statistical analysis identified four differentially expressed miRNAs during the discovery phase. During the validation, linear regression (p = 0.029; fold change = 2.10) and logistic regression (p = 0.048; odds ratio = 1.38) showed that hsa-miR-29c-5p was upregulated in extremely high MD-derived NAF. Identified candidate mRNA targets of hsa-miR-29c-5p are CFLAR, DNMT3A, and PTEN. Further validation and exploration of targets and downstream pathways of has-miR-29c-5p will provide better insight into the processes involved in the development of high MD and in the associated increased risk of breast cancer.


Introduction
Mammography is currently the primary screening method for timely detection and treatment of breast cancer. The extent of fibro-glandular versus fatty tissue on a mammogram can classify breast tissue into four categories of density, from A to D. A higher A total of 110 healthy women with extremely high MD breasts (classified as an ACR4 or 'D' on a mammogram using Volpara software (Volpara Solutions)) and 101 healthy women with very low MD breasts (ACR1 or 'A'), representing the two extreme categories within the four-category MD spectrum (A to D [1]), were included in this study after informed consent was obtained. Sample size was determined by available samples. These women were participants of the larger, nationwide DENSE trial (NCT01315015; [15,16]) and the DENSE-on low biobanking study, approved by the Institutional Review Boards within the participating hospitals, and the UMC Utrecht Biobank Research Ethics Committee (TCBio; biobank study number 14-467 approved on 17 June 2015). All women were aged between 50-74 years and had a negative screening mammography and/or an MRI without abnormalities at the time of NAF acquisition. Median age was 55 years old, with an interquartile range (IQR) of 52-58. 5. In the group of women with entirely fatty breasts, age was restricted to 50-60 years to prevent the effect of age on mammographic density. Information on lifestyle, reproductive characteristics and medical history were obtained by means of a questionnaire.
Bilateral NAF samples were acquired between June 2015 and March 2016 at several Dutch hospitals or at the home of the study subjects by trained research nurses. Collection was performed using a vacuum device after nasal oxytocin administration, as described earlier by Suijkerbuijk et al. [14]. The collected fluid was conserved in a buffer solution (50 mM Tris pH 8.0, 150 mM NaCl, 2 mM EDTA) at −80 • C until required for analysis. Details concerning sample collection success, sample volume (between 5-100 µL) and sample color were registered. Analyses were conducted in two main phases, a discovery phase (21 extremely high and 20 very low MD NAF samples) and a validation phase (89 extremely high and 81 very low MD NAF samples). The Biospecimen Reporting for Improved Study Quality (BRISQ) guidelines were taken into account [17].

RNA Isolation, Reverse Transcription, and Pre-Amplification
Total RNA was extracted from pooled NAF samples (intra-individual samples from left and right breast combined, for discovery) or unilateral NAF samples (left or right breast at random, for validation) according to the manufacturer's protocol using the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen, Hilden, Germany). All isolations were performed starting with 20 µL of pooled NAF or 10 µL of unilateral NAF if available. Pooled NAF was used in the discovery phase to ensure sufficient RNA for profiling. Non-human synthetic ath-miR159a (with a 5 phosphate) was spiked in as procedural control at 300 pg by pre-mixing with RLT plus lysis buffer. Total RNA was eluted in 30 µL RNAse-free water. RNA concentrations were determined with the Qubit RNA HS Assay Kit (Invitrogen, Q32852, Waltham, MA, USA) measured by Qubit 3.0 (ThermoFisher Scientific, Waltham, MA, USA) fluorometric quantification.
A uniform RNA sample concentration was obtained by dilution in RNAse-free water to 4 ng/µL or to 2.5 ng/µL for the discovery phase and the validation phase, respectively. First, according to the manufacturer's instructions, 8 ng (discovery phase) or 5 ng (validation phase) of total RNA was poly-A tailed. After adaptor ligation and reverse transcription, cDNA was pre-amplified for nineteen cycles using the TaqMan Advanced miRNA cDNA Synthesis Kit (ThermoFisher Scientific, Waltham, MA, USA) on a Veriti 96-well thermal cycler (ThermoFisher Scientific, Waltham, MA, USA). The pre-amplification product was subsequently diluted 20× in 0.1× Tris buffer, pH 8.0 and stored at −20 • C until qPCR.

Discovery Phase: Taqman OpenArray Profiling Analysis of Nipple Aspirate Fluid
The expression levels of 754 human mature miRNAs (Supplementary Table S1), that were functionally validated with miRNA artificial templates, were profiled using the fixedcontent TaqMan OpenArray Human Advanced MicroRNA Panel on a QuantStudio 12K Flex system (ThermoFisher Scientific, Waltham, MA, USA). The samples were loaded from a 384-well plate onto the TaqMan OpenArray Human Advanced MicroRNA Panel array slide using the OpenArray AccuFill system. Relative threshold values (CRT [18]) were automatically generated using the ThermoFisher Cloud system (https://apps.thermofisher. com/ (accessed on 24 November 2020)). These were proven to be more robust than baseline threshold values for analysing data generated using nanolitre fluidics based OpenArray plates. Analysis settings included the following restrictions: a minimum CRT of 10, a minimum AMPSCORE (low signal in linear phase) of 1, a minimum calculated confidence in the quantification cycle (CQCONF) (Cq) value of 0.6, a maximum CRT of 30 with inclusion of maximum CRT in calculations. Additionally, all miRNA amplification plots were visually inspected on curve shape and signal timing. The global mean was used for normalization [19,20].

Validation Phase: Individual TaqMan Advanced miRNA qPCR Assays
Validation of differentially expressed candidate miRNAs identified in the discovery phase was performed in a larger, independent NAF cohort (n = 170, 89 extremely high and 81 very low MD samples). Individual TaqMan advanced mature miRNA assays for hsa-miR-92a-3p, hsa-miR-92b-3p, hsa-miR-22-5p, hsa-miR-29c-5p, hsa-miR-125a-5p (endogenous control identified by GeNorm analysis (https://genorm.cmgg.be/ (accessed on 24 November 2020)) [20]), and ath-miR159a (technical control) (for assay IDs see Supplementary Table S2) were used with TaqMan Fast Advanced Master Mix (ThermoFisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions on a ViiA7 real-time PCR system (ThermoFisher Scientific, Waltham, MA, USA). All qPCR reactions were performed in duplicate and interplate calibrator samples and non-template controls were used throughout. Cycle threshold (CT) values >35 were considered undetermined. The data were analyzed using comparative quantification [22]. Delta CT values (CT target miRNA minus CT endogenous control miRNA) were calculated and used for further statistical analysis with hsa-miR-125a-5p as endogenous control. Supplementary Table S3 shows raw CT values per sample. Analyses were performed after replacement of missing miRNA values (missings in <3/4 miRNAs per sample) by the maximum CT value obtained for the particular miRNA +1 [23].

Statistical Analysis
Established risk factors for breast density and breast cancer were included in the statistical analysis. Body mass index (BMI) values, calculated based on self-reported height and weight, were derived from the moment of inclusion in the DENSE(-on) study. The baseline characteristics and reproductive information of the study subjects were analyzed as continuous data (for age, BMI, age at menarche and age at first live birth) or as dichotomous categorical data (for parity (nulliparous or parous) or having a first-degree family member with breast cancer (yes/no)).
Statistical analyses were performed using IBM SPSS Statistics for Windows version 27.0 (IBM Corp., Orchard Road Armonk, NY, USA) and RStudio version 1.3.1093 (Public-benefit corporation, Vienna, Austria). Normal distribution of the data was analyzed by visual inspection (histograms, residual, and Q-Q plots) and by means of a Kolmogorov-Smirnov test. Differences between patient characteristics were assessed by ANOVA and independent sample t-tests for normally distributed data, whereas Mann-Whitney U-tests were used for non-normally distributed data. For dichotomous data comparison, Pearson chi-square or Fisher's exact tests were applied. A Pearson's correlation test was performed to compare profiling and individual assay results within the discovery set, and to explore collinearity.
To identify differentially expressed miRNAs between extremely high and very low MD cohorts within the discovery set, a two-step approach was used. First, using the delta CT (DCT) values of 118 miRNAs expressed in >80% of discovery samples (>32/41), a linear regression was performed in SPSS with miRNA DCT value as outcome variable and BMI, dense category (extremely high or very low), parity (parous or nulliparous), age, and NAF color class as predictor variables (enter method, with variable entry at p < 0.05 and removal at p < 0.10). The latter (NAF color class) was included in the linear regression model because previous in-house data showed a significant association between these parameters and miRNA expression [24]. NAF color was classified into four categories (clear white/yellow, turbid white/yellow, bloody/orange/pink and green/brown), of which dummy variables were made. miRNAs with p < 0.20 for breast density as predictor were noted. Second, using the same 118 miRNAs, a Firth's bias-reduced logistic regression was performed in RStudio (package logistf, version 1.24) using dense category as outcome, and the other variables (BMI, age, and NAF color) including miRNA DCT values as predictors. In this analysis, the NAF color variable was reduced to three categories: white/yellow, bloody/orange/pink and green/brown. Again, miRNAs with p < 0.20 for predicting dense category and their coefficients were noted. All miRNAs overlapping between both analyses were selected for subsequent validation experiments. To identify differentially expressed miRNAs between extremely high and very low MD cohorts within the validation set, we used the same approach. Fold change in extremely high versus very low MD and odds ratios with 95% confidence intervals (CI) were calculated for, respectively, the linear regression and the Firth's bias-reduced logistic regression approach [22]. Here, miRNAs with p < 0.05 were considered significant. GraphPad Prism 8.3 for Windows (San Diego, CA, USA) was used for graphical visualization of the results.

Study Subject and NAF Characteristics Per Cohort
Baseline characteristics of extremely high and very low MD discovery and validation cohorts are shown in Table 1. Between the two extreme MD categories, significant differences were observed for the variables body mass index (BMI) (discovery and validation cohorts showed lower BMI in extremely high MD, both p < 0.0001), age at menarche (validation cohort older menarche age in extremely high MD cohort, p = 0.001) and NAF color (discovery cohort more bloody/orange/pink NAF samples in extremely high MD, p = 0.028). Table 1. Baseline characteristics of extremely high ("high") and very low ("low") mammographic density cohorts. (a) Discovery cohort and (b) validation cohort. Cohort size per baseline characteristic can differ due to missing values. Bold p-values indicate significant difference.  No significant differences were observed between the discovery and validation cohorts concerning age at the time of NAF collection, BMI, age at first live birth, age at menarche, parity, and having a first-degree family member with breast cancer (Supplementary  Table S4). NAF color was however significantly different (p < 0.0001) between discovery and validation cohorts. Total RNA concentrations in both cohorts varied between 2.5 and 134 ng/µL, with a median RNA concentration of 33 ng/µL in the discovery cohort and 5 ng/µL in the validation cohort (p < 0.0001).
The reliability of the microfluidics-based discovery was further technically validated with regular-volume qPCR, showing a fair to high correlation for all 14 randomly selected miRNAs (Pearson correlation coefficient range 0.37-0.93) (Supplementary Figure S3).

Validation: hsa-miR-29c-5p Is Differentially Expressed between Extremely High and Very Low MD
Eleven NAF samples (6.5% of the validation cohort) showed undetermined CT values for all four candidate miRNAs. Another three samples (1.8%) showed undetermined or late CT (>31) for at least three of the interrogated miRNAs, and these were excluded Unsupervised hierarchical clustering of NAF samples based on the expression level of the four candidate miRNAs (without correction for confounders) revealed two clusters, one cluster containing relatively more extremely high-density samples and one cluster containing relatively more very low-density samples (Supplementary Figure S2).
The reliability of the microfluidics-based discovery was further technically validated with regular-volume qPCR, showing a fair to high correlation for all 14 randomly selected miRNAs (Pearson correlation coefficient range 0.37-0.93) (Supplementary Figure S3).

Validation: Hsa-miR-29c-5p Is Differentially Expressed between Extremely High and Very Low MD
Eleven NAF samples (6.5% of the validation cohort) showed undetermined CT values for all four candidate miRNAs. Another three samples (1.8%) showed undetermined or late CT (>31) for at least three of the interrogated miRNAs, and these were excluded from further analysis. Statistical analysis was performed on the remaining 156 samples, and sporadic missings were replaced by the maximum CT value +1 per miRNA. Based on 146 samples (86% of total NAF validation cohort) in the complete regression model including BMI, age, parity and NAF color class, linear regression confirmed that the expression of one of the candidate miRNAs, hsa-miR-29c-5p, was significantly increased in an extremely high MD context (FC 2.1; p = 0.029) (Figure 2a). Firth's corrected logistic regression supported the linear regression (OR = 1.38 (95% CI 1.00-1.98); p = 0.048) (Figure 2b). All four miRNAs were differentially expressed between NAF colors. Hsa-miR-29c-5p and hsa-miR-92b-3p were downregulated in parous versus nulliparous women (FC 0.46; p = 0.019, and FC 0.53; p = 0.016, respectively) while BMI and age had no influence on miRNA expression. from further analysis. Statistical analysis was performed on the remaining 156 samples, and sporadic missings were replaced by the maximum CT value +1 per miRNA. Based on 146 samples (86% of total NAF validation cohort) in the complete regression model including BMI, age, parity and NAF color class, linear regression confirmed that the expression of one of the candidate miRNAs, hsa-miR-29c-5p, was significantly increased in an extremely high MD context (FC 2.1; p = 0.029) (Figure 2a). Firth's corrected logistic regression supported the linear regression (OR = 1.38 (95% CI 1.00-1.98); p = 0.048) (Figure 2b). All four miRNAs were differentially expressed between NAF colors. Hsa-miR-29c-5p and hsa-miR-92b-3p were downregulated in parous versus nulliparous women (FC 0.46; p = 0.019, and FC 0.53; p = 0.016, respectively) while BMI and age had no influence on miRNA expression.
(a) (b) mRNA target analysis was performed to investigate the functional role of hsa-miR-29c-5p (Table 2). Hsa-miR-29c-5p has few established targets as most research so far has focused on the dominant miRNA arm, hsa-miR-29c-3p, and other family members such as miR-29a and miR-29b sharing the same seed region [29]. Established miR-29c-5p targets so far include DNMT3A [30], TMEM98 [31], and CPEB4 [32]. CFLAR, YY1, PTEN, and CD36 are potentially interesting targets in relation to MD based on computational methods. Considering the small number of established targets for this miRNA, no enrichment analysis could be performed. A graphical overview of the validated miRNA and its potential relation to high MD reflected by some of its established and predicted targets is provided in Figure 3. mRNA target analysis was performed to investigate the functional role of hsa-miR-29c-5p (Table 2). Hsa-miR-29c-5p has few established targets as most research so far has focused on the dominant miRNA arm, hsa-miR-29c-3p, and other family members such as miR-29a and miR-29b sharing the same seed region [29]. Established miR-29c-5p targets so far include DNMT3A [30], TMEM98 [31], and CPEB4 [32]. CFLAR, YY1, PTEN, and CD36 are potentially interesting targets in relation to MD based on computational methods. Considering the small number of established targets for this miRNA, no enrichment analysis could be performed. A graphical overview of the validated miRNA and its potential relation to high MD reflected by some of its established and predicted targets is provided in Figure 3. Table 2. Relevant targets of breast density associated hsa-miR-29c-5p. Targets are confirmed by strong (luciferase reporter assays, Western blot, or qPCR) or weaker experimental evidence (next generation sequencing, microarray, or others, annotated with *). Protein class, Gene Ontology (GO) biological process (BP) and Reactome pathways were explored via PANTHER 16.0 pathway analysis (http://www.pantherdb.org/ (accessed on 23 February 2022)) and summarized.

Discussion
Mammographic breast density is one of the strongest, independent risk facto breast cancer [1]. MiRNAs may be involved in the development of extremely high and can be studied in NAF as a representative liquid biopsy of the breast microenv ment. Here, we demonstrated elevated hsa-miR-29c-5p expression in two indepe NAF cohorts acquired from women with extremely high MD, compared to very low Hsa-miR-29c-5p has been identified as one of the most significantly differential pressed miRNAs between ER-positive and -negative breast tumors, demonstrat strong partly GATA3-dependent upregulation in ER-positive/luminal subtype tu compared to normal breast tissue and ER-negative tumors, and which is already obs in pre-invasive breast lesions [33,34]. The fact that miR-29c-5p dysregulation occurs on during carcinogenesis is in line with its elevated expression in the context of extre high MD. Other studies have, however, suggested a tumor suppressor rather than genic role of miR-29-5p in (breast) cancer as this miRNA is generally more high pressed in less aggressive/better prognosis subtypes, and overexpression in cance lines typically inhibits metastasis, proliferation, and migration while increasing apo [31,32,35]. Additionally, many of its established targets so far seem to contribute to a cancerous phenotype [36][37][38][39][40]. Nevertheless, our finding of a higher NAF miR-29c-5 pression in the context of extremely high MD could be a consequence rather than a c aiming to suppress breast cancer in individuals with a higher breast cancer risk d extremely high MD. Or, this miRNA could be selectively released into NAF leadi differing cellular and extracellular miRNA profiles. Interestingly, the miR-29 famil been reported to frequently reside in extracellular exosomes and upregulation of p exosomal miR-29c-5p has recently been suggested as a diagnostic biomarker for early carcinoma [41]. Furthermore, miR-29c-5p was upregulated in bladder serum cancer ples compared to normal samples and its dysregulation was correlated to advanced and poor outcome in bladder cancer patients [42].
To better understand the role of miR-29c-5p in MD, a target analysis was perfo Few robustly validated targets have been described for miR-29c-5p. Its expressi

Discussion
Mammographic breast density is one of the strongest, independent risk factors for breast cancer [1]. MiRNAs may be involved in the development of extremely high MD and can be studied in NAF as a representative liquid biopsy of the breast microenvironment. Here, we demonstrated elevated hsa-miR-29c-5p expression in two independent NAF cohorts acquired from women with extremely high MD, compared to very low MD.
Hsa-miR-29c-5p has been identified as one of the most significantly differentially expressed miRNAs between ER-positive and -negative breast tumors, demonstrating a strong partly GATA3-dependent upregulation in ER-positive/luminal subtype tumors compared to normal breast tissue and ER-negative tumors, and which is already observed in pre-invasive breast lesions [33,34]. The fact that miR-29c-5p dysregulation occurs early on during carcinogenesis is in line with its elevated expression in the context of extremely high MD. Other studies have, however, suggested a tumor suppressor rather than oncogenic role of miR-29-5p in (breast) cancer as this miRNA is generally more highly expressed in less aggressive/better prognosis subtypes, and overexpression in cancer cell lines typically inhibits metastasis, proliferation, and migration while increasing apoptosis [31,32,35]. Additionally, many of its established targets so far seem to contribute to a more cancerous phenotype [36][37][38][39][40]. Nevertheless, our finding of a higher NAF miR-29c-5p expression in the context of extremely high MD could be a consequence rather than a cause, aiming to suppress breast cancer in individuals with a higher breast cancer risk due to extremely high MD. Or, this miRNA could be selectively released into NAF leading to differing cellular and extracellular miRNA profiles. Interestingly, the miR-29 family has been reported to frequently reside in extracellular exosomes and upregulation of plasma exosomal miR-29c-5p has recently been suggested as a diagnostic biomarker for early lung carcinoma [41]. Furthermore, miR-29c-5p was upregulated in bladder serum cancer samples compared to normal samples and its dysregulation was correlated to advanced stage and poor outcome in bladder cancer patients [42].
To better understand the role of miR-29c-5p in MD, a target analysis was performed. Few robustly validated targets have been described for miR-29c-5p. Its expression in breast cancer was found negatively correlated with the mRNA and protein expression of the DNA methyltransferase DNMT3A, which can alter global DNA methylation levels and, amongst others, result in increased collagen type I (COL1A1) promoter activity, contributing to high MD [36]. An interesting predicted miR-29c-5p target based on computational algorithms, assuming an oncogenic role, is the apoptosis regulator protein CFLAR, which was previously identified as one of the plasma proteins negatively associated with area-based breast density [37]. CFLAR mRNA expression is also lower in breast cancer when compared to normal counterpart tissue [43]. MiR-29c-5p might also target PTEN. Low stromal PTEN expression is associated with multiple cancers, but also with high mammographic density through collagen deposition [44]. Murine PTEN knockout models demonstrate more collagen secretion, fiber formation and increased expression of the ECM regulator SPARC [44]. Furthermore, a mouse model with fibroblast-specific PTEN deletion resulted in an alteration of collagen organization, a promotion of collagen orientation in surrounding tumors, and an enhancement of highly aligned matrices when comparing with a wildtype model [45].
Yin Yang 1 (YY1), involved in transcriptional control, chromatin remodeling, and DNA damage repair, might also be targeted by miR-29c-5p. Although the role of YY1 in cancer is controversial, with both oncogenic and tumor suppressor roles [46], YY1 acts as a strong negative regulator of periostin expression [47]. The glycoprotein periostin regulates collagen fibril morphology and potentially LOX-mediated fibril crosslinking. Periostin is overexpressed in raised mammographic density and in most breast cancers, where it enhances angiogenesis and tumor progression, and recruits Wnt ligands to maintain cancer stem cell maintenance [48]. YY1 was however overexpressed in two 3D mammary epithelial cell models that mimic high MD [38].
Mostly based on research that focused on other miR-29 family members having an identical miRNA seed region (the miRNA region that is essential for the binding of the miRNA to the mRNA), one of the most interesting potential targets of hsa-miR-29c-5p in relation to high MD development is CD36. This transmembrane receptor expressed on the surface of breast stromal cells modulates multiple MD-related processes including adipocyte differentiation, fibroblast activation, matrix accumulation, cell-ECM interactions, and immune signaling [2,9,[49][50][51]. CD36 downregulation in vitro causes accumulation of various ECM proteins including collagen type I and fibronectin. In vivo, CD36 knockout mice show less fat accumulation and more ECM accumulation, two prominent phenotypes observed in desmoplasia and high MD [8]. Consequently, low stromal CD36 expression has an established relation to high MD [8].
Future functional studies are needed to confirm suggested targets, to explore the associations with MD phenotype and breast cancer, to determine whether miRNA dysregulation is a cause or a consequence of high MD, or whether miRNA dysregulation and high MD are a consequence of another common factor. This study, for example, also showed elevated miR-29c-5p expression in nulliparous women. It has been demonstrated by other research that parity can affect miRNA expression [52], however there is no in depth reason as to why and how this takes place. Estrogen changes due to pregnancy likely influence the miRNA expression due to genomic and non-genomic mechanisms of action [53], but whether these miRNA changes are temporary has not been investigated.
From a research point of view, there are some limitations to this study. Firstly, only 754 miRNAs out of more than 2600 established human mature miRNAs were investigated [54], and only 118 of these 754 miRNAs were reliably expressed in NAF. Other potentially relevant miRNAs may thus not have been discovered here. In addition, these miRNAs were measured in NAF only, so defining their origin cells was not possible and future studies are necessary to discover the source of the miRNAs. Furthermore, statistical analysis of the discovery cohort was hampered by insufficient power, leading to adjusted strategies that may have affected the results. Future studies therefore require careful external and functional validation of the candidate miRNAs. Lastly, ethnicity and menopausal influences were not considered in this study [55].

Conclusions
We identified and independently validated increased miR-29c-5p expression in extremely high MD NAF. Comparison of the two extreme density categories was hypothesized to result in the biggest difference in miRNA expression. The next step will be to functionally investigate its targets, known or expected, for their involvement in high MD and subsequent breast cancer development. Our results provide new insights into how high MD might arise and why it is associated with an increased breast cancer risk. An ongoing Dutch trial (NTR6162/NL6031) is now examining miR-29c-5p in NAF from breast cancer patients, in relation to MD. Ultimately, this might enable the development of a breast cancer risk classifier for women with high MD, and of targeted therapies to modify this risk.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14153805/s1, Figure S1: Logistic regression analysis with the four candidate differentially expressed human mature miRNAs in the discovery cohort; Figure S2: Unsupervised hierarchical clustering of 41 nipple aspirate fluid samples based on the expression pattern (delta CT) of four miRNAs differentially expressed between extremely high mammographic density (MD) and very low MD categories in the discovery cohort; Figure S3: Pearson correlation analysis between microfluidics-based RT-qPCR profiling results (X axis; CRT value) and regular-volume RT-qPCR individual assay results (Y axis; CT value) (quality control for profiling) for 14 selected miRNAs showing acceptable concordance. Table S1: List of Taqman Advanced assays present in the 754-miRNA panel used for NAF sample profiling; Table S2: Taqman advanced miRNA assays and associated assay IDs used in the study; Table S3: Raw CT values for 4 candidate miRNAs and endogenous control miRNA hsa-miR-125a-5p in the validation cohort; Table S4: Baseline table comparing