Transcriptome analysis to identify candidate genes related to mammary gland development of Bactrian camel (Camelus bactrianus)

Introduction: The demand for camel milk, which has unique therapeutic properties, is increasing. The mammary gland is the organ in mammals responsible for the production and quality of milk. However, few studies have investigated the genes or pathways related to mammary gland growth and development in Bactrian camels. This study aimed to compare the morphological changes in mammary gland tissue and transcriptome expression profiles between young and adult female Bactrian camels and to explore the potential candidate genes and signaling pathways related to mammary gland development.

Methods: Three 2  years-old female camels and three 5  years-old adult female camels were maintained in the same environment. The parenchyma of the mammary gland tissue was sampled from the camels using percutaneous needle biopsy. Morphological changes were observed using hematoxylin-eosin staining. High-throughput RNA sequencing was performed using the Illumina HiSeq platform to analyze changes in the transcriptome between young and adult camels. Functional enrichment, pathway enrichment, and protein–protein interaction networks were also analyzed. Gene expression was verified using quantitative real-time polymerase chain reaction (qRT-PCR).

Results: Histomorphological analysis showed that the mammary ducts and mammary epithelial cells in adult female camels were greatly developed and differentiated from those in young camels. Transcriptome analysis showed that 2,851 differentially expressed genes were obtained in the adult camel group compared to the young camel group, of which 1,420 were upregulated, 1,431 were downregulated, and 2,419 encoded proteins. Functional enrichment analysis revealed that the upregulated genes were significantly enriched for 24 pathways, including the Hedgehog signaling pathway which is closely related to mammary gland development. The downregulated genes were significantly enriched for seven pathways, among these the Wnt signaling pathway was significantly related to mammary gland development. The protein–protein interaction network sorted the nodes according to the degree of gene interaction and identified nine candidate genes: PRKAB2, PRKAG3, PLCB4, BTRC, GLI1, WIF1, DKK2, FZD3, and WNT4. The expression of fifteen genes randomly detected by qRT-PCR showed results consistent with those of the transcriptome analysis.

Discussion: Preliminary findings indicate that the Hedgehog, Wnt, oxytocin, insulin, and steroid biosynthesis signaling pathways have important effects on mammary gland development in dairy camels. Given the importance of these pathways and the interconnections of the involved genes, the genes in these pathways should be considered potential candidate genes. This study provides a theoretical basis for elucidating the molecular mechanisms associated with mammary gland development and milk production in Bactrian camels.

1. Introduction

Bactrian camels (two-humped camels), are mammals of the genus Camelidae of the family Mammalia that mainly inhabit cooler areas in Asia (1). There are five local species of Bactrian camels in China, namely Junggar, Tarim, Alashan, Sunit, and Qinghai, which are distributed in the arid and semi-arid regions of the Xinjiang, Inner Mongolia, Gansu, and Qinghai provinces in northwestern China. The Bactrian camel grazes a wide range of desert plants, can tolerate extreme environments, and has a specialized immune system. It is a distinctly important economic animal resource in desert areas, providing livestock products such as milk, meat, hides, and wool to herders in remote areas. In addition, camel milk has a high nutritional value, is hypoallergenic, and contains small-molecule nano-antibodies that can assist in the treatment of multifunctional disorders and autoimmune diseases such as lung cancer and diabetes (2, 3). However, camels grow slowly, reaching puberty at a later age than other livestock species. Sexual maturity is probably reached at 3–4 years. Female Bactrian camels give birth to their first calves at the age of 4 or 5 and produce offspring once every two years with a gestation period of 13–14 months. This slow reproductive rate, combined with the small geographical distribution, free-range grazing, and traditional breeding, results in low milk production and a shortage of camel milk supply to the market. Camels are crucial to the economies of many countries in arid and semi-arid regions of the world and camel breeding is thus increasing annually worldwide. With this increase, the metabolic processes and molecular mechanisms related to adaptation to extreme desert environments have been identified by genome sequencing and comparative genome analysis of Bactrian camels and dromedaries (4). There is a crucial need for developing knowledge of the genetic and molecular mechanisms related to milk production traits in Bactrian camels. Elucidating the mechanisms of mammary gland development is an initial key step to achieve this.

The mammary gland develops and produces milk under the regulation of systemic hormones, which is a dynamic and highly complex multistep process involving the periodic cycling of mammary epithelial proliferation, differentiation, and apoptosis (5). The components of milk are derived from blood-circulating nutrients, and they are synthesized in the epithelial cells of the mammary gland in lactating animals. Almost all studies on milk production traits have been conducted based on analysis of mammary glands and animal hematology (6). Mammary tissues undergo rapid growth during adolescence and develop again during pregnancy until lactation is complete (7). Animal studies have revealed apparent anatomical and physiological differences between juvenile and adult mammary glands. Therefore, the development and function of the mammary glands are essential for the provision of nutrition to offspring and the study of milk production traits. However, there are limitations in the analysis of milk performance traits using mammary tissues, such as difficulties in sampling, damage to animas, and ethical considerations.

Functional genomic tools have been widely used to study the molecular mechanisms of growth, development, and production traits in livestock (8, 9). Transcriptome analysis has been used to identify the molecular information that regulates economically important traits of an organism. Recently, researchers have studied the molecular mechanisms underlying mammary gland development in various species, such as cattle (10, 11), sheep (12), cats (13), rabbits (13), and rats (14). However, to date, the mechanisms underlying mammary gland development and the initiation of lactation have not been reported in Bactrian camels. The effects of specific candidate genes on physiology, mammary gland development and milk production traits remain unknown. Therefore, the current study uses transcriptome analysis of the mammary gland tissues of young and adult camels to investigate important candidate genes and pathways involved in mammary gland development in Bactrian camels.

