Breast cancer is a heterogeneous disease. The identification of molecular subtypes (Perou et al., 2000) was a major breakthrough in order to categorize this heterogeneity, made possible by the emergence of whole-genome microarray technology (Rueda, 2014). Molecular subtypes associate expression of certain markers to phenotypical differences in cancer, pointing to different prognosis as well as distinct therapies (Liu et al., 2014).

The identification of breast cancer molecular subtypes is a paramount example of the impact of high throughput technologies in the study of cancer. These technologies are able to provide even deeper biological understanding, when analyzed by an integrative, systemic framework. In this regard, network theory has emerged as a major tool to achieve this goal (Newman, 2010; Chen et al., 2015) under a systems biology view of matters.

A network is a mathematical construct composed by a set of nodes or vertices, and a set of links that represent a relation between them. A biological network, is a network where nodes represent any kind of biological molecules: genes, transcripts, proteins, metabolites, etc., and links represent physical or chemical interactions between those molecules (Jeong et al., 2000; Hasty et al., 2001; Jeong et al., 2001; Thattai and Van Oudenaarden, 2001; Lee et al., 2002; Maslov and Sneppen, 2002; Davidson and Levin, 2005; Guimera and Amaral, 2005; Levine and Davidson, 2005; Davidson and Erwin, 2006). With gene expression microarray technologies, it is factible to construct transcriptional networks where nodes are transcribed genes, and links represent a correlation between expression values of said genes, which point to a possible interaction between them at the transcriptional level (Tovar et al., 2015).

Many correlation measures have been implemented in order to construct biologically meaningful transcriptional interaction networks based on the inference of statistical dependency (Friedman et al., 2000; Gardner et al., 2003; Giuliani et al., 2004; Wang et al., 2005; Cowell, 2006; Nielsen and Jensen, 2009). This is especially fitting in view of the isomorphism existing between a network structure and a correlation matrix whose elements are the strength of the interaction between the intervening nodes. We must notice, however, that this correlation structure is usually given in the presence of an accompanying variance structure among gene expression levels (Giuliani et al., 2004). It has long been proven that the best estimator of statistic dependency is mutual information (MI) (Basso et al., 2005; Margolin et al., 2006; Hernández-Lemus and Rangel-Escareño, 2011; de Matos Simoes and Emmert-Streib, 2012; Hernández-Lemus and Siqueiros-García, 2013). These statistically inferred networks provide a deeper level of biological understanding in two main directions: giving support to previously identified biological observations, and giving new insights regarding novel biological interactions.

