LncRNA as ceRNAs may be involved in lactation process

The main function of the mammary gland is to secret milk for newborn growth. Milk production process is regulated by hormones, growth factors, noncoding RNAs and other factors locally. Long non-coding RNAs (lncRNAs), one type of recently discovered non-coding RNA, have been found in mammary gland and some studies suggested lncRNA may play important roles in mammary gland development. Competing endogenous RNAs (ceRNAs) are emerging to compete for miRNA binding and, in turn, regulate each other. In the current study, we sequenced mRNA, miRNA and lncRNA in goat mammary tissue at 2 points in lactation (early and mature). All data were co-expressed together from the same samples. Our data showed that the ceRNAs up-regulated during the mature lactation phase were associated with lipid, protein, carbon and amino acid synthesis and metabolism. This correlates with the function of the mature lactation phase: i.e. the continuous production of large amounts of milk, rich in proteins, lipids, amino acids and other nutrients. Alternately, the ceRNAs up-regulated during early lactation were associated with PI3K-AKT pathways and ECM-receptor interactions; these fulfil the functional role of preparing the mammary gland for full lactation. Therefore, the results suggest that ceRNAs work synergistically during different developmental stages to regulate specific functions associated with lactation control. This study suggests that ceRNAs (lncRNA-mRNA) may be involved in lactation process.


INTRODUCTION
Milk contains crucial biologically active components for infant growth and development, including immunological, gastrointestinal, neural and even intellectual development [1]. Mammary gland is a complex gland that originates during the embryonic stage and develops quickly from the pubertal stage with ductal growth and early alveolar development during menstrual cycles; mammary tissue then undergoes proliferation, differentiation and death during pregnancy, lactation and involution stages, respectively [2]. The main function of the mammary gland is to secret milk for infant nutrition [3][4][5][6][7][8][9][10].
The regulation of milk synthesis impacts on the health of both the mother and the neonate [11][12][13]. Anderson et al. [13] reports four phases for the mammary gland according to its function: proliferation phase during early pregnancy; secretory differentiation phase starting from midpregnancy; secretory activation phase beginning around parturition; and lactation phase associated with continued www.impactjournals.com/oncotarget/ Oncotarget, 2017, Vol. 8, (No. 58), pp: 98014-98028 Research Paper milk production. Lactation itself includes 2 phases: early phase (un-mature phase) when milk is composed of large amounts of immuno-factors however with less volume; and the mature secretion phase when large volumes of milk are continuously produced to support growth and development of the newborn [13]. Milk production is an evolutionary necessity and this process is tightly regulated at a local level, i.e. within the mammary gland itself, by hormones and other factors [7][8][9][10].
Elucidation of the molecular events of mammary gland development during lactation can aid understanding of molecular development from theearly phase (un-mature phase) to the mature secretion phase and could help to improve milk production [14]. Therefore, identifying the regulatory principles that govern mammary gland development during lactation has attracted the attention of developmental, molecular biologists and even clinicians [14]. MicroRNAs (miRNAs), a class of small noncoding RNAs, have been broadly investigated in mammary gland development due to their inhibitory effects on their target genes [15][16][17][18]. Recently, competing endogenous RNAs (ceRNAs), sharing miRNA recognition elements (MREs), are emerging to compete for miRNA binding and, in turn, regulate each other. Interactions of ceRNAs have potential roles in gene expression and cell phenotypes [19,20]. Nowadays, long non-coding RNAs (lncRNAs), one type of recently discovered non-coding RNA, act as ceRNAs to regulate gene expression [21][22][23]. Furthermore, the important roles of lncRNAs in mammary gland development are becoming increasingly evident [24][25][26][27][28]. Standaert et al. [28] identified that the lncRNA Neat1 plays an important role in mammary gland cell proliferation during lobular-alveolar development. However, these studies simply investigated the expression of lncRNA in the mammary gland but no ceRNAs of lncRNA-mRNA (based on miRNA) were explored. Moreover, the roles of lncRNA associated-ceRNAs in lactation, particularly in the transition from the early (un-mature) phase to the mature secretion phase, are not fully understood. Moreover, the dynamic patterns of ceRNA interactions during this developmental process remain unknown. To better understand the dynamic process of mammary gland development during early lactation, RNA-sequencing of mammary tissue at different time periods during lactation was performed to analyse mRNA, lncRNA and miRNA expression. Furthermore, ceRNA networks (ceR NETs) were established for protein coding genes and lncRNAs.

