Transcriptome analysis of SerpinB2-deficient breast tumors provides insight into deciphering SerpinB2-mediated roles in breast cancer progression

SerpinB2 is highly expressed in immune and tumor cells and is involved in multiple biological functions, including cell survival and remodeling for disease progression. This study prepared SerpinB2-deficient mice and analyzed the differentially expressed genes (DEGs) to determine if loss of this protein delays mammary tumor progression. A total of 305 DEGs (75 upregulated and 230 downregulated; > 1.5-fold difference, P < 0.05) were identified in SB2−/−;PyMT tumors compared with PyMT tumors. The DEGs were mainly involved in immune and inflammatory responses related to T cell differentiation, IFN-γ production, and lymphocyte chemotaxis based on 61 enriched GO terms, hierarchical clustering, KEGG pathways, and a functionally grouped annotation network. The significantly changed DEGs (Anxa3, Ccl17, Cxcl13, Cxcr3, IFN-γ, Nr4a1, and Sema3a) annotated with at least two GO categories in SB2−/−;PyMT tumors was validated by qRT-PCR. SerpinB2 deficiency alters the expression of multiple genes in mammary tumors, which might cause a delay in PyMT-induced mammary tumor progression.

We generated SerpinB2-deficient MMTV-PyMT (SB2−/−;PyMT) mice by intercrossing SerpinB2-deficient (SB2−/−) mice with C57BL/6 strain background MMTV-PyMT (PyMT) mice widely used to study human breast cancer. Mammary tumor onset and progression were significantly delayed in SB2−/−;PyMT mice compared with those in PyMT mice. We aimed to analyze transcriptome profiles and networks of mammary tumor tissue samples collected from age-matched PyMT and SB2−/−;PyMT mice using RNA-Sequencing (RNA-Seq) technology to understand the underlying mechanism by which SerpinB2 deficiency is involved in delayed breast cancer development and metastasis. Gene Ontology (GO) terms, pathway enrichment, and a functionally organized GO/pathway network of differentially expressed genes (DEGs; 1.5-fold change (FC) and P < 0.05) identified in SB2−/−;PyMT tumors relative to PyMT tumors were analyzed. This work provides insights into the biological functions and regulatory mechanisms of SerpinB2 in mammary tumors.

SB2−/−;PyMT mice exhibit delayed onset and growth of mammary tumors
To investigate the potential in vivo role of SerpinB2 in breast cancer development and progression, SerpinB2deficient PyMT (SB2−/−;PyMT) mice were produced by crossing SB2−/− female mice with C57BL/6 strain background PyMT males, and the PyMT transgene and SerpinB2-deficient genotype in all female offspring were analyzed by PCR (Fig. 1A). The protein expression level of ERα was not detected in 25-week-old SB2−/−;PyMT and PyMT mice tumors, while the HER2 protein level was comparably lower in SB2−/−;PyMT tumors (1.31 ± 0.17) than in PyMT tumors (2.15 ± 0.28), suggesting that SerpinB2 deficiency may be associated with reduced HER2 expression status (Fig. 1B). The raw data of Fig. 1 A, B are shown in Supplementary data 1. The time-to-first appearance of a palpable tumor per mouse (mean ± standard error), number of palpable tumorbearing mice per group (%), as well as multifocal tumor numbers per mouse per time point (mean ± standard error) were investigated. The first appearance of palpable tumors in PyMT mice and SB2−/−;PyMT mice was observed at 79.53 ± 2.98 days and 92.47 ± 3.75 days after birth, respectively (P = 0.011, Fig. 1 C). Palpable mammary tumors developed within 9-14 weeks of birth in 20-100% of PyMT mice and 11-50% of SB2−/−;PyMT mice. At 20 weeks of age, 100% of the PyMT mice and 88.24% of the SB2−/−;PyMT mice had palpable mammary tumors (P < 0.0001, Fig. 1 D). Prior to 12 weeks of age, multifocal tumors arose more slowly in the mammary glands of SB2−/−;PyMT mice than in those of PyMT mice, and the number of multifocal palpable primary tumors in PyMT mice and SB2−/−;PyMT mice at 20 weeks of age were 7.05 ± 0.34 and 4.6 ± 0.60, respectively (Fig. 1E). Among 10 mammary glands, SB2−/−;PyMT mice exhibited a remarkably slower tumor growth rate in 4th and 5th mammary glands, and the volume of tumors generated from the 4th and 5th mammary glands of SB2−/−;PyMT mice was smaller at 20 weeks compared with those of PyMT mice (Fig.1 F).