Network structural properties have been related to features in the biological context (Kitano, 2002; Albert, 2005; Serrano, 2007; Hakes et al., 2008). Therefore, analyzing these properties in transcriptional networks may provides us better understanding of the underlying biological phenomena. Global network metrics often provide information regarding the system as a whole; while local parameters provide information regarding the relevance of particular nodes (Barabasi and Oltvai, 2004; Newman, 2010; Barabási et al., 2011; Biane et al., 2016; Robinson and Nielsen, 2016). The transcriptional network approach has proven be useful to unveil transcriptional regulation in cancer (Carro et al., 2010; House et al., 2010; Pe'er and Hacohen, 2011; Madhamshettiwar et al., 2012) and in particular in breast cancer (Van De Vijver et al., 2002; Lim et al., 2009; Cicatiello et al., 2010; Gu et al., 2010; Tovar et al., 2015; Castro et al., 2016).

Transcriptional networks are representations of the regulatory programs behind phenotypes. Given the intrinsic heterogeneity of breast cancer molecular subtypes, a fundamental question which remains unsolved is whether the transcriptional architecture of these subtypes is different. To answer this we constructed transcriptional networks for breast cancer molecular subtypes based on mutual information of genome-wide gene expression. We compared them to a network of healthy mammary tissue.

We identified differences in network architecture between phenotypes. We observed major differences between the cancer subtype networks and the non-tumor network. Particular architectural features were associated to the different molecular subtypes. We find that in these networks, the connectivity of particular genes may indicate differences of their role in the transcriptional program of each subtype. Identifying such differences may be key to understand how the specific transcriptional program shapes a particular phenotype. This in turn, will enhance our insight on the nature of molecular subtypes, with basic and clinical implications.

We inferred transcriptional networks for each breast cancer molecular subtype and for non-tumor mammary tissue. Network approaches are highly relevant for the understanding of the connection between sample state-variability and gene expression patterns, essential to elucidate the role that such expression patterns play in the establishment of cellular phenotypes. In this work we decided to pursue this by calculating mutual information correlation measures rather that parametric correlation coefficients due to the higher generality of the former.

Based on the evaluation of network threshold (Supplementary Material 2), we will present the results for networks constructed from the 10,000 interactions ranked with highest MI values, which are also, as previously mentioned, those with the highest statistical significance (p-value < 10⁻¹⁰). This cut-off value lies in the identified range (between the 10³ and 10⁵ highest ranked interactions) in which the node-to-link ratio is consistent with the expected values for biological complex networks.

Deconvoluting transcriptional network architecture is a step toward uncovering differences between health and disease, as well as between different manifestations of disease. Doing this is relevant for two main lines of research: basic biological knowledge and clinical therapeutical interventions. On the one hand, unveiling the hidden relationship between system's level transcriptional programmes—as sketched by gene regulatory networks—and molecular and physiological phenotypes—given by the different breast cancer subtypes—is of great importance on the way to understand to a deeper level the elusive genotype/phenotype relationships (Barabási et al., 2011; Censi et al., 2011; Ahmad et al., 2012). On the other hand, detailed and specific knowledge about the differences and similarities of phenotype-specific gene regulatory networks, and their underlying biochemical processes as represented by metabolic and signaling pathways, may pave the way to find specific pharmacological targets useful in the clinical and therapeutical management of different breast cancer subtypes.

In this work, by using a computational approach based on information theory, we inferred transcriptional networks from genome-wide gene expression microarrays of the following breast cancer subtypes: Luminal A, Luminal B, Basal, and HER2-enriched; these were contrasted with comparable non-tumor breast tissue networks. We showed that the architecture of these transcriptional networks is unique to each phenotype, and emphasize the major differences between non-tumor and cancer networks, as well as the differences between different subtypes of breast cancer. We propose that these transcriptional networks reflect the transcriptional programs behind each phenotype, which is involved in the physiopathological differences observed in breast cancer. To our knowledge, this is the first time that transcriptional network structure has been characterized at any level of description for these clinically important breast cancer subtypes.

Breast cancer heterogeneity is a widely acknowledged fact, with molecular, histological, and clinical manifestations. In this work, we have shown breast cancer heterogeneity manifested in transcriptional network architectures; since molecular subtypes are defined by differences in gene expression profiles (see Supplementary Material 2), we could expect that transcriptional networks, derived from such data, would reflect this. It is well-known that network structure is intimately linked to functionality in biological networks (Kitano, 2002; Albert, 2005; Serrano, 2007; Hakes et al., 2008). In this case, we propose that our inferred networks are a representation of an underlying transcriptional program associated to each of the studied phenotypes.

We expected a drastic difference between non-tumor and tumor networks, just as there is a drastic difference between non-tumor and tumor tissue. Indeed, we found a marked difference in structure between breast cancer networks and a non-cancer network; this is in agreement with the evident differences in gene expression profiles between these two states (see Supplementary Material 2). The non-tumor transcriptional network is dominated by a giant component, while the networks for each of the breast cancer subtypes studied present a larger number of disconnected components. This suggests the existence of generalized transcriptional communication in healthy cells, which is lost and supplanted with a fractured, more autonomous regulation in different cancer manifestations. We consider that more research is needed in order to find mechanistic causes to this phenomenon.

The differences between cancer and non-cancer networks are also seen in the genetic composition of these networks, suggesting differences in the importance of genes in the regulatory programs of health and disease. There are only 52 genes which are part of both the non-tumor and all the cancer networks. Most of these are located in small “islands” in the network, where they interact with few other genes, rendering their impact to the global network topology negligible. Only 5 out of these 52 genes belong to the larger network components. These genes seem to be associated to general biological functions:

Breast cancer heterogeneity, as reflected by molecular subtypes, is mostly related to the cellular composition of the normal mammary tissue: the mammary epithelium has an inner layer composed by luminal cells, and an outer layer formed by basal cells (Skibinski and Kuperwasser, 2015); the luminal and basal subtypes show phenotypical similarities to these cells. HER2-enriched breast cancer is mostly determined by an overexpression of the ERBB2 receptor, regardless of whether its cell exhibits luminal or basal characteristics. With this in mind, we expected transcriptional networks to reflect this. It would be expected to find common features in transcriptional networks for the two luminal subtypes, and also to think that these networks would be somewhat different from the basal network, with the HER2-enriched network perhaps exhibiting commonalities with all networks. Our results show the expected similarity between the luminal A and B transcriptional networks. However, it was surprising to find that the highest similarity in both nodes and links was not between them, but rather between the luminal B and Basal networks. Meanwhile, the HER2-enriched subtype network showed was found to have the most unique architecture, as well as being the most different of all breast cancer networks. It is worth noting that we were able to recover a genetic network composition that captures the differences between normal breast and breast cancer tissue, as well as the differences between the various molecular subtypes, at the gene expression level (see Supplementary Material 2).

We believe that these results suggest that the oncogenic processes behind breast cancer originate unique transcriptional programs which drive each molecular subtype. The similarities of the luminal B subtype network to both the luminal A and the Basal networks are in line with the idea of a common breast cancer progenitor cell for these subtypes; at the same time, the nature of the HER2-enriched network could be indicative of a different molecular origin. (Skibinski and Kuperwasser, 2015).

The existence of a “core set” of genes which are shared among all breast cancer subtype networks was noteworthy. More interesting was the fact that even though all subtypes share these genes, the wiring structure of these is not the same in each phenotype, which again suggests that differences in the regulatory program are driven not only by the genes that participate, but by the relationships between them. The transcriptional networks define a landscape in which different elements may be playing distinct roles. Therefore, our breast cancer transcriptional networks not only provide us with global network features, but are also useful in order to identify the influence that any given gene may have in the transcriptional program. In order to illustrate how the same genes may have different roles in different breast cancer manifestations, we selected, based on their degree in each breast cancer network, the CNR2, LCK, and LUZP4 genes.

CNR2 is a gene highly connected in all breast cancer networks. CNR2 codifies the cannabinoid receptor 2 (CB2). This receptor is associated to immunomodulation and related processes by endocannabinoids (Munro et al., 1993). CB2 alterations have been found in different types of cancer (Guida et al., 2010; Jha et al., 2012; Pisanti et al., 2013), including breast cancer (Nasser et al., 2011; Pérez-Gómez et al., 2015; Sophocleous et al., 2015). Our analysis identifies CNR2 as one of the top 30 highest degree nodes in our breast cancer subtype transcriptional networks. Based on this, CNR2 might be pointed as an important gene in the general breast cancer transcriptional architecture. Cannabinoid receptors have been previously proposed as pharmacological targets for cancer (Chakravarti et al., 2014; Velasco et al., 2016), including breast cancer (Qamri et al., 2009; Morales et al., 2015). Our transcriptional network findings suggest that, if therapeutic benefits to this type of treatment are found, they may be of use to all types of breast cancer.

LCK gene was identified as a key player in HER2-enriched subtype while having a low degree in the other molecular subtypes network. LCK codifies the LCK proto-oncogene, Src family tyrosine kinase, a protein involved in signal transduction. LCK has been found expressed in breast cancer (Köster et al., 1990). Furthermore, a role in breast cancer progression and angiogenesis has been identified (Chakraborty et al., 2006). Interestingly enough, there are reports of a LCK-associated molecular signature with prognostic utility in HER2-enriched breast tumors (Rody et al., 2009). The result presented here reinforces the functional implication of LCK in the context of HER2-enriched breast cancer and emphasizes the necessity of further, focused studies.

LUZP4 was identified as a high degree gene in all breast cancer networks, except in the HER2-enriched. LUZP4 codifies a leucine zipper protein. This protein has not been extensively explored, however, recently, this leucine zipper has been identified as involved in mRNA exporting in cancer cells (Viphakone et al., 2015). Our results indicate that this gene may be an important player in Luminal A, Luminal B, or Basal. Furthermore, unlike the two previously discussed genes, LUZP4 has not been extensively studied in the context of breast cancer. This gene is an example of the type of new biological information that can be recovered from existing data only with the use of network-based approaches. We believe that further experimental exploration of this molecule may be of interest in the future.

As it has been largely discussed (Margolin et al., 2006; Baca-López et al., 2012), the probabilistic inference of gene regulatory networks is fundamental (and in some cases almost mandatory) to unveil complex transcriptional interactions that would be otherwise extremely difficult to notice. However, as it has also been noticed (Margolin et al., 2006), reverse engineering methods have also a number of limitations. In the particular case of transcriptional networks inferred via Mutual Information calculations, one important aspect is that regulatory interactions may include not only canonical transcription factor to target interactions, but also a number of indirect and more complex relationships. In this sense, MI-inferred networks (such as the ones discussed here) provide us with information about statistical dependency in the transcriptional profiles (Hernández-Lemus et al., 2009; Baca-López et al., 2012; Creixell et al., 2015).

A major open question in the field of network reconstruction is still how to establish a proper cut-off value to include a link in a network (Serrano et al., 2009). We decided to make this choice based on two constraints: on one hand, to include links with a high statistical confidence, as measured with the associated p-value provided by ARACNE, and on the other one, to obtain a network that belongs to the complex network regime, as described by Albert and Barabási (2002), which exhibit a node-to-link ratio in which nodes << links.

It is evident that the threshold choice will modify network parameters. We identified that, in an interval surrounding our selected node-to-link ratio (0.1), network parameters are generally stable (Supplementary Material 3). More importantly, even though the actual numeric values may be modified, is important to notice that the behaviors such as similarity between networks (Supplementary Material 4) and the relative influence of genes, based on their node degree (Supplementary Material 5) are preserved in this interval. Importantly, this shows robustness of the transcriptional regulatory program, while at the same time allowing us to observe structural differences between phenotypes along the interval. It is worth noticing that the network structures obtained with our methodology are quite different to networks generated at random (see Supplementary Material 6).

We have shown how a network approach can be useful to understand the heterogeneity of breast cancer. We were able to infer and compare the transcriptional programs of breast cancer molecular subtypes, and contrast them with that of healthy mammary tissue. We showed how this paradigm can help to identify novel roles of molecules in different manifestations of breast cancer, with potential biomedical implications.

GdA devised the experiment, implemented the methodologies and provided this manuscript's first draft. TV contributed to the discussion and visualization of results and helped in the revision of the manuscript. JE contributed to the analysis of results, discussed biological implications, and helped in the writing of the manuscript. EH proposed the experiment, participated in its design and coordination and helped in the writing of the manuscript. All authors read and approved the final manuscript.