2. Materials and methods

2.1. Experimental camels and sample collection

The transcriptome sequencing sample population consisted of 6 domestic Junggar Bactrian camels with no consanguineous relationships among them, including adolescent young camels (YTR; n = 3; 2 years old) and lactating adult female camels (CNR; n = 3; 5 years old) from Fuhai County, Altay Prefecture, northwest China’s Xinjiang Uygur Autonomous Region (87° 35′25′′E, 46° 86′6′′N). The mammary glands of 2 years-old adolescent female camels are in the developmental stage and do not lactate. At 4 or 5 years of age, female camels begin to sexually mature until the end of their first pregnancy, when their mammary glands are fully developed and they have the ability to secrete milk after giving birth. All experimental camels were fed under the same conditions to reduce the environmental effects on gene expression.

A specialist veterinarian conducted sampling to minimize harm to the camels. The parenchyma of the mammary gland tissue, which contains lobules that synthesize milk (15), was harvested from camels by percutaneous needle biopsy. Part of the tissue samples were immediately frozen in liquid nitrogen then transferred to the laboratory and stored at −80°C for subsequent RNA extraction. The remaining mammary gland tissues were washed with sterile saline and preserved in a 2.5% glutaraldehyde solution for later use in hematoxylin-eosin (H&E) staining (Figure 1). All the experimental procedures were performed in accordance with the guidelines of the Laboratory Animal Administration Regulations issued by the National Science and Technology Committee (China). All experimental animal studies were reviewed and approved by the Animal Welfare Committee of Xinjiang University (approval ID: XJU2019012).

FIGURE 1

Figure 1. Experimental design of the study.

2.2. Micromorphological examination

The mammary gland tissues obtained from the experimental camels were fixed in 10% neutral buffered formalin, dehydrated in a graded alcohol series, cleared in xylol, embedded in paraffin, and sectioned to a thickness of 4–5 μm. Tissue sections were prepared and stained with H&E following the protocols of a previous study (16) and photographed under a light microscope (Eclipse E100 Nikon, Japan).

2.3. RNA extraction, quality, and integrity determination

Total RNA was extracted from the camel mammary tissue using a TRIZOL RNA Extraction Kit (Invitrogen, Carlsbad, CA, United States). RNA integrity was evaluated using agarose 1.0% gel electrophoresis and the RNA Nano 6,000 Assay Kit of the Agilent 2,100 Bioanalyzer (Agilent, Santa Clara, CA, United States). RNA samples that passed quality control were used for subsequent sequencing library construction.

2.4. Library construction and transcriptome sequencing

Library preparation and transcriptome sequencing were performed by Novogene Co., Ltd. (Tianjin, China). The mRNA was purified from 1 μg of total RNA using oligo dT magnetic beads. mRNA was fragmented using divalent cations at high temperatures. After fragmentation, the cDNA library was synthesized using random hexamer primers, M-MuLV Reverse Transcriptase, and paired-end sequencing using the HiSeq 6000 platform (Illumina) (17). Reads containing adapters, poly-N, and low-quality reads were removed to obtain clean reads by Fastp (v 0.19.7). The Q20, Q30, and GC contents of the clean data were calculated and all downstream analyses were performed using high-quality clean data. A set of genomic index files of Camelus bactrianus reference genome was built using Hisat2 (version 2.0.5) and clean paired-end reads were mapped to the reference genome using Hisat2. Gene expression levels were quantified using featureCounts (version 1.5.0-p3) and the fragments per kilobase of transcript per million mapped reads (FPKM) of each annotated gene were counted based on the length of the gene and the read counts mapped to this gene (18).

2.5. Differentially expressed genes identification

Principal component analysis (PCA) was used to visually identify differences between the different sequencing sample groups in the transcriptome data. Differential expression analysis was performed between the two groups of camels using the DESeq2 package in R (version 4.2.1), based on negative binomial distribution. Considering previous studies (19, 20) and our specific experimental situation, adjusted p-values (p_adjust) < 0.01 and |log2FoldChange| >2 were established as the thresholds for significance. p-values were adjusted for multiple testing using the Benjamini and Hochberg methods. Genes with p_adjust < 0.01 and log2FoldChange >2 by DESeq2 were defined as upregulated differentially expressed genes (DEGs); genes with p_adjust < 0.01 and log2FoldChange <−2 were regarded as downregulated DEGs.

2.6. Gene functional enrichment analysis of DEGs

To explore functions and pathways associated with the DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the NovoMagic Cloud Platform.¹ A cut-off of p < 0.05 was used to screen significant functions and pathways. We further screened the core pathways and candidate genes by systematically reviewing the literature.

2.7. Protein–protein interaction network of DEGs

Gene regulatory networks have elucidated the regulation of mammary gland development in animals (21, 22). Therefore, we uploaded all identified DEGs to the Search Tool for the Retrieval of Interacting Genes/Proteins database² to create a protein–protein interaction (PPI) network (23). PPI networks were established using Cytoscape software (version 4.8.0). Hub genes were identified and screened by network analysis using Cytoscape and its plugin (CytoHubba) (24).

2.8. Quantitative real-time polymerase chain reaction validation of DEGs