Morphologic changes in mammary gland tissue and changes in expression of mRNAs, miRNAs and lncRNAs during lactation
Histological changes in goat mammary gland morphology during early lactation (L-5 d) and peak lactation (L-30 d) are shown in Figure 1A. At L-5 d, which represented the colostral phase, the luminal space of the lobuloalveolar structures expanded and the epithelial cell layer was prominent against the adipocytes ( Figure 1A left). At L-30 d, which represented the mature secretion phase, the goat mammary gland was stably producing copious amounts of milk and exhibited prominent luminal structures and ducts with only a few adipocytes ( Figure 1A right). Adipocytes are known to become de-lipided with no change in cell number [13,29]. Goat milk contains about 3.5% protein (caseins, whey acidic protein (WAP), lactoferrin, secretory immunoglobulin A, and others), 3.8% lipid and 4.1% lactose, a disaccharide that is unique to milk [30,31].
To systematically detect the expression of coding genes (mRNA), microRNAs and lncRNAs during goat lactation, we first analysed the expression profiles of RNAs by RNA-seq and small RNA sequencing for goat mammary gland tissues from L-5 d and L-30 d (Supplementary Table 1, Supplementary Table 2,  Supplementary Table 3 and Supplementary Table 4). In order to search for the genes, miRNAs or lncRNAs related to lactation, the data at L-5 d was considered as basal (control) and compared with data at L-30 d. In total, 42 690 mRNAs were found in mammary gland tissues of which 1150 genes were differentially expressed at L-30 d (549 down-and 601 up-regulated; Figure 1B). The expression of 18 genes was confirmed by real time quantitative RT-PCR ( Figure 1E). The q-RT-PCR data and the RNA-seq data were in the similar trend. In total, 391 known miRNAs were detected in the mammary tissues of which 24 known miRNAs were differentially expressed (15 down-and 9 up-regulated; Figure 1C; Table 1). The expression of 4 miRNAs was confirmed by real time quantitative RT-PCR ( Figure 1F). The q-RT-PCR data and the RNA-seq data were similar. In addition, 985 novel miRNAs were found in the mammary tissue of which 7 novel miRNAs were down-and 10 differentially upregulated (Supplementary Table 5, Supplementary Table  6). In total, 1451 lncRNAs were detected in the mammary tissue of which 33 were differentially expressed at L-30 d (5 down-and 28 up-regulated; Figure 1D; Table 2). The expression of 5 lncRNAs was confirmed by real time quantitative RT-PCR ( Figure 1F). The q-RT-PCR data for lncRNAs was similar to the RNA-seq data.
Moreover, we performed GO and KEGG pathway analyses for the differentially expressed mRNAs for downor up-regulated genes separately ( Figure 2; Supplementary Table 7, Supplementary Table 8, Supplementary Table 9  and Supplementary Table 10). In the up-regulated gene cluster, the PPAR signalling pathway, biosynthesis of amino acid, fatty acid biosynthesis, fatty acid metabolism, steroid biosynthesis, and TCA cycle were the most enriched pathways, with a greater number of genes in the top 20 KEGG pathways ( Figure 2B; Supplementary Table 7). In the GO biological process category, the top 10 classes were the cholesterol biosynthetic process, fatty acid biosynthetic process, isoprenoid biosynthetic process, malonyl-CoA biosynthetic process, positive regulation of protein polymerization, 2-oxoglutarate metabolic process, acetyl-CoA metabolic process, sterol biosynthetic process, geranyl diphosphate biosynthetic process and farnesyl diphosphate biosynthetic process, which matched well with the KEGG pathways ( Figure 2C; Supplementary  Table 9). However, the PI3K-AKT signalling pathway, MAPK signalling pathway, EMC receptor interaction and p53 signalling pathway were the most enriched pathways with a greater number of genes in the top 20 KEGG pathways of the down-regulated genes ( Figure 2D; Supplementary Table 8). Similarly, the top 10 classes in the GO biological process category of the down-regulated genes include regulation of mRNA stability, skeletal system morphogenesis, sphingoid catabolic process, defence response to bacterium, ventricular zone neuroblast division, fibroblast growth factor receptor signalling pathway involved in orbitofrontal cortex development, dendrite regeneration, sphingosine metabolic process, ureteric bud development and cyclooxygenase pathway ( Figure 2D; Supplementary Table 10). Our data partially matched early data from bovine mammary gland lactation stages [32], in both bovine milk [33] and goat milk [34]; this may be because of the different stages or different species of animals investigated in the different studies.

