Systematic and Integrated Analysis Approach to Prioritize Mastitis Resistant Genes

Iqra Mahmood1, Asif Nadeem1*, Masroor Ellahi Babar2, Muhammad Muddassir Ali1, Maryam Javed1, Aisha Siddiqa1, Tanveer Hussain2 and Muhammad Tariq Pervez2

1Institute of Biochemistry and Biotechnology, University of Veterinary and Animal Sciences, Lahore 54000, Pakistan

2Virtual University of Pakistan, Lahore 54000, Pakistan

ABSTRACT

Mastitis (mammary gland inflammation) is one of the most important diseases causing economic losses to dairy farm owners. In order to develop new strategies to prevent mastitis, a detailed understanding of the molecular mechanisms underlying the host immune response to the infectious agents is necessary. Identification of the genes, and their variants, involved in innate immune responses is essential for the understanding of this inflammatory disease and to identify potential genetic markers for resistance to mastitis. This article presents a systematic integration of complex biological interactions of 226 mammary gland genes to uncover underlying regulatory networks. Data were collected from various databases after having an exhaustive literature study. Network, functional and pathway analysis of the mammary gland genes were performed to find mastitis resistant genes, which helped us to identify prime genes of mammary gland that are more resistant to mastitis, and analysis of those resistant genes helped to find the solution of mastitis and immunity related diseases.

Article Information

Received 03 August 2016

Revised 28 August 2016

Accepted 07 September

Available online 02 Novemebr 2016

Authors’ Contributions

AN and MMA conceived the study; IM, AN, MEB, MMA, MJ, AS, TH and MTP designed and performed the experiments. AN, MJ, IM drafted the article.

Key words

Systematics, Inflammatory diseases, Intra-mammary infection, Mammary glands, Immunity, Bovine genome, Networks, Pathways, Mastitis resistance genes

DOI: http://dx.doi.org/10.17582/journal.pjz/2017.49.1.101.106

* Corresponding author: asifnadeem@uvas.edu.pk

0030-9923/2016/0004-1161 $ 8.00/0

Copyright 2016 Zoological Society of Pakistan

Introduction

Worldwide human consumption of bovine milk and by- products has increased in the last decades, especially in countries where this food was not part of the traditional diet (Wiley, 2007). The bovine mammary gland is an extraordinary organ able to produce 30,000 kg of milk in a complete lactation cycle (Bionaz et al., 2012). Major advances have been seen in the past several decades in understanding the physiology of lactating mammary gland driven partly by the desire to improve the efficiency of production and partly by the desire to increase quality of the milk (Bauman et al., 2006). The development of high throughput tools such as microarray analysis has provided the unique feature to uncover the molecular networks in mammary tissue during the course of pregnancy, milk synthesis, involution and diseases related to mammary gland e.g., mastitis (Rudolph et al., 2003). However, interpretation of this large amount of data was extremely challenging. The advent of bioinformatics tools has allowed summarizing large microarray data sets, revealing an integral view of cellular phenomena (Lemay et al., 2007). Mastitis is an inflammatory disease of the bovine mammary gland that occurs in response to physical damage or infection with pathogenic microorganisms. Infection of the mammary gland causes macrophages and epithelial cells to release cytokines which cause polymorpho-nuclear neutrophils (PMNs), monocytes, and other leukocytes to migrate from the blood to the site of infection in the mammary tissue (Schukken et al., 2011; Paape et al., 2003). This influx of leukocytes results in an increased level of somatic cells in the milk. It costs more than $2 billion alone in USA (Viguier et al., 2009). In Pakistan only in Punjab Rs. 240 million is the average loss due to mastitis per year (Khan and Khan, 2006). It usually occurs primarily in response to intra-mammary bacterial infection but mycoplasmal, fungal, or algal infections are also causes of mastitis among which Staphylococcus aureus and Escherichia coli are the most important gram-positive and gram-negative bacteria, respectively (Schukken et al., 2011). The prevalence of mastitis in dairy cattle is relatively high.