The Prime Script RT Reagent Kit (Takara Biotechnology Co., Ltd., Shiga, Japan) was used to produce cDNA from RNA. The obtained cDNA was analyzed immediately or stored at −20°C. Camel beta-actin (β-actin) (NCBI accession no. XM_010965866.2) was used as an internal reference gene (25). Fifteen genes were selected from the DEGs between the young camel and adult camel groups at another sampling site. Polymerase chain reaction (PCR) was performed using the following cycling conditions: 95°C for 5 min; 40 cycles of 95°C for 10 s and 60°C for 30 s; 95°C for 15 s, 60°C for 60 s, and 95°C for 15 s. Specific primers were designed using the Primer-BLAST online website. Three replicates of each sample were used to calculate the relative expression by the 2^−ΔΔct method (26). The identities of these DEGs and specific primer sequences are listed in Table 1.

TABLE 1

Table 1. The primers used for RT-qPCR are listed (F = forward, R = reverse).

3. Results

3.1. Histomorphometric analysis of mammary gland

We first examined the morphological and microscopic characteristics of camel mammary glands according to published literature on other species. The morphological differences at different stages of mammary gland development are shown in Figure 2. As reported for other young animals (27–29), the mammary gland tissue of young camels was incompletely developed and consisted of a large amount of connective tissue (Figures 2A,B). Furthermore, the terminal duct lobular units (TDLU) in the mammary gland were incompletely developed, with inner epithelial tissue surrounded by lobular connective tissue and an outer overlay of interlobular connective tissue, and with fewer and irregularly shaped ducts and epithelial cells. The mammary glands of adult female camels were more fully developed (Figures 2C,D), with a large number of mammary epithelial cells and uniform chromatin distribution, adipose tissue, a small amount of thin connective tissue, and fat pads scattered with alveoli, ducts, and alveolar cavities.

FIGURE 2

Figure 2. Mammary gland morphology in young and adult female camels. (A) Morphological structure of the mammary glands of young camels (200×). (B) Morphological structure of the mammary glands of young camels (400×). (C) Morphological structure of the mammary glands of adult camels (200×). (D) Morphological structure of the mammary glands of adult camels (400×).

3.2. RNA-sequencing reads and mapping to the reference genome

Raw reads were obtained for each library and between 43 and 59 million clean reads were obtained after quality control procedures and applied to the YTR and CNR libraries (Table 2). The Q30 reads (quality score >97.69%) were >95% for all samples (Supplementary Table S1). Between 38 and 53 million reads were mapped to the Bactrian camel genome for the YTR and CNR groups, respectively. Approximately 92% of the clean reads mapped to the camel reference genome, of which nearly 89% were uniquely mapped. The detailed statistical information is provided in Supplementary Table S2. These data provided a solid basis for subsequent analyses.

TABLE 2

Table 2. Summary of the RNA-seq analyses after mapping to the reference genome.

3.3. Analysis of DEGs

The number of genes at different expression intervals was determined (Supplementary Table S3). The first principal component (PCA1), which had the largest variance (41.43%), distinctly clustered the samples into two groups, with three biological replicates within each group clustered together (Figure 3A). A total of 2,851 genes were identified as differentially expressed: 1,420 (49.81%) were upregulated and 1,431 (50.19%) were downregulated (YTR vs. CNR) (Figure 3B and Supplementary Table S4). The expressed genes were subjected to hierarchical clustering to comprehensively and intuitively display the differences in gene expression between the two groups (Figure 3C). Statistical analysis of the gene biotypes was performed on the DEGs to identify the composition of gene categories. Of these, 2,419 (84.8%) were protein-coding genes, followed by unknown type genes (10.8%) and 93 (3.3%) long non-coding RNA genes (Figure 3D).

FIGURE 3

Figure 3. Analysis of differentially expressed genes among sample groups. (A) Principal-component analysis comparing the two sample groups. (B) Volcano plot of differentially expressed genes. (C) Heatmap and the hierarchical cluster analysis of the differentially expressed genes. (D) Biotype of genes statistical analysis.

3.4. Go enrichment analysis

To obtain a more comprehensive understanding of the DEGs, we performed GO functional enrichment analysis using DAVID. These DEGs were enriched in 827 GO terms, of which 434 were related to biological processes (BP), 279 to molecular functions (MF), and 114 to cellular components (CC). Moreover, collagen trimers and extracellular matrix structural constituents were significantly enriched (p_adjust < 0.05). The complete GO analysis results are provided in Supplementary Table S5. Significance analysis of GO enrichment of upregulated DEGs revealed 33 significantly enriched entries (p_adjust < 0.05), including 29 BP terms, 3 MF terms, and 1 CC term (Supplementary Table S6). These terms were related to cytoskeletal organization, organelle regulation, actin regulation, protein organization, cytoskeletal and component organization regulation, and calcium ion transport across membranes (Figure 4A).

FIGURE 4

Figure 4. Bubble diagram of Gene Ontology (GO) enrichment result. (A) GO enrichment data for genes that were up-regulated. (B) GO enrichment data for genes that were down-regulated.

GO enrichment of the downregulated DEGs revealed 13 significantly enriched entries (p_adjust < 0.05), including five BP terms, five MF terms, and three CC terms (Supplementary Table S7). These terms were involved in molecular functional regulators, riboprotein complexes, cytoplasmic fraction, peptidase regulator activity and inhibitor activity, ribosomes, translation, peptide metabolism and peptide biosynthesis processes, cellular amide metabolism processes, amide biosynthesis processes, ribosome structural components, and structural molecular activity (Figure 4B). In the reference GO terms reported in previous studies, these significantly enriched GO entries were intimately correlated with the maintenance of cell morphology and structure, organismal tissue growth and development, and response to cell growth. Among these, peptidase regulator activity and inhibitor activity entries can regulate aminopeptidase N to promote mammary gland development in animals.

