Identification of the molecular relationship between intravenous leiomyomatosis and uterine myoma using RNA sequencing

The purpose of this study was to explore the potential relationship between intravenous leiomyomatosis (IVL) and uterine myoma (UM) at the molecular level. RNA-sequencing was performed on IVL tumours, UM tumours, and adjacent normal uterine muscle. We compared the gene expression levels between IVL and normal uterine muscle, UM and normal uterine muscle, to identify differentially expressed genes (DEGs). Then we used Gene Ontology Enrichment Analysis to determine the functions of the DEGs and performed specimen cluster analysis. We obtained 98 DEGs between IVL and adjacent normal uterine muscle, and 61 DEGs between UM and adjacent normal uterine muscle. Functional enrichment of both IVL and UM DEGs showed that they are associated with hormone stimulus, extracellular matrix, and cell adhesion. Unsupervised clustering analysis showed that IVL and UM could not be separated completely. Among these dysregulated genes, we found that HOXA13 showed a distinct dysregulated status between IVL and UM. HOXA13 may therefore serves as a biomarker to distinguish IVL and UM. Our results showed that IVL and UM may have similar dysregulated gene networks. They may be closely related, and HOXA13 may serves as a biomarker to distinguish between IVL and UM.


Results
Heterogeneous gene expression in IVL. The morphology of IVL is often irregular. In addition, a common phenomenon of IVL is intracardiac extension with the potential for lung metastasis. This phenomenon indicates that the tumour may evolve while growing along the postcava. To address this issue, we conducted multi-position sampling for each tumour (the distal end, the midpiece, and the proximal end) and then conducted RNA sequencing for each sample. Compared with the paired myometrial tissues, we obtained 98 DEGs (Supplement materials) from the three parts of IVL (Fig. 2). Among the DEGs, 24.5% were shared between the proximal end and the midpiece of IVL, while only 5.1% were shared between the distal end and the midpiece (Fig. 3A). These results indicated the existence of extensive heterogeneous gene expression in IVL.

Dysregulated gene networks in IVL.
The results of Gene Ontology analysis of the differentially expressed genes showed that the products of these genes were significantly associated with the extracellular region, acting in response to hormone stimulus, cell-cell adhesion, and in the regulation of cell growth and programmed cell death. Specifically, both the up-and down-regulated genes were particularly associated with response to hormone stimulus, extracellular matrix, cell adhesion, programmed cell death (Fig. 4A). This indicated that IVL  Dysregulated gene networks in uterine myoma. Owing to similarities in the clinical features of uterine myoma and IVL, we explored the dysregulated gene networks of uterine myoma. We conducted RNA sequencing for uterine myoma and paired myometrium, and obtained 61 DEGs (Fig. 5). Gene Ontology analysis results showed that these DEGs were particularly associated with response to steroid hormone stimulus, extracellular regions, and cell adhesion (Fig. 4B).
IVL could be a subclass of uterine myoma. To further explore the relationship between IVL and uterine myoma, we compared the DEGs found in two diseases. As shown in Fig. 3A, only seven dysregulated genes were shared between IVL and uterine myoma and the IVL proximal end shared two DEGs with UM: CYR61 and HOXA13. When we took the dysregulation status into consideration, we found that HOXA13 is not included among the up-regulated DEGs (Fig. 3B) or down-regulated DEGs (Fig. 3C), showing that IVL and uterine myoma exhibited a difference in their dysregulation status. Next, we carried out real-time PCR of HOXA13 (Sangon Biotech; HOXA13 forward primer sequence: 5′-CAC GAA CCC TTG GGT CTT C-3′; HOXA13 reverse primer sequence: 5′-TCT TTG GGG CAG TAC ATT TGG-3′) to ascertain its expression level in IVL, UM, and normal uterine muscle (Fig. 6). The obtained findings revealed that HOXA13 may serve as a biomarker to distinguish IVL and uterine myoma. We next explored the functional characteristics of the genes dysregulated in IVL and uterine myoma. We found that the differentially expressed genes of IVL and uterine myoma are both related to the dysregulation of steroid hormone stimulus, extracellular matrix, and cell adhesion. This indicated that IVL and uterine myoma might have similar dysregulated gene networks.
We extracted the expression profile of dysregulated genes in both IVL and uterine myoma and performed cluster analysis. As shown in Fig. 7, all myometrial samples were clustered together. However, IVL and UM could hardly be distinguished from each other. These results indicated that some IVLs are similar to UM.