Competitive interactions among coding-RNAs and lncRNAs during lactation
To reduce background and non-specific interactions, only differentially expressed mRNAs (1150 genes), miRNAs (n =24) and lncRNAs (n = 33) were used to explore target competition during lactation. We used the miRanda algorithm to obtain a prediction of miRNA regulation of protein-coding RNAs (mRNAs) and lncRNAs (see Materials and Methods). It was found that all the differentially expressed genes (1150) were predicted to be regulated by the 24 differentially expressed miRNAs with 17 989 of miRNA recognition elements (MREs; Figure 3A). Then, interactions between miRNAs and mRNAs were predicted based on their correlation in gene expression and 10 859 MREs were found. The shared MREs between predicted MREs based on miRNA binding sites and MREs predicted from the expression were considered as valuable MREs for further ceRNAs searches (see Materials and Methods). The shared number of MREs for miRNAs-mRNAs was 4032 ( Figure 3A). Similarly, the shared number of MREs between the predicted MREs (415) and calculated MREs (269) for miRNAs-lncRNAs was 59 ( Figure 3A). The shared 4032 miRNAs-mRNAs and 59 miRNAs-lncRNAs were then used to build the ceRNAs. In this investigation, we aimed to investigate lncRNAs as ceRNAs (lncRNA-mRNA pair) in the lactation process. The number of predicted ceRNA based on ceRNA score was 242 from using the shared miRNA-mRNAs and shared miRNA-lncRNAs (Supplementary Table 11). The number of calculated ceRNAs based on the expression of miRNAs, mRNAs and lncRNAs was 5712. The number of shared ceRNA between predicted ceRNAs and calculated ceRNAs was 164, and these were the ceRNAs (lncRNA-mRNA) used in this investigation for further analysis ( Figure 3A) [35][36][37]. In mammary gland cells, lncRNAs are dynamically expressed and they function during mammary gland development and even in breast cancer [24][25][26][27][28]. However, lncRNAs that act as ceRNAs have not yet been implicated in mammary gland development, especially during the lactation process. In our study, at the transit stage of lactation, 33 lncRNAs were differentially expressed and   most of the lncRNAs (n = 28) were up-regulated at L-30 d. As demonstrated by the clustering in Figure 3B, 1/3 of the 164 ceRNAs were up-regulated at L-30 d and 2/3 of the ceRNAs were down-regulated in this phase ( Figure  3B). To predict the functions of ceRNAs during lactation, we determined the enrichment of GO terms and KEGG pathway analysis for coding genes directly interacting with lncRNAs. In the GO biological process category, the most significant classes were fatty acid biosynthesis, steroid biosynthesis, lateral sprouting from an epithelium; and in KEGG pathway analysis the most significant pathways included ECM-receptor interaction, steroid biosynthesis, protein digestion and absorption, glycerolipid metabolism, fatty acid biosynthesis and carbon metabolism, which support the findings that lipid metabolism, steroid metabolism and EMC-receptor interaction are important in the lactation process ( Figure 3C & 3D) [11][12][13].
In order to search the interaction of ceRNAs in the ceRNA network, we further analysed ceRNAs by comparing expression levels. We found that lncRNAs TCONS-00055666, TCONS-00144689 and TCONS-00108242 were up-regulated at L-30 d and their ceRNAs paired mRNAs (genes) THBS1, FGFR2, CACNA1C, LAMA2, COL4A5 and PPKAA1 were down-regulated during the same phase. Functions of these genes include EMC-receptor interaction, PI3K-Akt signalling pathway and MAPK signalling pathway ( Figure 4A & 4C). This agrees with the mammary gland morphology change stated earlier and also matches the functional changes of these 2 phases [11][12][13]. On the other hand, lncRNAs TCONS-00144434, TCONS-00148514 and TCONS-00055666 were up-regulated at L-30 d and their ceRNAs paired mRNAs (genes) FASN, LPL, GPAM SC5D and MSMO1 were also up-regulated in this phase. The functions of these genes include fatty acid biosynthesis, sterol biosynthesis and positive regulation of cholesterol storage ( Figure 4B & 4C). This agrees with previous findings that the mammary gland produces a continuous supply of milk for the infant and that fatty acids are one of the most important components of milk in this phase [11][12][13]. The ceRNA network of lncRNA-mRNA was present in Figure 4C, where the nodes represented mRNAs or lncRNAs and the edges represented their competing interactions. LncRNA TCONS_00055666 had the most number of mRNAs as ceRNA, then less number of mRNAs were as the ceRNAs for lncRNAs TCONS_0088402, TCONS_00108242, TCONS_00144689, TCONS_00158176 and TCONS_ 00148514. There were just 2 mRNAs as ceRNAs for lncRNA TCONS_00144434 or TCONS_00146083. The network of ceRNA (lncRNA-mRNA) and miRNA was shown in Figure 5A, where the nodes represented mRNAs, lncRNAs or miRNAs and the edges represented their competing interactions. miR-200b and miR-200c formed more interaction with mRNA or lncRNA than other miRNAs. These networks suggested that miRNAs, lncRNAs and miRNAs interact together in the lactating mammary gland. They might be involved in the lactation process.
Since, ceRNAs play important roles in lactation, next we further analysed protein levels of a few important regulators in lactation: FASN, INSIG1 and COL4A5, which were found to be differentially expressed. In agreement with ceRNA (lncRNA-mRNA) expression patterns, COL4A5 was down-regulated at L-30 d ( Figure  5A & 5B), while FASN and INSIG1 were up-regulated at L-30 d ( Figure 5A & 5C) [13].