3.5. KEGG pathway enrichment analysis

Signaling pathways are essential for mammary gland development. We performed KEGG pathway enrichment analysis to investigate the gene-enriched pathways. Overall, 2,851 DEGs were assigned to 318 KEGG pathways, of which 11 were significantly enriched (p_adjust < 0.05). The main pathways were insulin secretion, calcium and glucagon signaling pathway, steroid biosynthesis, HIF-1 signaling pathway, protein digestion and absorption, adrenergic signaling in cardiomyocytes, and tight junction (Supplementary Table S8).

The upregulated and downregulated genes were subjected to KEGG pathway enrichment analysis. The enrichment analysis results indicated that 24 pathways were significantly up-regulated and 7 were significantly down-regulated (p_adjust < 0.05) (Figure 5). Upregulated genes were significantly enriched in the Hedgehog, insulin, and oxytocin signaling pathways (Supplementary Table S9). Downregulated genes were significantly enriched in the Wingless-Type MMTV Integration Site Family (Wnt) signaling pathway, steroid biosynthesis, and tight junctions (Supplementary Table S10). As previously reported, the Hedgehog signaling pathway (30, 31), Wnt signaling pathway (32, 33), hormone-related pathways (34–36), and tight junctions (37) are closely associated with mammary gland development in animals. Collectively, these results provide a basis for selecting candidate genes related to mammary gland growth and development.

FIGURE 5

Figure 5. Dot plot of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis results.

3.6. PPI network of the DEGs

We explored the biological and regulatory functions of hub DEGs at the protein level. Next, we collated the genes in important pathways for the (PPI) network analysis and sorted the nodes according to the degree of interaction (Figure 6 and Supplementary Table S11). Among the top 10 central genes, PRKAB2 and PRKAG3 are associated with animal growth, slaughter, and meat quality traits (38); PLCB4 genes are important candidate genes for pig growth and development and average daily weight gain (39). BTRC affects spermatogenesis and mammary gland development in mice (40). The four hub genes showed a high degree of interaction. Notably, some non-hub genes, such as GLI1, WIF1, DKK2, FZD3, and WNT4 genes, have been reported to be associated with mammary gland development. This demonstrated that these genes play important roles in animal growth and development.

FIGURE 6

Figure 6. Networks analysis between protein–protein interaction (PPI) network and KEGG pathway related to mammary gland development. Non-hub genes are represented as grey circle while yellow octagons represent pathways. Hub genes are indicated by rectangles. The important candiate genes are highlighted in red. The light red solid line indicates the high degree of gene linkage. The light red dashed line indicates the low level of gene connectivity.

3.7. Validation of gene expression

For the preliminary validation of the RNA-sequencing data, fifteen DEGs were selected for experimental validation and statistical analysis in another random camel population, including nine upregulated genes (CAMK2D, CAMK2A, BTRC, PRKAB2, PLCB4, PPP1R3A, GLI1, DKK2, and KCNJ12) and six downregulated genes (SLC1A1, LTF, WIF1, FZD3, WNT4, and PIGR). The gene expression levels were normalized relative to β-actin. The expression trends of fifteen genes in the random camel population were consistent with the transcriptome results (Figure 7 and Supplementary Table S12).

FIGURE 7

Figure 7. Validation of RNA-seq differential gene expression levels using RT-qPCR. The relative expression levels of both the RT-qPCR and RNA-seq were expressed as log2FoldChange (FC), which is defined as the ratio of the expression level in the young camels to that in the adult camels for a specific gene. Error bars are standard deviations across three camels per group evaluated.

4. Discussion

Camels, like cattle, have four mammary glands in the inguinal region responsible for the synthesis and secretion of milk proteins (41). Because camels do not have gland cisterns for storage, milk production is a reflex secretion (42). Previous studies have demonstrated that milkability traits can be assessed using udder and teat morphological traits of camels, contributing to the selection of high-producing camels (43). Camel teats have a double ductal anatomy with a promising prophylactic effect against mastitis-causing pathogens (44). Although whole-genome sequencing of dromedaries and Bactrian camels has been completed (45, 46), gene annotation remains incomplete. It is necessary to obtain more data for the further annotation of camelids. To the best of our knowledge, the present study is the first RNA-sequencing study of mammary gland developmental gene expression in young and adult camels. Although the sample sizes in the current study were relatively limited, significant differences were found between the morphology and transcription profiles of mammary gland tissue of young and adult camels. Additionally, the three mammary gland samples in each group showed a high correlation, proving that the biological samples within each group had good repeatability.

The morphology of mammary glands undergoes significant changes during development. A distinctive feature of the mammary gland development is its extensive proliferative and differentiation potential, particularly during puberty and pregnancy. Adolescence is the ideal stage for studying mammary gland development. The mammary gland develops its adult form through a biological process referred to as branching morphogenesis (47). Our observations of the mammary gland are consistent with those of previous studies (44), showing that age and lactation strongly influence the microscopic anatomy of the camel mammary gland.