Keeping in view the importance of milk products and loss due to mastitis, there is a need to develop new and improved/control strategies. For this purpose a detailed understanding of the molecular mechanisms underlying the host immune response to the infectious agents is necessary (Jensen et al., 2013). Although the differential expression of genes studied through microarray technique has helped to understand the gene networks and functional pathways affected most in mammary tissue in response to an intra-mammary infection with S. uberis (Moyes et al., 2009), identification of the genes and their variants involved in innate immune responses is essential for the understanding of this inflammatory disease and to identify potential genetic markers for resistance to mastitis in bovine. This article presents analytical study of bovine mammary gland genes to uncover underlying regulatory networks.

MATERIALS AND METHODS

Literature survey and date collection

This article presents analytical study of bovine mammary gland genes to uncover underlying regulatory networks. Network, functional and pathway analysis of 226 mammary gland genes were performed to find mastitis resistant genes. Data was collected from UniGene, NCBI, Ensembl and bovine genome databases. Network analysis was performed by STRING (Search Tool for the Retrieval of Interacting Genes/Proteins). For functional analysis DAVID (the Database for Annotation, Visualization and Integrated Discovery) was used and pathway analysis was performed through KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa, 2013).

A list of genes representing part of the mammalian mammary gland transcriptome was assembled by an extensive literature survey and searching various databases. The UniGene which is an NCBI database of the transcriptome was searched to identify genes expressed in bovine mammary gland. Literature published was reviewed by searching for the relevant publications through PubMed using key phrases: association, gene candidates, epigenetics, genetics, mammary gland, mastitis, methylation, milk, miRNA, QTL, SNP. A systematic review of both original research articles and reviews was performed for searching literature associated with bovine mammary gland genes. Data of genes was retrieved from the NCBI, Ensembl and bovine genome databases.

Network analysis

Complete knowledge of all direct and indirect interactions between proteins in a given cell would represent an important milestone towards a comprehensive description of cellular mechanisms and functions. Network analysis of the list of mammary gland genes was completed by STRING. It sorted out our data according to predicted functional association among different genes and proteins. We used four views for our analysis which were, Co-expression, database, experiment and text mining at high confidence score of 0.700. This analysis helped in systematically making the group of genes performing similar functions and showing similar expression pattern.

Functional analysis

DAVID was used to perform analysis of functionally related genes. Functional annotation and gene functional classification tools were used for the functional analysis of mammary gland genes at p <0.05. This gave us an in-depth understanding of the biological themes in list of genes that are enriched in genome-scale studies.

Pathway analysis

Pathway analysis of mastitis resistant genes was performed through KEGG which gave an insight into pathways where the genes were involved and it helped to interpret the data in the context of biological processes, pathways and networks. This analysis helped to understand the mechanisms of mastitis, the unaffected pathways and the pathways that were are affected most. These findings ultimately helped to find resistance mechanism and genes involved in resistance.

Fig. 1. The network showing 3 clusters: In the center there is the largest network of the dataset; in the middle there is a second largest network of genes and at the sides scattered genes are shown that are not interconnected to any other gene in the provided dataset.

Table I.- Gene functional classification.