DISCUSSION
In this investigation, we determined the ceRNAs of lncRNA-mRNA during different phases of lactation through the MREs in miRNAs. Over the past decades, miRNAs have been extensively explored in many organisms, tissues and even in mammary gland development. Studies suggest that miRNAs play important roles in different stages of mammary gland development [38,39]. It was found that some important miRNAs were up-regulated during lactation [38][39][40]. In our current investigation, we found 9 miRNAs were up-regulated at L-30 d compared with early lactation (L-5 d).
Compared to miRNAs, only a few studies have explored lncRNAs in mammary gland development. First reported was lncRNA Pinc, which is a pregnancy-induced lncRNA [41]. Other lncRNAs related to mammary gland development include Zfas [42], HOTAIR and GAS5-AS1 [40], Ror and Neat1 [28]. Taken together, these studies suggest that lncRNAs as well as miRNAs play a crucial role in regulating mammary epithelial cells. It has been reported that the expression of lncRNAs is more tissue-specific than protein coding genes [43], and the competitive interactions of lncRNAs are more tissuespecific [44]. These were the fundamental principles used in our studies and we identified 33 lncRNAs that were differentially expressed in the 2 phases of lactation and 28 of them were up-regulated at L-30 d.
MicroRNAs are commonly considered as active regulatory elements to decrease the stability of target RNAs or to block their translation [35]. However, Poliseno et al., demonstrated that pseudogenes acting as microRNA competitors thereby actively competed with their ancestral protein-coding genes for the same pool of microRNAs through sets of conserved MREs [45]. This suggests that non-coding RNA targets can cross-talk through their ability to compete for microRNA binding [35,45]. Therefore, MREs can be used by the transcripts to actively communicate with each other and to regulate their respective expression levels [35]; furthermore, MREs may functionalize the entire mRNA dimension through cross-talking ceRNAs and even ceRNA networks [35].
In previous studies, only lncRNA and mRNA were analysed in the same samples, and miRNAs were computed from other databases to determine the ceRNA of lncRNA-mRNA, mRNA-mRNA and lncRNA-lncRNA [44]. In the current study, we sequenced mRNA, miRNA and lncRNA in the same samples at 2 phases of lactation. All data were co-expressed together from the same samples. Furthermore, we constructed the ceRNA and ceRNA network using the expression data for the same samples. This is the first study to be performed in this way. Even though a lower number of ceRNAs were found in our study, the data is more accurate and more reliable because all the expression of mRNA, miRNA and lncRNA were correlated. In total, 164 ceRNAs were constructed in this study and further functional analysis suggested that the dynamic expression and regulation pattern were related to lactation. As early studies indicate that highly expressed genes in particular developmental stages or cell types are functionally active, our data showed that 1/3 of the ceRNAs were up regulated during the mature lactation phase with the functions of lipid, protein, carbon, amino acid synthesis and metabolism. This correlates with the function of the mature lactation phase, i.e. the continuous production of a large amount of milk rich in protein, lipids, amino acids and other nutrients. Furthermore, 2/3 of ceRNAs were up-regulated in early lactation with functions in PI3K-AKT pathways or ECM-receptor interactions. This correlated with the function of this lactation phase, i.e. to prepare the mammary gland for full lactation. Therefore, the results indicated that ceRNAs work synergistically during different developmental stages to regulate specific functions, particularly to control lactation in the mammary gland.
In summary, by applying ceRNA network analysis of transcriptomes (mRNA, miRNA, lncRNA) obtained during lactation, this study provides a novel approach for understanding gene functionality in this biological process. The investigation yielded many new insights into the structure of molecular pathways in the lactation process. Analysing ceRNA interactions in the context of lactation may provide insights into the regulation of lactation. Moreover, this study suggests that ceRNAs (lncRNA-mRNA) may be involved in lactation process.