Extensive experimental efforts have been made to investigate the endocrine signaling and pathways that control the proliferation and differentiation of mammary epithelial cells. Such studies have reported the regulatory roles of multiple signaling pathways in mammary gland proliferation and development. Several studies have revealed that the Wnt pathway is required at the earliest stage of mammary development in the embryo for specification of the mammary placode and initiation of mammary morphogenesis (48) and in postnatal development for proper stem cell maintenance, branching morphogenesis, and mammary alveolar development (49, 50). The Hedgehog signaling pathway is also necessary for mammary gland development and distinguishes the mammary gland from other epidermal appendages (51, 52). Mammary gland development involves complex interactions among various hormones. The mammary gland is composed of a branched ductal system consisting of milk-producing epithelial cells that form ductile tubules surrounded by a layer of myoepithelial cells that contract in response to oxytocin stimulation, thereby releasing milk (53). After the embryonic and prepubertal stages, further development of the mammary glands is highly dependent on hormones such as steroids, which control the development of breast ducts and alveoli (53). The results of RNA sequencing and histomorphological identification of the mammary glands of young and adult camels lend further support to these pathways and hormones regulating mammary gland development in animals.

The functions of genes involved in these pathways require further investigation. Comprehensive PPI network analysis identified ten hub genes. Among these, variations in the BTRC gene have been reported to influence spermatogenesis and mammary gland development in mice (40). PRKAB2 and PRKAG3 have been implicated in animal growth, slaughter, and meat quality regulation (38). In pigs, PLCB4 is a candidate gene for growth, development, and meat production (39). Among these non-core genes, GLI1 regulates the glioma-associated oncogene homolog 1 (GLI1) activator and causes mammary bud formation failure in mice (54). The WIF1 gene is expressed at high levels in normal human breast cells and mammary tissues and maintains mammary gland development (55). DKK2 and FZD3 are closely related to the Wnt signaling pathway during embryonic mammary gland development (56). The Wnt ligand WNT4 plays a critical role in normal mammary gland development. Our study identified a set of crucial genes that might exert their functions by regulating the Wnt, Hedgehog, and sex hormone signaling pathways. However, additional functional experiments are required to confirm this hypothesis.

In summary, this study is of great value for identifying critical genes involved in mammary gland development and can inform subsequent studies on the progression mechanism of mammary gland development in camels. However, there were some limitations. The present study is limited by the small samples size of each group. Three camels were used as replicates at each stage. The findings of the study need to be verified in the future through studies with an expanded sample size. Additionally, there are differences in genetic information between individual animals, which may affect the results of the differential expression analysis of the mammary glands (57). Although distinct differences in the function and structure of mammary glands have been found between young and adult camels, mammary gland growth and development are complex biological processes, and inclusion of more sampling and sequencing time points might help to elucidate the molecular change mechanisms of the underlying biological processes and better reflect the continuous characteristics of dynamic gene changes during mammary gland development.

5. Conclusion

The Bactrian camel is a versatile animal with high ecological and economic value in remote desert areas. Milk synthesis is an essential function of the mammary glands in mammals. This study described the microstructure and transcriptome profiles of the mammary gland parenchyma of young and adult camels. Differential expression analysis revealed 2,851 DEGs, of which 1,420 were upregulated and 1,431 were downregulated. Functional enrichment analysis revealed significant enrichment of the Wnt, Hedgehog, oxytocin, and insulin signaling pathways and steroid biosynthesis. The PRKAB2, PRKAG3, PLCB4, BTRC, GLI1, WIF1, DKK2, FZD3, and WNT4 genes involved in these pathways should be considered as important candidate genes for mammary gland development. Nevertheless, the potential functions of these genes still require further verification. This study revealed DEGs and signaling pathways related to mammary gland development in Bactrian camels and laid a theoretical foundation for the improvement of milk production traits in camels.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/bioproject; PRJNA946168.

Ethics statement

The animal study was reviewed and approved by Ethics Committee of Xinjiang University. Written informed consent was obtained from the owners for the participation of their animals in this study.

Author contributions

JY conceived the study. H Yao writing-original draft and writing-review and editing. H Yao, ZD, and XL conducted conceptualization and writing-review and editing. ZZ, WM, ZH, H Yan, YW, ZH, ZW, and GC collected the experimental samples. All authors contributed to the article and approved the submitted version.

Funding

This work was financially supported by the Key Technology Research and Development Program in Xinjiang Uygur Autonomous Region (2018B01003), the National Key Research and Development Projects of China (2019YFC1606103), and the Postgraduate Scientific Research Innovation Program of Xinjiang Uygur Autonomous Region (XJ2019G026).

Acknowledgments

The authors would like to extend their special thanks to the camel breeders and camel farm managers, as well as to the experts, veterinarians, and laboratories who lent their assistance.

Conflict of interest

GC was employed by Xinjiang Wangyuan Camel Milk Limited Company.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2023.1196950/full#supplementary-material

Footnotes

1. ^https://magic.novogene.com/customer/main#/loginNew

2. ^http://string-db.org/

References

1. Faye, B, and Konuspayeva, G. The encounter between bactrian and dromedary camels in central Asia, vol. 451 Verlag der Österreichischen Akademie der Wissenschaften (2012) 27–33.

Google Scholar