Gene Group 1 | Enrichment Score: 14.712468824680236
OFFICIAL_GENE_SYMBOL | Gene Name
CCL11 | chemokine (C-C motif) ligand 11
IL13 | interleukin 13
IL8 | interleukin 8
CCL2 | chemokine (C-C motif) ligand 2
CCL3 | chemokine (C-C motif) ligand 3
CCL5 | chemokine (C-C motif) ligand 5
CXCL16 | chemokine (C-X-C motif) ligand 16
CXCL5 | chemokine (C-X-C motif) ligand 5
CSF3 | colony stimulating factor 3 (granulocyte)
CCL4 | chemokine (C-C motif) ligand 4
IL1RN | interleukin 1 receptor antagonist
CCL20 | chemokine (C-C motif) ligand 20
CXCL10 | chemokine (C-X-C motif) ligand 10
Gene Group 2 | Enrichment Score: 6.197470798536051
OFFICIAL_GENE_SYMBOL | Gene Name
IFNG | interferon, gamma
IL10 | interleukin 10
IL1B | interleukin 1, beta
IL2 | interleukin 2
TNF | tumor necrosis factor (TNF superfamily, member 2)
Gene Group 3 | Enrichment Score: 3.3291366663923707
OFFICIAL_GENE_SYMBOL | Gene Name
RELA | v-relreticuloendotheliosis viral oncogene homolog A (avian); similar to v-relreticuloendotheliosis viral oncogene homolog A, nuclear factor of kappa light polypeptide gene enhancer in B-cells 3, p65
NFKB2 | nuclear factor of kappa light polypeptide gene enhancer in B-cells 2 (p49/p100)
NFKB1 | nuclear factor of kappa light polypeptide gene enhancer in B-cells 1
RELB | v-relreticuloendotheliosis viral oncogene homolog B
Gene Group 4 | Enrichment Score: 2.1814765223414283
OFFICIAL_GENE_SYMBOL | Gene Name
IGF1R | insulin-like growth factor 1 receptor
YES1 | v-yes-1 Yamaguchi sarcoma viral oncogene homolog 1
FGFR1 | fibroblast growth factor receptor 1
MAP3K8 | mitogen-activated protein kinase kinasekinase 8
MET | met proto-oncogene (hepatocyte growth factor receptor)
HCK | hemopoietic cell kinase
KIT | v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog
Gene Group 5 | Enrichment Score: 1.8117924213734606
OFFICIAL_GENE_SYMBOL | Gene Name
S100A2 | S100 calcium binding protein A2
S100A8 | S100 calcium binding protein A8
S100A9 | S100 calcium binding protein A9
S100A12 | S100 calcium binding protein A12 (calgranulin C)
Gene Group 6 | Enrichment Score: 1.6666772544237716
OFFICIAL_GENE_SYMBOL | Gene Name
CX3CR1 | chemokine (C-X3-C motif) receptor 1
C5AR1 | complement component 5a receptor 1
CCR2 | chemokine (C-C motif) receptor 2
CCR1 | chemokine (C-C motif) receptor 1

Fig. 2. Co-expressed genes in the dataset. The paler red color represent the less co-expressed genes, the dark red color shows the genes are more co-expressed.

Table II.- KEGG pathway annotation.

Category | Term
KEGG_PATHWAY | bta04060:Cytokine-cytokine receptor interaction
KEGG_PATHWAY | bta04620:Toll-like receptor signaling pathway
KEGG_PATHWAY | bta04621:NOD-like receptor signaling pathway
KEGG_PATHWAY | bta04660:T cell receptor signaling pathway
KEGG_PATHWAY | bta05200:Pathways in cancer
KEGG_PATHWAY | bta04010:MAPK signaling pathway
KEGG_PATHWAY | bta04062:Chemokine signaling pathway
KEGG_PATHWAY | bta04510:Focal adhesion
KEGG_PATHWAY | bta04640:Hematopoietic cell lineage
KEGG_PATHWAY | bta04210:Apoptosis
KEGG_PATHWAY | bta05222:Small cell lung cancer
KEGG_PATHWAY | bta04920:Adipocytokine signaling pathway
KEGG_PATHWAY | bta05215:Prostate cancer
KEGG_PATHWAY | bta05218:Melanoma
KEGG_PATHWAY | bta04623:Cytosolic DNA-sensing pathway
KEGG_PATHWAY | bta05330:Allograft rejection
KEGG_PATHWAY | bta04940:Type I diabetes mellitus
KEGG_PATHWAY | bta04622:RIG-I-like receptor signaling pathway
KEGG_PATHWAY | bta04630:Jak-STAT signaling pathway
KEGG_PATHWAY | bta04520:Adherens junction
KEGG_PATHWAY | bta05332:Graft-versus-host disease
KEGG_PATHWAY | bta05020:Prion diseases
KEGG_PATHWAY | bta05210:Colorectal cancer
KEGG_PATHWAY | bta05310:Asthma
KEGG_PATHWAY | bta04662:B cell receptor signaling pathway
KEGG_PATHWAY | bta04722:Neurotrophin signaling pathway
KEGG_PATHWAY | bta05214:Glioma
KEGG_PATHWAY | bta04810:Regulation of actin cytoskeleton
KEGG_PATHWAY | bta05219:Bladder cancer
KEGG_PATHWAY | bta04610:Complement and coagulation cascades
KEGG_PATHWAY | bta04650:Natural killer cell mediated cytotoxicity
KEGG_PATHWAY | bta05220:Chronic myeloid leukemia
KEGG_PATHWAY | bta04512:ECM-receptor interaction
KEGG_PATHWAY | bta05221:Acute myeloid leukemia
KEGG_PATHWAY | bta00590:Arachidonic acid metabolism

RESULTS

Network analysis