Discussion
IVL is a rare, benign, smooth muscle cell tumour. Few studies have been conducted to investigate its molecular changes or its relationship with uterine myoma. In this study, we employed RNA sequencing to explore the dysregulated gene networks in IVL and uterine myoma, to shed light on their relationship. Since first IVL was identified, 200 cases have been reported in the literature 14 . IVL patients are typically middle-aged women, of whom approximately 50% have undergone a hysterectomy 15 . IVL is pathologically classified as a benign tumour, but behaves like a malignant lesion. IVL can extend into the atrium, leading to fatality 16 , or metastasize to the pulmonary artery 17 . However, with complete excision, the probability of recurrence is close to zero 15,18 .
The diagnosis and treatment strategy of IVL have been standardized, but its pathogenesis remains uncharacterized. Moreover, the origin of the tumour and the relationship between IVL and uterine myoma remains unknown. There are two main hypotheses about the origin of IVL. One suggests that it originates from the smooth muscle of the vessel wall, while the other suggests that it occurs due to uterine myoma extending through a vein 19 . For each of these theories, there are lots of supportive evidences. One report of a pathological examination of the tumour revealed that it is benign and consists of proliferating smooth muscle fibres (rate of proliferation < 1%) with normal mitotic activity. Immunohistochemical analysis showed low levels of oestrogen receptor expression (around 30-40%), while progesterone receptor was strongly expressed (80%). These results suggest that IVL growth is connected to hormone levels 20 , which resembles the pathogenesis of uterine myoma.
Cytogenetic studies carried out to investigate IVL and uterine myoma suggested that spontaneous chromosome rearrangements may be responsible for the growth of uterine leiomyoma 21 . Leiomyosarcomas present  Bradley et al. studied the molecular cytogenetics of IVL, comparing it with uterine leiomyoma. They demonstrated that the pathogenesis of IVL is related to typical uterine leiomyoma and has a similar chromosomal mechanism. About half of uterine leiomyomas have structural chromosomal abnormalities, 20% of which show rearrangements of 12q14-15 targeting the gene encoding the high-mobility group AT-hook 2 (HMGA2) 24 . As this chromosomal rearrangement is also seen in IVL, some researchers think that uterine leiomyoma and IVL may represent different stages of a single disease. However, the uniclonal origin of IVL distinguishes it from multiclonal uterine leiomyoma 25 .
Another study used PCR sequencing to analyse nine IVL cases. The results demonstrated that all nine cases of IVL had a wild-type MED12 gene at codon 44, suggesting that IVL is a distinct smooth muscle tumour with a unique pathogenesis different from that of conventional leiomyomas. (The latter are associated with MED12 mutation.) Sekine et al. used Sanger sequencing to determine the prevalence of MED12 mutations in smooth muscle tumour of different organs. They found that, in contrast to uterine tumours, extrauterine smooth muscle tumours including leiomyomas, leiomyosarcomas, and angioleiomyomas do not have MED12 mutations [26][27][28] .
In summary, the relationship between IVL and uterine myoma has not been characterized, with no reported studies involving RNA sequencing. The current study may be the first to use RNA sequencing to investigate the origin of IVL. Our results showed that IVL and uterine myoma have similar dysregulated gene networks; both diseases are related to the dysregulation of steroid hormone stimulus, extracellular matrix, and cell adhesion. Moreover, clustering analysis showed that it is difficult to separate IVL from uterine myoma. These results indicated that IVL and uterine myoma might have similar origins. However, we found that only seven dysregulated genes were shared between IVL and uterine myoma. In addition, some genes playing a key role in regulating programmed cell death are significantly upregulated in IVL, but not in uterine myoma. Among these dysregulated genes, HOXA13 was shown to have a distinct dysregulated status between IVL and uterine myoma. These results indicated that IVL is not identical to conventional uterine myoma. In addition, HOXA13 may serve as a biomarker to distinguish IVL and uterine myoma.

Conclusion
We concluded from this study that IVL is closely related to uterine myoma. While other study indicated that IVL cases may share some molecular cytogenetic characteristics with uterine leiomyoma 29 , IVL may be one subtype of uterine myoma, it is not identical to conventional uterine myoma. The HOXA13 gene showed a distinct dysregulated status between IVL and uterine myoma and may serves as a biomarker to distinguish these two diseases.
The results of this study demonstrated that there is a close relationship between IVL and uterine myoma, but it still requires further exploration.

Materials and Methods
Patients and specimens. Fresh tumour samples and paired adjacent normal tissue were obtained from patients who underwent surgery at Peking Union Medical College Hospital between 2014 to 2016. Fresh samples were frozen with liquid nitrogen and stored at −80 °C. All patients were clinically diagnosed and the results were confirmed by two pathologists. IVL specimens were divided into three segments according to their morphology, namely, the proximal end, midpiece, and distal end, and then numbered sequentially ( Data analysis. Clean reads were mapped to human reference genome version 19 (hg19) using TopHat2 30 .
Gene expression levels were calculated by Cufflinks-2.0.2 31 , and differentially expressed genes (P value < 0.05, q value < 0.05, fold change 2) were identified in accordance with the recommended procedure by Cuffdiff-2.0.2 32 . The heatmap of DEG was constructed using R software-3.4.3. We performed Gene Ontology enrichment analysis using The Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.7 (https://david-d.ncifcrf. gov/home.jsp) 33 and the results were visualized using Cytoscape-2.8.0 34 . Cluster analysis was performed using the Cluster R package. RT-PCR data were analysed using GraphPad software. Results are expressed as mean ± SEM. A value of P < 0.05 was considered to be statistically significant. There were no non-default parameters were used for the analyses. Be aware that the software is updated on schedule, please use the latest version.

Data Availability
The datasets generated during and analysed during the current study are available in The European Genome-phenome Archive.The accession ID:EGAS00001002504.