2. Yang, J, Dou, Z, Peng, X, Wang, H, Shen, T, Liu, J, et al. Transcriptomics and proteomics analyses of anti-cancer mechanisms of TR35–an active fraction from Xinjiang Bactrian camel milk in esophageal carcinoma cell. Clin Nutr. (2019) 38:2349–59. doi: 10.1016/j.clnu.2018.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Abdalla, E, Anis Ashmawy, AEH, Farouk, MH, Abd El-Rahman Salama, O, Khalil, F, and Seioudy, A. Milk production potential in Maghrebi she-camels. Small Rumin Res. (2015) 123:129–35. doi: 10.1016/j.smallrumres.2014.11.004

CrossRef Full Text | Google Scholar

4. Sun, L, Fu, Y, Yang, Y, Wang, X, Cui, W, Li, D, et al. Genomic analyses reveal evidence of independent evolution, demographic history, and extreme environment adaptation of Tibetan plateau Agaricus bisporus. Front Microbiol. (2019) 10:1786. doi: 10.3389/fmicb.2019.01786

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Xuan, R, Chao, T, Zhao, X, Wang, A, Chu, Y, Li, Q, et al. Transcriptome profiling of the nonlactating mammary glands of dairy goats reveals the molecular genetic mechanism of mammary cell remodeling. J Dairy Sci. (2022) 105:5238–60. doi: 10.3168/jds.2021-21039

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bai, X, Zheng, Z, Liu, B, Ji, X, Bai, Y, and Zhang, W. Whole blood transcriptional profiling comparison between different milk yield of Chinese Holstein cows using RNA-seq data. BMC Genomics. (2016) 17:512. doi: 10.1186/s12864-016-2901-1

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Harris, HR, Willett, WC, Vaidya, RL, and Michels, KB. Adolescent dietary patterns and premenopausal breast cancer incidence. Carcinogenesis. (2016) 37:376–84. doi: 10.1093/carcin/bgw023

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Arora, R, Siddaraju, NK, Manjunatha, SS, Sudarshan, S, Fairoze, MN, Kumar, A, et al. Muscle transcriptome provides the first insight into the dynamics of gene expression with progression of age in sheep. Sci Rep. (2021) 11:22360. doi: 10.1038/s41598-021-01848-5

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Lian, Y, Gòdia, M, Castelló, A, Rodríguez-Gil, J, Maestre, S, Sanchez, A, et al. Characterization of the impact of density gradient centrifugation on the profile of the pig sperm transcriptome by RNA-seq. Front Vet Sci. (2021) 8:8. doi: 10.3389/fvets.2021.668158

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tucker, HLM, Parsons, CLM, Ellis, S, Rhoads, ML, and Akers, RM. Tamoxifen impairs prepubertal mammary development and alters expression of estrogen receptor α (ESR1) and progesterone receptors (PGR). Domest Anim Endocrinol. (2016) 54:95–105. doi: 10.1016/j.domaniend.2015.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Velayudhan, BT, Huderson, BP, McGilliard, ML, Jiang, H, Ellis, SE, and Akers, RM. Effect of staged ovariectomy on measures of mammary growth and development in prepubertal dairy heifers. Animal. (2012) 6:941–51. doi: 10.1017/s1751731111002333

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Nørgaard, JV, Nielsen, MO, Theil, PK, Sørensen, MT, Safayi, S, and Sejrsen, K. Development of mammary glands of fat sheep submitted to restricted feeding during late pregnancy. Small Rumin Res. (2008) 76:155–65. doi: 10.1016/j.smallrumres.2007.11.001

CrossRef Full Text | Google Scholar

13. Hughes, K. Comparative mammary gland postnatal development and tumourigenesis in the sheep, cow, cat and rabbit: exploring the menagerie. Semin Cell Dev Biol. (2020) 114:186–95. doi: 10.1016/j.semcdb.2020.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Williams, SA, Harata-Lee, Y, Comerford, I, Anderson, RL, Smyth, MJ, and McColl, SR. Multiple functions of CXCL12 in a syngeneic model of breast cancer. Mol Cancer. (2010) 9:250. doi: 10.1186/1476-4598-9-250

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hao, Z, Zhou, H, Hickford, JGH, Gong, H, Wang, J, Hu, J, et al. Transcriptome profile analysis of mammary gland tissue from two breeds of lactating sheep. Genes. (2019) 10:781. doi: 10.3390/genes10100781

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yao, H, Liu, M, Ma, W, Yue, H, Su, Z, Song, R, et al. Prevalence and pathology of Cephalopina titillator infestation in Camelus bactrianus from Xinjiang, China. BMC Vet Res. (2022) 18:360. doi: 10.1186/s12917-022-03464-5

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhu, H, Blum, RH, Bernareggi, D, Ask, EH, Wu, Z, Hoel, HJ, et al. Metabolic reprograming via deletion of CISH in human iPSC-derived NK cells promotes in vivo persistence and enhances anti-tumor activity. Cell Stem Cell. (2020) 27:224–237.e6. doi: 10.1016/j.stem.2020.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lan, Y, Sun, J, Chen, C, Sun, Y, Zhou, Y, Yang, Y, et al. Hologenome analysis reveals dual symbiosis in the deep-sea hydrothermal vent snail Gigantopelta aegis. Nat Commun. (2021) 12:1165. doi: 10.1038/s41467-021-21450-7

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu, G, Yang, G, Zhao, G, Guo, C, Zeng, Y, Xue, Y, et al. Spatial transcriptomic profiling to identify mesoderm progenitors with precision genomic screening and functional confirmation. Cell Prolif. (2022) 55:e13298. doi: 10.1111/cpr.13298

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Peng, H, Guo, Q, Xiao, Y, Su, T, Jiang, TJ, Guo, LJ, et al. ASPH regulates osteogenic differentiation and cellular senescence of BMSCs. Front Cell Dev Biol. (2020) 8:872. doi: 10.3389/fcell.2020.00872

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Rudolph, MC, McManaman, JL, Hunter, L, Phang, T, and Neville, MC. Functional development of the mammary gland: use of expression profiling and trajectory clustering to reveal changes in gene expression during pregnancy, lactation, and involution. J Mammary Gland Biol Neoplasia. (2003) 8:287–307. doi: 10.1023/b:jomg.0000010030.73983.57

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Loor, JJ, Batistel, F, Bionaz, M, Hurley, WL, and Vargas-Bello-Pérez, E. Mammary gland: gene networks controlling development and involution^☆ In: PLH McSweeney and JP McNamara, editors. Encyclopedia of dairy sciences. 3rd ed. Oxford: Academic Press (2022). 167–74.

