Single-cell transcriptional landscape of MG in dairy cows
A total of 28,393 cells (13,463 in NB and 14,840 in AD) were captured and 24,472 high-quality cells (12,432 in NB and 12,040 in AD) were analyzed after quality control (Fig. 1a). Twenty-four clusters were identified based on the marker genes, with the composition of the cell lineages being very different between NB and AD cows (Fig. 1b and c; Table S1). Based on the gene expression profiles, the 24 clusters were categorized into four major groups as follows: immune cells (clusters 0, 1, 6, 8, 14, 15, 18, 21, and 23); epithelial cells (clusters 16, 18, and 20); fibroblast cells (clusters 2, 3, 5, 7, 9, 10, 11, 13, and 20); and endothelial cells (clusters 4, 12, and 22) (Fig. 1d and e).
All immune cells were determined based on the high expression of PTPRC (Fig. 1e). We observed that the immune cells consisted of CD8T, CD4T, NK, Mast, T, macrophages, neutrophils, B, and activated T cells (Table S2). The cells from the MG epithelium were identified based on the high expression of EPCAM (Fig. 1e) and co-expression of KRT19, KRT8, KRT7, and KRT15 (Table S1). We observed clear differences in gene expression profiles between these three epithelial cell types. Specifically, cluster 16 (C16) showed patterns of luminal cells (high expression of KRT19, LTF), basal cells (high expression of KRT5 and KRT15), and their progenitors (ELF5) at the same time (Table S1), indicating that C16 may represent a mixed linage that contains several epithelial cell types. Cluster 17 (C17) was interpreted as mature luminal cells based on the high expression of KRT8, KRT18, and AREG (Table S1). Cluster 19 (C19) expressed not only the marker gene of progenitor cells ELF5 but also the genes associated with milk protein synthesis, such as LTF, CSN2, CSN3, and LALBA (Table S1). Thus, C19 may represent an alveolar luminal progenitor cell or even a mixed lineage. Fibroblast cells expressed genes encoding collagen and extracellular matrix-related molecules, such as COL1A1, DPT, COL3A1, and DLK1 (Fig. 1e). Endothelial cells were identified due to the high expression of PECAM1, APOLD1, CLDN5, ECSCR, and EGFL7 (Fig. 1d, e).
Analyzing these major cell types, we found that the distribution of epithelial and fibroblast cells differed significantly between NB and AD cows (Fig. 1f), which that epithelial cells were mainly observed in the AD, fibroblast cells were mainly identified in the NB. Due to the diversity of epithelial and fibroblast cell types, additional re-clustering was explored.
Reconstruction of mammary epithelium differentiation hierarchy
To better evaluate epithelial cells (C16, 17, and 19), additional re-clustering was performed, and six independent sub-clusters were generated (L0–L5; Fig. 2a and b). All of the sub-clusters featured high expression of the luminal epithelial cell marker genes KRT18, KRT19, and KRT8 (Fig. 2c). Moreover, the expression profiles of the top abundant genes were very similar between L2 and L3. Parts of the top abundantly expressed genes were shared between L4 and L0/L1 (Fig. 2b). Both L2 and L3 were the luminal alveolar progenitor cell subtypes because of the high expression of ELF5, CSN2, and CSN3 (Table S3). L4 showed characteristics of both progenitor cells (ALDH1A3) and mature luminal cells (KRT18, KRT8, AREG; Table S3); thus, it was considered an intermediate luminal cell subtype. The L0, L1, and L5 genes expressed the unique marker gene APOA1, STC2, and PTX3, respectively (Fig. 2c; Table S3).
To further investigate the characteristics of L0, L1, and L5 and understand how these observed cell subtypes and states are related, the pseudotemporal analysis of single cells was performed using Monocle2. We found one tightly connected differentiation trajectory that was separated into two main branches, of which the start points of pseudotime were located within L2 and L3 cells (Fig. 2d). L0, L1, and L5 were present at the end of the two branches (Fig. 2d), indicating these clusters are different subtypes of mature luminal epithelial cells. Except for expressing the luminal cell marker gene KRT19, the L0 subtype also expressed LTF, as well as KRT5 and KRT15 (marker genes of basal cells) (Table S3). Therefore, the L0 subtypes showed the properties of differentiated luminal cells. L1, located at the end of the other branch of the epithelial differentiation hierarchy, showed the typical markers of mature luminal cells such as KRT18, KRT8, AREG, and PRLR (Table S3). It was considered a mature luminal cell and was defined as a secretory alveolar cell. L5 expressed PTX3, which responds to inflammatory stimuli in epithelial cells and regulates the innate response to pathogens and inflammatory reactions. Thus, we defined this subtype as an immune-sensing luminal cell. Additionally, the immunofluorescence analysis was used to validate the three subtypes at the protein level, to integrate the newly discovered cell subtypes spatially, and to create the single-cell atlas of the cattle MG (Fig. 2e).
Characteristics and functions of epithelial cell subtypes
Next, we identified 682 genes with branch-specific expression patterns, as well as the key genes and biological processes changing over pseudotime (Fig. 3a). Towards the L1 branch, we found three gene clusters that increased in expression during differentiation. These clusters consisted of key genes that were involved in the development of the MG alveolus and responded to amino acid, mineral, and hormone synthesis (Fig. 3a, Table S4). Towards the L0 branch, one cluster of highly expressed genes was involved in lipid function (e.g., isoprenoid, sterol, steroid, and cholesterol), biosynthetic processes, and regulation of epithelial cell differentiation and proliferation (Fig. 3a, Table S4). We also observed one cluster of genes that gradually decreased expression in both branches, which were involved in response to progesterone and estradiol and the prolactin signaling pathway (Fig. 3a, Table S4).
Five other co-expressed gene modules were generated using WGCNA to evaluate the gene co-expression patterns of L0, L1, and L5. There were 412, 209, 1682, 503, and 21 genes in the blue, brown, grey, turquoise, and yellow modules, respectively (Fig. 3b, Table S5). We found that the turquoise, brown, and yellow modules showed a highly positive correlation with the L0 (R = 0.81), L1 (R = 0.88), and L5 (R = 0.44) subtypes, respectively (Fig. 3b). Next, we performed GO term enrichment analysis to explore the biological processes of the genes in these three key modules (Fig. 3c). The results showed that genes from the yellow module were involved in cytoskeleton organization and immune response. The genes from the brown module were mainly associated with MG epithelial cell differentiation and alveolus development, lactation, and protein, hormone, and mineral metabolism. Finally, the turquoise module genes were mainly involved in epithelium development, epithelial cell differentiation, and the regulation of lipid biosynthetic and metabolic processes, including sterol, steroid, cholesterol, and fatty acid metabolism (Fig. 3c).
Reconstruction of fibroblast cell subtypes
To further reveal the subtypes of fibroblast cells, additional clustering was performed and 14 independent sub-clusters were generated (F0–F13; Fig. 4a and b). Based on the marker gene expression profiles, we defined F0, 3, 6, 7, and 10 as lipofibroblasts and F4 as the matrix fibroblast cell. F5, 8, 13 showed the characteristics of myofibroblasts. However, the other five fibroblast cell subtypes could not be defined (Table S6). Interestingly, we found that F13 showed strong characteristics of epithelial and progenitor cells, which was reflected by the highly expressed marker genes such as EPCAM, KRT18, KRT8, ELF5, and CSN3 (Fig. 4c). Immunofluorescence staining confirmed the existence of F13 with its marker genes COL1A1 and CSN3 (Fig. 4d). In addition, F13 was mainly involved in cell differentiation; response to estradiol, progesterone, and estrogen; cell migration; apoptosis; and growth according to the GO analysis (Fig. 4e).
Investigation of luminal and fibroblast cell differentiation hierarchy
For a comprehensive understanding of the relationships between fibroblasts and epithelial cells, a diffusion map was used to computationally reconstruct differentiation stages and the transitions between luminal epithelial and fibroblast cell subtypes. We found that the epithelium and fibroblasts were not completely separated, and ~60% of F13 cells were clustered with epithelial cells (Fig. 5a). The F13 subtype showed a strong relationship with L1 (R = 0.871) and L5 (R = 0.928) based on Pearson correlation (Fig. 5b). In addition, the epithelial cells (L1–L5) and fibroblasts (F12) interacted with F13 through the ligand-receptor COL1A1-ITGA2B, IGF2-IGF1R/INSR, HMGB1/THBS1-SDC1, and RBP4-STRA6 (Fig. 5c).