A list of total 226 mammary gland genes was created by systematic and integrated review of literature and provided to STRING. Two hundred and five genes matched with STRING database for Bos Taurus and a network was generated of these 205 genes at high confidence score of 0.700 (Fig. 1). One hundred and thirty four genes form the largest network. Figure 2 shows the genes that are co-expressed.

Functional analysis

DAVID was used for functional analysis of the short listed genes by network analysis. 134 genes that were the part of the largest network were selected for functional analysis to uncover the meaning of interaction of these genes. Functional annotation clustering resulted in a total of 72 clusters. The majority of gene ontology (GO) terms were related to immune response, for example, GO identifiers were related to inflammatory response (GO:0006954), chemotaxis (GO:0006935), immune response (GO:0006955), leukocyte migration (GO:0050900), response to lipopolysaccharide (GO:0032496), and TLR signalling pathway (GO:0002224). To get significant results, we applied an enrichment score cutoff criterion of >2. The gene list was further sorted out by using another DAVID tool of gene functional classification which sorted the gene list into 6 gene groups as shown in Table I.

Pathway analysis

Pathway analysis of the functionally related genes was conducted by using KEGG database. Thirty five significant pathways related to these genes were identified among which 11 pathways were directly related to the immunity related pathways. Pathways other than immunity related pathways were apoptosis, small cell lung cancer, adipocytokine signaling pathway, prostate cancer, melanoma, cytosolic DNA-sensing pathway, allograft rejection, type I diabetes mellitus, adherens junction, graft-versus-host disease, prion diseases, colorectal cancer, asthma, glioma, regulation of actin cytoskeleton, bladder cancer, complement and coagulation cascades, chronic myeloid leukemia, ECM-receptor interaction, acute myeloid leukemia, arachidonic acid metabolism. Total pathways associated with the provided gene list are represented in Table II.

DISCUSSION

The biological complexity of animals unavoidably requires a systems biology approach, i.e., a way to systematically study the complex biological interactions using a method of integration instead of reduction (Loor and Cohick, 2009). Important goals of systematic and integrated analysis are to uncover the underlying links (pathways, regulatory networks, and functional organization), and also to discover new emergent properties that may arise from examining the interactions between all components of a system (Bruggeman and Westerhoff, 2007). In the context of intra-mammary infection (IMI) and the mammary gland response, a “true” systems biology approach would encompass not only the intra-cellular networks of factors, but also communication between tissues (Piantoni et al., 2010) and between organs/systems.

Most of the genes were related to functions including Immune Response, Immune Disease, Connective Tissue Disorders, Lipid and Carbohydrate Metabolism, Molecular Transport, Cell-To-Cell Signaling, Tissue Development, Cellular Development, and Immune and Lymphatic System Development and Function. These genes encoded cytokines (CSF3, CSF2, CCL3, TNF, CCR1, IL13, KIT, CCL5, CCL4, IL10, CXCL10, IL12RB2, IFNG, IL1B, IL1A, IL8, MET, CD40, HGF, LEP, CXCL16, CX3CR1, CCR2, PDGFRA, IL12B, IL2), lipid metabolism related genes (CD36, GPX4, LPIN1, LPL, LPB), transcription regulators (BCL3, FOS and NFKBIA), receptors (TNF, IL8, RELA, TLR1, TIRAP, TLR2, NFKBIA, NFKB1, TLR4, CD40, CCL5, CXCL10, FOS, JUN, MAP3K8, IL1B, LBP, IL12B, CD14, SPP1), and others such as SELP, SELL, and SOD1, all play a role in some aspect of the immune response including cytokine activity (IL10, TNF, IL8, and IL1B), cell adhesion (SELL and SELP), immune activation (CD14 and TLR2), acute phase reaction (TNF, IL1B, and SAA3), apoptosis (BCL2, BAX). Similar results were observed by Moyes et al. (2009) in their study on bovine mammary tissue challenged with Streptococcus uberis.