Mammary gland tissue collection and morphological analysis
This investigation was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Committee on the Ethics of Animal Experiments of Qingdao Agricultural University IACUC (Institutional Animal Care and Use Committee). Six 2-year-old Laoshan dairy goats, one of 4 native Chinese dairy goat breeds, were used in this study. Goat mammary gland tissue from lactation day 5 (L-5 d; early lactation) and lactation day 30 (L-30 d; mature lactation) were surgically collected and frozen immediately in liquid nitrogen for further analysis (3 animals/group) [46]. Part of the mammary gland samples were fixed in 10% neutral formalin and then paraffin embedded. Subsequently, 5-μm sections were prepared and stained with hematoxylin and eosin (H&E) for the analysis of morphological changes [47,48].

Sequencing
Total RNA was extracted from individual mammary gland samples and cDNA libraries were constructed. For small RNA cDNA library, total RNA (about 1 μg) was ligated with RNA 3' Adapter and 5' Adapter, then the ligated RNAs were reversely transcribed to cDNAs using RT primers. The cDNAs were amplified by PCR, purified by gel electrophoresis. The quality was determined by Agilent 2100 chip. For RNA-seq analysis library (mRNA and lncRNA), total RNA was purified to remove rRNA by Ribo-Zero™ rRNA Removal Kit, then RNA was fragmented. The fragmented RNA was reversely transcribed into the first strand cDNA by TruSeq® Stranded kit, then DNA polymerase I and RNAase H were used to form the double strands cDNA. The double stranded cDNA was adenylate 3' Ends and ligated Adapters. Then the RNA was amplified by PCR and purified to form the cDNA library. Libraries were sequenced using the Illumina HiSeq 2500 platform with the 90-bp pair-end sequencing strategy for total RNA sequencing and Illumina-HiSeqXTen platform for small RNA sequencing. The original image data generated by the sequencing machine were converted into sequence data by Illumina pipeline CASAVA v1.8.2 and then subjected to standard QC criteria to remove the contaminant following parameters [44]: (1) Reads aligned to adaptors or primers with no more than 2 mismatches.