Google Scholar

23. Szklarczyk, D, Gable, AL, Lyon, D, Junge, A, Wyder, S, Huerta-Cepas, J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. (2019) 47:D607–d613. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Nangraj, AS, Selvaraj, G, Kaliamurthi, S, Kaushik, A, Cho, W, and Wei, D-Q. Integrated PPI and WGCNA retrieving shared gene signatures between Barrett’s esophagus and esophageal adenocarcinoma. Front Pharmacol. (2020) 11:11. doi: 10.3389/fphar.2020.00881

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yu, X, Wu, Y, Zhang, J, Jirimutu, ZA, and Chen, J. Pre-evaluation of humoral immune response of Bactrian camels by the quantification of Th2 cytokines using real-time PCR. J Biomed Res. (2020) 34:387–94. doi: 10.7555/jbr.34.20190035

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Liu, Z, An, L, Lin, S, Wu, T, Li, X, Tu, J, et al. Comparative physiological and transcriptomic analysis of pear leaves under distinct training systems. Sci Rep. (2020) 10:18892. doi: 10.1038/s41598-020-75794-z

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Hindman, AR, Mo, XM, Helber, HL, Kovalchin, CE, Ravichandran, N, Murphy, AR, et al. Varying susceptibility of the female mammary gland to in utero windows of BPA exposure. Endocrinology. (2017) 158:3435–47. doi: 10.1210/en.2017-00116

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Finot, L, Chanat, E, and Dessauge, F. Molecular signature of the putative stem/progenitor cells committed to the development of the bovine mammary gland at puberty. Sci Rep. (2018) 8:16194. doi: 10.1038/s41598-018-34691-2

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Chatterjee, SJ, Halaoui, R, Deagle, RC, Rejon, C, and McCaffrey, L. Numb regulates cell tension required for mammary duct elongation. Biol Open. (2019) 8:bio042341. doi: 10.1242/bio.042341

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Gerashchenko, TS, Zolotaryova, SY, Kiselev, AM, Tashireva, LA, Novikov, NM, Krakhmal, NV, et al. The activity of KIF14, Mieap, and EZR in a new type of the invasive component, Torpedo-like structures, predetermines the metastatic potential of breast Cancer. Cancers. (2020) 12:1909. doi: 10.3390/cancers12071909

CrossRef Full Text | Google Scholar

31. Xing, Z, Lin, A, Li, C, Liang, K, Wang, S, Liu, Y, et al. lncRNA directs cooperative epigenetic regulation downstream of chemokine signals. Cells. (2014) 159:1110–25. doi: 10.1016/j.cell.2014.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Shahzad, N, Munir, T, Javed, M, Tasneem, F, Aslam, B, Ali, M, et al. SHISA3, an antagonist of the Wnt/β-catenin signaling, is epigenetically silenced and its ectopic expression suppresses growth in breast cancer. PLoS One. (2020) 15:e0236192. doi: 10.1371/journal.pone.0236192

CrossRef Full Text | Google Scholar

33. Jia, S, Zhou, J, Wee, Y, Mikkola, ML, Schneider, P, and D’Souza, RN. Anti-EDAR agonist antibody therapy resolves palate defects in Pax9^−/− mice. J Dent Res. (2017) 96:1282–9. doi: 10.1177/0022034517726073

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Sun, G, Wang, C, Wang, S, Sun, H, Zeng, K, Zou, R, et al. An H3K4me3 reader, BAP18 as an adaptor of COMPASS-like core subunits co-activates ERα action and associates with the sensitivity of antiestrogen in breast cancer. Nucleic Acids Res. (2020) 48:10768–84. doi: 10.1093/nar/gkaa787

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Reynolds, P, Canchola, AJ, Duffy, CN, Hurley, S, Neuhausen, SL, Horn-Ross, PL, et al. Urinary cadmium and timing of menarche and pubertal development in girls. Environ Res. (2020) 183:109224. doi: 10.1016/j.envres.2020.109224

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang, Y, Gc, S, Patel, SB, Liu, Y, Paterson, AJ, Kappes, JC, et al. Growth hormone (GH) receptor (GHR)-specific inhibition of GH-induced signaling by soluble IGF-1 receptor (sol IGF-1R). Mol Cell Endocrinol. (2019) 492:110445. doi: 10.1016/j.mce.2019.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Dianati, E, Poiraud, J, Weber-Ouellette, A, and Plante, I. Connexins, E-cadherin, claudin-7 and β-catenin transiently form junctional nexuses during the post-natal mammary gland development. Dev Biol. (2016) 416:52–68. doi: 10.1016/j.ydbio.2016.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Milan, D, Jeon, J-T, Looft, C, Amarger, V, Robic, A, Thelander, M, et al. A mutation in PRKAG3 associated with excess glycogen content in pig skeletal muscle. Science. (2000) 288:1248–51. doi: 10.1126/science.288.5469.1248

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Edea, Z, Hong, JK, Jung, JH, Kim, DW, Kim, YM, Kim, ES, et al. Detecting selection signatures between Duroc and Duroc synthetic pig populations using high-density SNP chip. Anim Genet. (2017) 48:473–7. doi: 10.1111/age.12559

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cortellari, M, Barbato, M, Talenti, A, Bionda, A, Carta, A, Ciampolini, R, et al. The climatic and genetic heritage of Italian goat breeds with genomic SNP data. Sci Rep. (2021) 11:10986. doi: 10.1038/s41598-021-89900-2

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Yang, Y, Fang, X, Yang, R, Yu, H, Jiang, P, Sun, B, et al. MiR-152 regulates apoptosis and triglyceride production in MECs via targeting ACAA2 and HSD17B12 genes. Sci Rep. (2018) 8:417. doi: 10.1038/s41598-017-18804-x

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Alluwaimi, A. The Camel’s (Camelus dromedarius) mammary gland immune system in health and disease. Adv Dairy Res. (2017) 05:5. doi: 10.4172/2329-888X.1000171