Functionally grouped annotation network
Cytoscape software functionally organized the GO network for DEGs in SB2−/−;PyMT tumors. Overall, 51 GO terms in the BP, MF, and CC categories were significantly enriched, and these categories were organized into 18 annotation networks that reflected the relationships between the terms based on the integration of their associated genes (Fig. 4). The main networks were endopeptidase inhibitor activity, lymphocyte chemotaxis, proteinaceous extracellular matrix, CD8positive α-β T-cell differentiation, and positive regulation of interferon-γ production.

Discussion
SerpinB2 is expressed by a variety of cells and is implicated in various diseases including cancer, inflammation, and immune-associated diseases [8]. Dougherty et al. produced SerpinB2-deficient (SB2−/−) mice to examine the potential in vivo functions of SerpinB2 in murine development [9]. MMTV-PyMT transgenic (PyMT) mice are widely used for studying breast cancer development and metastasis, especially the luminal B subtype [10]. There is a lack of studies analyzing the in vivo role of SerpinB2 in transgenic and knockout animal models of breast cancer. We generated SerpinB2-deficient MMTV-PyMT (SB2−/−;PyMT) mice by crossing SB2−/− mice with PyMT mice to provide insights into the in vivo function of SerpinB2 in mammary tumor development and progression and found that SerpinB2 deficiency resulted in a delay in PyMT-induced mammary tumor initiation and progression. RNA-Seq is a valuable method for quantifying transcriptome signatures and deciphering the regulatory network underlying breast tumor progression from primary tumor development to metastasis [11,12]. Transcriptome profiling using RNA-Seq was performed in mammary tumors of PyMT and SB2−/−;PyMT mice to obtain the DEGs, the genetic network, and understand how SerpinB2 deficiency affects delayed breast cancer development and metastasis.
Three hundred and five DEGs were observed in SB2−/−;PyMT tumors, including 75 upregulated genes and 230 downregulated genes compared with PyMT tumors. SerpinB2 subfamily members (SerpinB1-13) inhibit activity of cytotoxic and apoptotic proteases and papain-like enzymes and non-inhibitory molecules [13]. SerpinB8, SerpinB11, and SerpinB5 were among the top 10 upregulated DEGs in SB2−/−;PyMT tumors. Particularly, SerpinB5 (maspin) is a tumor suppressor that inhibits metastasis in human mammary epithelial cells     mammary epithelium. Amtn (amelotin) may play a role in the controlled mineralization of hydroxyapatite that is an early diagnostic marker for breast cancer, is involved in enhancing breast cancer progression, and is associated with poor breast cancer prognosis [15]. Therefore, upregulated SerpinB5 and downregulated Amtn may be involved in suppressing mammary cancer development and lymph node metastasis in SB2−/−;PyMT mice. The role of SerpinB2 expression by diverse types of cells in cancer growth and metastasis is controversial. High SerpinB2 expression is noted in mouse melanoma cell line of B16 cells, and human melanoma cell line of M24met inhibiting metastasis [16,17]. In addition, perivascular expression of SerpinB2 in C8161 human melanoma cells promotes brain metastasis [18]. The Ser-pinB2-deficient human lung cancer cell line H2030-BrM3 and human breast cancer cell line MDA231-BrM2 do not form brain metastases because the deleterious effects of SerpinB2 on the vascular attachment and survival of cancer cells are inhibited [19]. Recently, Jin et al. showed that SerpinB2 knockdown in MDA-MB-231 breast cancer cells suppressed their migratory activity and lung   [20]. SerpinB2-deficient mouse embryonic fibroblasts exhibit increased pancreatic tumor growth and local invasion with reduced collagen deposition [4], whereas fibroblasts overexpressing SerpinB2 reduce human breast tumor cell apoptosis. MDA-MB-231 promotes breast tumor growth and lung metastasis and reduces tumor apoptosis [21]. SerpinB2 upregulation in monocytes/macrophages following infection or stimulation with inflammatory mediators is involved in a cellular response to delay cell death, possibly allowing cells to complete vital functions such as lymphocyte activation and antigen presentation [22]. SerpinB2-deficient mice produce more IgG2c and OVA-specific IFN-γ-secreting T cells than wild-type mice suggesting that SerpinB2 regulates adaptive immunity [3]. In the aforementioned studies [19,20], SerpinB2 expression substantially mediates physiological and pathological functions including inhibitory activity against serine protease plasminogen activators, ECM remodeling, cell survival and migration, inflammation, and immunity. Thus, the role of SerpinB2 in cancer development and aggressive progression may depend on the cancer type.
Our results showed that SerpinB2 plays roles in immunity, inflammation, chemotaxis, and ECM modulation. GO enrichment and functional network analysis of DEGs identified from SB2−/−;PyMT tumors revealed multiple biological processes including T cell differentiation, INF-γ production, lymphocyte chemotaxis, ECM, and peptidase inhibitor activity which is consistent with previous studies [2,3,13,23,24]. Eleven DEGs (Anxa3, Ccl17, Cxcl2, Cxcl13, Cxcr3, INF-γ, Itgad, Nr4a1 Sema3a, Tnfsf14, and Trem1) were annotated with at least two GO categories. We here validated quantitative expression levels of these genes between PyMT and SB2−/−;PyMT tumors using qRT-PCR. Seven of the 11 genes (Anxa3, Ccl17, Cxcl13, Cxcr3, IFN-γ, Nr4a1, and Sema3a) exhibited identical RNA-Seq and qRT-PCR results, whereas four genes (Cxcl2, Itgad, Tnfsf14, and Trem1) did not. Thus, the RNA-Seq and qRT-PCR findings agreed on 63.6% of the genes. RNA-Seq is a strong tool, however variations in starting quantities or cDNA quality can obscure quantitative differences between two samples [25]. As a result, more precise methods are needed to confirm the exact copy counts of specific mRNAs. The only approach to avoid this error is to use internal standard-normalized qRT-PCR.
Gene of Ccl17, Cxcl13, Cxcr3, IFN-γ, and Sema3a, were significantly decreased while the Anxa3 and Nr4a1 were significantly upregulated in SB2−/−;PyMT mammary tumors. High levels of Ccl17, Cxcl13, and Cxcr3 in mammary tumors or tumor-associated macrophages induce cancer cell migration which is associated with metastasis disease [26][27][28][29]. Ccl17 and Cxcr3 increased in mammary tumors that protects cancer cell survival and induce tumor growth to promote breast cancer lung metastasis [28,29] Hypoxia-induced chemokine Sema3a stimulates tumor-associated macrophages to the alternative type (M2), leading to pro-tumor immunity [30]. And the tumor suppresser gene of Nr4a1 in TNBC was decreased cancer cell proliferation and invasiveness [31]. These studies support our findings that low levels of Ccl17 and Cxcl13, Cxcr3, and Sema3a as well as high levels of Nr4a1 in SB2−/−;PyMT tumors may delay mammary tumor development and may lead to decreased metastasis. However, in SB2−/−;PyMT tumors, a high levels of IFN-γ, which activates the antitumor immune response [32], and a low levels of Anxa3, which controls cancer cell invasion and lung metastasis [33], do not support SerpinB2 deficiency causing a delay in PyMTinduced mammary tumor progression. Given these results, the crosstalk between the immune system and cancer cells mediated by these cytokines and chemokines in SB2−/−;PyMT tumors during tumor initiation and progression in breast cancer requires further study.