Read mapping and gene expression analysis
The reference genome, transcriptome and annotation reference for the goat is located on the NCBI platform (https://www.ncbi.nlm.nih.gov/genome/?term=Capra%20 aegagrus%20hircus/genome.fa.gz; https://www.ncbi. nlm.nih.gov/genome/?term=Capra%20aegagrus%20 hircus/transcriptome.fa.gz; https://www.ncbi.nlm.nih. gov/genome/?term=Capra%20aegagrus%20hircus/ gff.gz). Clean paired-end reads were aligned to the reference genome using TopHat (version 2.0.6) [49]. The transcriptome of each sample was constructed using Cufflinks (version 2.0). Transcripts >200 nt were identified as lncRNAs if they did not overlap with known genomic annotations from the database [49]. miRNAs were analysed using Bowtie software and the miRBase database. Next, the reads per kilo base of model per million base pairs sequenced (RPKM) was used to quantify the expression levels of a gene or lncRNA [44] and transcript per million (TPM) was used to determine the expression levels of miRNA [37]. The difference www.impactjournals.com/oncotarget between different group was determined by DEGseq software (Fold Change≥2.00 and FDR≤ 0.001). miRNA target prediction miRanda database was used with the default parameters to identify conserved miRNA target sites in the 3'UTR of coding transcripts and full-length lncRNA transcripts [37]. If multiple annotated 3'UTR/lncRNA sequences were found for a coding/lncRNA gene, the longest was used in the analysis. These were the predicted pairs for miRNA target of the differentially expressed miRNAs, mRNAs and lncRNAs in the 2 goat mammary gland phases studied.
Construction of the ceRNA NET work related to lactation phase [44,50] Overview of the processes used to identify ceRNA interaction pairs Based on the expression levels of mRNAs, miRNAs and lncRNAs, Pearson's correlation coefficient and P-value were calculated for miRNA-target. Negatively correlated pairs with a P-value <0.05 and Pearson's correlation coefficient >0.7 were selected for further analysis. These were the predicted pairs for miRNA-mRNA, miRNA-lncRNA and mRNA-lncRNA of the differentially expressed miRNAs, mRNAs and lncRNAs in the 2 goat mammary gland phases. Subsequently, shared pairs from the predicted pairs from binding sites and the predicted pairs from the expression of mRNA, lncRNA and miRNA were used for further analysis.
Based on the principle for ceRNA prediction, shared pairs of miRNA-mRNA and miRNA-lncRNA were used to predict ceRNA score according to the following formula [50]. ceRNA_score #MRE_for_share_miRNA #MRE_for_lncRNA_miRNA = P-value was calculated as follows [50]: The shared pairs from the predicted pairs of lncRNA-mRNA based on the expression of lncRNA and mRNA (Pearson's correlation coefficient) and the predicted pairs of lncRNA-mRNA based on ceRNA score principle were then considered as the true ceRNAs (lncRNA-mRNA).

Enriched functional categories for ceRNA [44]
The functional annotations of genes for ceRNAs (lncRNA-mRNA) were obtained from the GO and KEGG databases

Construction of the ceRNA network
The ceRNA network was constructed by assembling all the ceRNA pairs, where the nodes represented mRNAs or lncRNAs and the edges represented their competing interactions. The network for ceRNA and miRNAs was also formed, where the nodes represented miRNAs and the edges represented the ceRNA and miRNA interactions

Real-time quantitative RT-PCR
The procedure for mRNA q-RT-PCR was reported in our early publications [51]. MiRNA q-RT-PCR was performed using the Hairpin-it™ miRNA RT-PCR (probe) Quantitation kit from GenePharma Co., Ltd (Shanghai, China) following the manufacturer's instructions, as described in our recent publication [52]. Similarly, lncRNA q-RT-PCR was performed using the Hairpin-it™ lncRNA RT-PCR (probe) Quantitation kit from GenePharma Co., Ltd following the manufacturer's instructions. RNA from mammary gland tissues was extracted using TRIzol Reagent (Invitrogen Corp., Carlsbad, CA, USA) and purified using an RT2 qPCR-Grade RNA Isolation Kit from SABiosciences Co., Ltd (MD, USA). Total RNA was quantified using a Nanodrop 3300 (ThermoScientific, DE, USA). The quality of RNA was controlled by the A260:A280 ratio being greater than 2.0 and confirmed by electrophoreses, with a fraction of each total RNA sample with sharp 18S and 28S ribosomal RNA (rRNA) bands as reported in our recent publication [52]. One microgram of total RNA was used to make the first strand cDNA in 20 μl. The program for the reaction of miRNA and lncRNA was 25 °C for 30 min, 42 °C for 30 min, 85 °C for 5 min, then 4 °C or on ice. The qPCR was performed with the Roche LightCycler 480 (Roche, Germany) and the reaction was as follows: Step 1, 95 °C for 3 min; Step 2, 40 cycles of 95 °C for 12 s; 62 °C for 40 s. The primer sets for mRNA, miRNA and lncRNA are given in Supplementary Table 12, 13 and 14. Three independent experimental samples were analysed [52]. q-RT-PCR was statistically analysed using proprietary software from SABiosciences online support (www. SABiosciences.com).

Immunofluorescence staining
Mammary gland sections (5 μm) were prepared and subjected to antigen retrieval and immunostaining as previously described [46]. Briefly, sections were first blocked with normal goat serum in PBS, followed by incubation (1:150 in PBS-1% BSA) with rabbit anti-FASN Ab (Cat no.: bs-1498R) and rabbit anti-INSIG1 Ab (Cat no.: bs-5074R) from Beijing Biosynthesis Biotechnology Co., Ltd (Beijing, China) and rabbit anti-COL4A5 Ab (Cat no.: sc-9302) from Santa Cruz Co., Ltd (CA, USA) overnight at 4 °C. The following morning, after 3 washes with PBS Tween 20 (0.5%), the slides were incubated with Alexa Fluor 546 goat anti-rabbit IgG (1:200) for 30 min in darkness at RT. The negative control samples were incubated with secondary antibody and without primary antibody. Slides were washed thrice with PBS Tween-20 and then incubated with DAPI (4.6-diamidino-2-phenylindole hydrochloride, 100 ng/ml) as a nuclear stain for 5 min. After a brief wash with ddH 2 O, the slides were covered with an anti-fading mounting medium from Vector Co., Ltd (CA, USA). Fluorescent images were obtained using a Leica Laser Scanning Confocal Microscope (LEICA TCS SP5 II, Germany) [46].

Data access
The raw sequencing data generated in this study has been uploaded to the NCBI SRA database with the accession number: PRJNA361394 (http://www.ncbi.nlm. nih.gov/bioproject/PRJNA361394/).

Author contributions
YZ and LM provided key intellectual input in the conception and design of these studies and aided in the writing of this manuscript. SY, FL, MChu, YH and YF performed the experiments. HZ, JL, MC, LL, and WS provided expertise for data explanation and contributed to the writing of the manuscript. All authors read and approved the final manuscript.