CrossRef Full Text | Google Scholar

43. Kumar, M, Nehara, M, Prakash, V, Pannu, U, and Jyotsana, B. Udder, teat, and milk vein measurements of Indian dromedary camel and its relationship with milkability traits. Trop Anim Health Prod. (2023) 55:36. doi: 10.1007/s11250-023-03457-y

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Kausar, R. Gross and microscopic anatomy of mammary gland of dromedaries under different physiological conditions. Pak Vet J. (2001) 21:189–93.

Google Scholar

45. Jirimutu, WZ, Ding, G, Chen, G, Sun, Y, Sun, Z, Zhang, H, et al. The Bactrian camels genome sequencing and analysis consortium. correction: corrigendum: genome sequences of wild and domestic bactrian camels. Nat Commun. (2013) 4:2089. doi: 10.1038/ncomms3089

CrossRef Full Text | Google Scholar

46. Fitak, RR, Mohandesan, E, Corander, J, and Burger, PA. The de novo genome assembly and annotation of a female domestic dromedary of North African origin. Mol Ecol Resour. (2016) 16:314–24. doi: 10.1111/1755-0998.12443

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Sieh, W, Rothstein, JH, Klein, RJ, Alexeeff, SE, Sakoda, LC, Jorgenson, E, et al. Identification of 31 loci for mammographic density phenotypes and their associations with breast cancer risk. Nat Commun. (2020) 11:5116. doi: 10.1038/s41467-020-18883-x

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Roarty, K, Shore, AN, Creighton, CJ, and Rosen, JM. Ror2 regulates branching, differentiation, and actin-cytoskeletal dynamics within the mammary epithelium. J Cell Biol. (2015) 208:351–66. doi: 10.1083/jcb.201408058

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Watson, CJ, and Khaled, WT. Mammary development in the embryo and adult: a journey of morphogenesis and commitment. Development. (2008) 135:995–1003. doi: 10.1242/dev.005439

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Szmatoła, T, Gurgul, A, Jasielczuk, I, Ząbek, T, Ropka-Molik, K, Litwińczuk, Z, et al. A comprehensive analysis of runs of homozygosity of eleven cattle breeds representing different production types. Animals. (2019) 9:1024. doi: 10.3390/ani9121024

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Hatsell, SJ, and Cowin, P. Gli3-mediated repression of hedgehog targets is required for normal mammary development. Development. (2006) 133:3661–70. doi: 10.1242/dev.02542

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Behbod, F, and Rosen, JM. Will cancer stem cells provide new therapeutic targets? Carcinogenesis. (2005) 26:703–11. doi: 10.1093/carcin/bgh293

CrossRef Full Text | Google Scholar

53. Hennighausen, L, and Robinson, GW. Signaling pathways in mammary gland development. Dev Cell. (2001) 1:467–75. doi: 10.1016/S1534-5807(01)00064-8

CrossRef Full Text | Google Scholar

54. Niyaz, M, Khan, MS, and Mudassar, S. Hedgehog signaling: An Achilles’ heel in Cancer. Transl Oncol. (2019) 12:1334–44. doi: 10.1016/j.tranon.2019.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Wang, N, Wang, Z, Wang, Y, Xie, X, Shen, J, Peng, C, et al. Dietary compound isoliquiritigenin prevents mammary carcinogenesis by inhibiting breast cancer stem cells through WIF1 demethylation. Oncotarget. (2015) 6:9854–76. doi: 10.18632/oncotarget.3396

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Pausch, H, Jung, S, Edel, C, Emmerling, R, Krogmeier, D, Götz, KU, et al. Genome-wide association study uncovers four QTL predisposing to supernumerary teats in cattle. Anim Genet. (2012) 43:689–95. doi: 10.1111/j.1365-2052.2012.02340.x

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Conesa, A, Madrigal, P, Tarazona, S, Gomez-Cabrero, D, Cervera, A, McPherson, A, et al. A survey of best practices for RNA-seq data analysis. Genome Biol. (2016) 17:13. doi: 10.1186/s13059-016-0881-8

PubMed Abstract | CrossRef Full Text | Google Scholar