A French group examined the milk somatic cell transcriptome (predominantly PMN) of mastitis susceptible and mastitis-resistant sheep infected successively by Streptococcus epidermidis (S. epidermidis) and Staphylococcus aureus (S. aureus) (Bonnefont et al., 2011). The challenges per se (i.e., S. aureus vs. S. epidermis) resulted in >5,000 DEG (FDR= 0.05 and fold change >1.5 ratio between S. aureus and S. epidermis), but the authors performed a functional analysis on those DEG with a >5-fold difference. The functional analysis of genes more expressed in S. aureus uncovered significantly-enrichment (EASE score ≤0.05) (da Huang et al., 2009a), among DEG of functions related to the immune response; whereas, genes more expressed in S. epidermidis resulted in over representation of functions related to metabolism and cell growth, indicating that more T cells were recruited upon inoculation by S. aureus than S. epidermidis. Few DEG (57 at FDR=0.05 and fold-change >1.5) were found between the mastitis-susceptible and mastitis-resistant sheep. The functional analysis indicated that mastitis-resistant sheep had a more active immune response, as indicated by the significant enrichment of “leukocyte adhesion and activation” among DEG more expressed in mastitis-resistant relative to the mastitis-susceptible sheep (Bonnefont et al., 2011).

Another German group performed direct transcriptomics comparison between primary bovine mammary epithelial cells isolated from cows with low- or high susceptibility to mastitis challenged either with E. coli or S. aureus (Brand et al., 2011). Similarly to the French study mentioned above (Bonnefont et al., 2011), the main finding from the ORA functional analysis was that the transcriptome of mammary epithelial cells from cows with low-susceptibility to mastitis challenged with either S. aureus or E. coli had a greater response to the bacterial challenge via up-regulation of genes associated with ‘acute phase response signaling’, i.e., interleukin-6 (IL6) and nuclear factor NF-kappa-B p100 subunit (NFKB2), than cows with high-susceptibility to mastitis (Brand et al., 2011). Both the German and the French studies appear to support the notion that the sensitivity of the mammary gland to bacteria plays an important role in the resistance to invading microorganisms. Perhaps the evolutionary transition of the mammary gland from a mucus-secreting surface epithelium serving a protective role to a milk-secreting gland with both protective and nutritional roles (Vorbach et al., 2006) has, as a consequence, a negative impact on the importance of the tissue’s innate immune system.

Conclusion

It was concluded that the stronger immune system of an animal would be more resistant to mastitis. In future, the analysis of this study can be used to further our knowledge about this disease which is causing a huge economical lose worldwide, however, a deeper understanding of its fundaments and development of more adequate bioinformatics tools is critically needed in order to really benefit from such an approach.

In conclusion, our network-based gene prioritization approach provides a general framework for identifying and ranking genes associated with complex diseases. To our knowledge this is the first time that protein interactions, orthologue mapping, gene expression, and literature mining were integrated for ranking candidate genes in any livestock species. The approaches and techniques that were systematically implemented in the present study are general and not confined by specific trait or species and can be applied to various complex diseases in different organisms.

Conflict of interest statement

We declare that we have no conflict of interest.

References

Bauman, D.E., Mather, I.H., Wall, R.J. and Lock, A.L., 2006. Major advances associated with the biosynthesis of milk. J. Dairy Sci., 89: 1235-1243. http://dx.doi.org/10.3168/jds.S0022-0302(06)72192-0

Bionaz, M., Periasamy. K., Rodriguez-Zas, S.L., Everts, R.E., Lewin, H.A., Hurley, W.L. and Loor, J.J., 2012. Old and new stories: Revelations from functional analysis of the bovine mammary transcriptome during the lactation cycle. PLoS One, 7: e33268. http://dx.doi.org/10.1371/journal.pone.0033268

Bonnefont, C.M.D., Toufeer, M., Caubet, C., Foulon, E, Tasca, C., Aurel, M., Bergonier, D., Boullier, S., Robert-Granie, C., Foucras, G. and Rupp, R., 2011. Transcriptomic analysis of milk somatic cells in mastitis resistant and susceptible sheep upon challenge with Staphylococcus epidermidis and Staphylococcus aureus. BMC Genomics, 12: 208. http://dx.doi.org/10.1186/1471-2164-12-208

Brand, B., Hartmann, A., Repsilber, D., Griesbeck-Zilch, B., Wellnitz, O., Kuhn, C., Ponsuksili, S., Meyer, H.H. and Schwerin, M., 2011. Comparative expression profiling of E. coli and S. aureus inoculated primary mammary gland cells sampled from cows with different genetic predispositions for somatic cell score. Genet. Sel. Evol., 43: 24. http://dx.doi.org/10.1186/1297-9686-43-24