Conclusion
SerpinB2-deficient MMTV-PyMT mice (SB2−/−;PyMT) exhibited altered expression of multiple genes related to many aspects of immunity, inflammation, chemotaxis, cell adhesion, ECM modulation, peptidase activity, and cell proliferation suggesting that they delay PyMTinduced mammary cancer development and metastasis. Further studies are required to understand the complex roles of multiple genes identified in our study in tumor microenvironments. University (SNU-150210-3-4). All methods are reported in accordance with ARRIVE guidelines (https:// arriv eguid elines. org) for the reporting of animal experiments.

RNA extraction from tumor tissues
Mammary tumor tissues were collected from age-matched SB2−/−;PyMT (n = 10) and PyMT (n = 10) mice (20 weeks old [n = 2], 22 weeks old [n = 2], and 25 weeks old [n = 6] per group). Total RNA content was extracted from a 20-mg piece of the tumor using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity number (RIN) was determined using an Agilent RNA 6000 Nano kit following the manufacturer's protocol on an Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA).

RNA-Seq and DEG analysis
Sequencing libraries were constructed using the Quant-Seq 3′ mRNA-Seq library prep kit (Lexogen, South Morang Victoria, Austria) according to the manufacturer's instructions. High-throughput RNA-Seq was performed for 75 single-end sequences using the Next-Seq 500 system (Illumina, San Diego, CA, USA). DEGs between SB2−/−;PyMT, and PyMT tumors were determined based on counts from unique and multiple alignments using coverage in BEDTools [35]. The read count data were processed based on the quantile normalization method using EdgeR within R (R Development Core Team, 2016) by Bioconductor [36]. Genes with a P value < 0.05 were ranked by the log10 P value and plotted against the log2 (FC) in a "volcano" plot. Genes that were upregulated and downregulated with a P value < 0.05 and log ratio >1.5 were considered DEGs.

GO, KEGG pathway, enrichment, and functional network analysis
Biological functional categories enriched in DEGs were identified using the functional annotation and clustering tool of the Database for Annotation, Visualization, and Integrated Discovery v6.7 (https:// david. ncifc rf. gov/) [37][38][39]. Significant GO terms were identified after multiple testing adjustments using the Benjamini-Hochberg method; Benjamini < 0.05, which indicated a statistically significant difference and GO term lists BP, MF, and CC were matched [40]. KEGG pathway enrichment analyses was performed using the KEGG pathway software (https:// www. kegg. jp/ kegg/) from the Kanehisa laboratory and the threshold value was set at P < 0.05. The functionally organized GO/pathway network was created using Cytoscape software version 2.6.2 (http:// www. cytos cape. org/) with the ClueGO plugin (http:// www. ici. upmc. fr/ cluego/ clueg oDown load. shtml) [41].

qRT-PCR
Experiments were performed using the same samples extracted from the tumor tissues of age-matched PyMT and SB2−/−;PyMT mice for RNA-Seq. qRT-PCR was performed using an ABI PRISM 7900 system, SYBR Green PCR master mix (Applied Biosystems, Foster City, CA), and specific primer sets for Anxa3, Ccl17, Cxcl2, Cxcl13, Cxcr3, INF-γ, Itgad, Nr4a1 Sema3a, Tnfsf14, and Trem1 (Supplementary Table 4). The results were analyzed using the comparative ∆Ct method with the relative gene expression normalized to the γ-actin housekeeping gene [42].

Statistical analysis
Statistical analyses were performed by GraphPad Prism 8.0 (GraphPad Software, Inc., La Jolla, CA). Statistical significance was determined by Student t-test to compare tumor size and number differences, and gene expression between two groups. For all tests, a P-value less than 0.05 was considered statistically significant.