Bruggeman, F.J. and Westerhoff, H.V., 2007. The nature of systems biology. Trends. Microbiol., 15: 45-50. http://dx.doi.org/10.1016/j.tim.2006.11.003

Huang da, W., Sherman, B.T. and Lempicki, R.A., 2009a. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protoc., 4: 44-57.

Jensen, K., Gunther, J., Talbot, R., Petzl, W., Zerbe, H., Schuberth, H., Seyfert, H. and Glass, E.J., 2013. Escherichia coli and Staphylococcus aureus induced mastitis differentially modulate transcriptional responses in neighboring uninfected bovine mammary gland quarters. BMC Genomics, 14: 36. http://dx.doi.org/10.1186/1471-2164-14-36

Kanehisa, M., 2013. Molecular network analysis of diseases and drugs in KEGG. Methods Mol. Biol., 939: 263-275. http://dx.doi.org/10.1007/978-1-62703-107-3_17

Khan, M.Z. and Khan, A., 2006. Basic facts of mastitis in dairy animals: a review. Pakistan Vet. J., 26: 204-208.

Lemay, D.G., Neville, M.C., Rudolph, M.C., Pollard, K.S. and German, J.B., 2007. Gene regulatory networks in lactation: identification of global principles using bioinformatics. BMC Syst. Biol., 1: 56. http://dx.doi.org/10.1186/1752-0509-1-56

Loor, J.J. and Cohick, W.S., 2009. ASAS centennial paper: lactation biology for the twenty-first century. J. Anim. Sci., 87: 813-824. http://dx.doi.org/10.2527/jas.2008-1375

Moyes, K.M., Drackley, J.K., Morin, D.E., Bionaz, M., Rodriguez-Zas, S.L., Everts, R.E., Lewin, H.A. and Loor, J.J., 2009. Gene network and pathway analysis of bovine mammary tissue challenged with Streptococcus uberis reveals induction of cell proliferation and inhibition of PPARγ signaling as potential mechanism for the negative relationships between immune response and lipid metabolism. BMC Genomics, 10: 542. http://dx.doi.org/10.1186/1471-2164-10-542

Paape, M.J, Bannerman, D.D., Zhao, X. and Lee, J., 2003. The bovine neutrophil: structure and function in blood and milk. Vet. Res., 34: 597-627. http://dx.doi.org/10.1051/vetres:2003024

Piantoni, P., Bionaz, M., Graugnard. D.E., Daniels, K.M., Everts, R.E., Rodriguez-Zas, S.L., Lewin, H.A, Hurley, H.L., Akers, M. and Loor, J.J., 2010. Functional and gene network analyses of transcriptional signatures characterizing preweaned bovine mammary parenchyma or fat pad uncovered novel inter-tissue signaling networks during development. BMC Genomics, 11: 331. http://dx.doi.org/10.1186/1471-2164-11-331

Rudolph, M.C., McManaman, J.L., Hunter, L., Phang, T. and Neville, M.C., 2003. 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, 8: 287-307. http://dx.doi.org/10.1023/B:JOMG.0000010030.73983.57

Schukken, Y.H., Günther, J., Fitzpatrick, J., Fontaine, M.C., Goetze, L., Holst, O., Leigh, J., Petzl, W., Schuberth, H.J., Sipka, A., Smith, D.G., Quesnell, R., Watts, J., Yancey, R., Zerbe, H., Gurjar, A., Zadoks, R.N. and Seyfert, H.M., 2011. Host-response patterns of intramammary infections in dairy cows. Vet. Immunol. Immunopathol., 144: 270-289. http://dx.doi.org/10.1016/j.vetimm.2011.08.022

Viguier, C., Arora, S., Gilmartin, N., Welbeck, K. and O’Kennedy, R., 2009. Mastitis detection: current trends and future perspectives. Trends Biotech., 27: 486-493. http://dx.doi.org/10.1016/j.tibtech.2009.05.004

Vorbach, C., Capecchi, M.R. and Penninger, J.M., 2006. Evolution of the mammary gland from the innate immune system? Bioassays, 28: 606-616. http://dx.doi.org/10.1002/bies.20423

Wiley, A.S., 2007. The globalization of cow’s milk production and consumption: Biocultural perspectives. Ecol. Fd. Nutri., 46: 281-312. http://dx.doi.org/10.1080/03670240701407657