Articles

Single-cell RNA sequencing highlights the functional role of

human endogenous retroviruses in gallbladder cancer
Jinghan Wang,a,f Meng Ren,b,c,f Jundan Yu,b,d Mingtai Hu,a Xiaojing Wang,b,d Wencong Ma,a Xiaoqing Jiang,e,** and
Jie Cuib,c,*
a
Department of Hepatobiliary and Pancreatic Surgery, Shanghai East Hospital, Tongji University School of Medicine, Shanghai 200120,
China
b
CAS Key Laboratory of Molecular Virology and Immunology, Institut Pasteur of Shanghai, Chinese Academy of Sciences, Shanghai
200031, China
c
Nanjing Advanced Academy of Life and Health, Nanjing 211135, China
d
University of Chinese Academy of Sciences, Beijing 100049, China
e
Department of Biliary Tract Surgery I, the Third Hospital of Naval Medical University, Shanghai 200438, China

Summary eBioMedicine
2022;85: 104319
Background Gallbladder cancer (GBC), the most common malignancy of the biliary tract, shows late diagnosis and
low survival rate and requires continued search for new diagnostic biomarkers and therapeutic targets. Human Published Online XXX
[Link]
endogenous retroviruses (HERVs) are specifically prone to be reactivated in diverse cancers and are implicated in
1016/[Link].2022.
cancer progression and immunotherapy. 104319

Methods Single-cell RNA sequencing was performed on tumor tissues and paired adjacent tissues from 4 GBC
patients. Dual-luciferase reporter assay was applied to measure enhancer activity of HERV sequences.

Findings We dissected the cellular diversity and described the HERV transcriptomic landscape for GBC. We found
that HERVs were transcribed in a cell type-specific manner and different HERV families were associated with diverse
biological effects. HERVs could function as enhancers, presumably causing altered expression of neighboring genes.
The transcription level of HERVH was gradually elevated with the malignant transformation of epithelial cells,
suggesting HERVH may be a potential early diagnostic biomarker of GBC. HHLA2, a newly emerging immune
checkpoint, was derived by HERVH, exhibited an expressional correlation with HERVH, and was identified as a
promising target for immunotherapy.

Interpretation Exploring the transcriptional landscape and potential functional impact of HERVs highlights the
important role of HERVs in GBC and provides a fresh perspective on managing GBC.

Funding This study was supported by the National Natural Science Foundation of China (31970176, 81972256) and
the research grants from the Innovation Capacity Building Project of Jiangsu province (BM2020019).

Copyright © 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND
license ([Link]

Keywords: Gallbladder cancer; Single-cell RNA sequencing; Human endogenous retrovirus; Enhancer; Immune
checkpoint; HERVH

Abbreviations: GBC, gallbladder cancer; HERV, human endogenous retrovirus; scRNA-seq, single-cell RNA sequencing; TME, tumor microenviron-
ment; WTA, whole transcriptome analysis; DEG, differentially expressed gene; CNV, copy number variation; GO, gene ontology; NK cell, natural killer
cell; NKT cell, natural killer T cell; DC, dendritic cell; ICS, intermediate cell state; HHLA2, human endogenous retrovirus-H long terminal repeat-
associating 2; CD4+ Th cell, CD4+ T helper cell; IgG, immunoglobulin G; cDC, conventional DC; mo-DC, monocyte-derived DC; CAF, cancer-
associated fibroblast; ECM, extracellular matrix; iCAF, inflammatory CAF; myoCAF, myo-cancer-associated fibroblast; TE, transposable element
*Corresponding author.
**Corresponding author.
E-mail addresses: jxq1225@[Link] (X. Jiang), jcui@[Link] (J. Cui).
f
These authors contributed equally: Jinghan Wang, Meng Ren.

[Link] Vol 85 November, 2022 1

 Articles

Research in context
Evidence before this study HERVs at the single-cell level. There are substantial differences
Gallbladder cancer (GBC) is characterized by late-stage in HERV activation among different cell types, which
diagnosis and poor prognosis, so research efforts should contributes to the intratumoral heterogeneity of GBC. HERVs
continue to find biomarkers for early detection and innovate were demonstrated to serve as enhancers, potentially
therapeutic approaches for improving clinical outcomes. regulating the expression of neighboring genes in cancer cells.
Human endogenous retroviruses (HERVs) are remnants of Biological functions that may be affected by those aberrantly
ancient exogenous retroviruses and make up ∼8% of the activated HERVs were also identified for each cell population.
human genome. HERVs are specifically reactivated in various HERV members (HERVH and HHLA2) were recognized as
cancers and are implicated in cancer development, which has promising targets to help achieve early diagnosis and advance
spurred many studies exploring HERVs as promising GBC immunotherapy.
diagnostic and therapeutic targets. The involvement of HERVs
Implications of all the available evidence
in GBC, which is essential for improved management of this
Exploring HERV transcription patterns and their potential
malignancy, has not yet been systematically elucidated.
impacts on cellular functions highlights the essential role of
Added value of this study HERVs in shaping the tumor microenvironment, provides
We revealed the cell-type diversity of GBC by single-cell RNA novel insights into GBC development, and offers a valuable
sequencing and delineated the transcriptional landscape of resource for better management of GBC.

Introduction antigens that can serve as targets of T cells.13–17 These

Gallbladder cancer (GBC) is a rare malignancy and the characteristics make them especially attractive as new
most common type of GBC is adenocarcinoma.1,2 GBC therapeutic possibilities for diverse malignancies.
is often diagnosed at a late stage because it usually Actually, HERVs have recently received much attention
doesn’t cause noticeable signs or symptoms until it’s in immunotherapy and the combined use of DNA
advanced. The progression of GBC is frequently fast and demethylation agents (resulting in HERV derepression)
the prognosis is very poor. Current therapeutic ap- and immune checkpoint therapy is being tested in
proaches, such as surgical resection, targeted therapy clinical trials in a variety of cancers.18–20
and immunotherapy, benefit a subset of patients with Single-cell RNA sequencing (scRNA-seq) has been
GBC.3,4 Innovative strategies to achieve early diagnosis widely applied to dissect the cellular composition of
and effective treatment are urgently needed. tissues and explore functional heterogeneity.21–23 It pro-
Human endogenous retroviruses (HERVs), which vides a convenient way to investigate how HERV
are estimated to account for 8% of the human genome, expression varies among different cell types. Given that
are relics of ancestral germ-line infections by exogenous HERVs are dispersed in multiple copies throughout
retroviruses and became stable during the course of human genome, achieving accurate estimation of locus-
evolution.5,6 The intact structure of HERV is the same as specific expression of HERVs with scRNA-seq data is
exogenous retroviruses: 5′ LTR-gag-pro-pol-env-3′ LTR, challenging.24 Tools tailored for scRNA-seq data to ach-
but they usually accumulated so many mutations and ieve locus quantification taking multi-mapped reads into
deletions that their viral replication is restricted in the consideration are lacking. Keeping only uniquely map-
human genome.7 Although most HERVs are silenced in ped reads is a commonly used strategy in RNA-seq
normal tissues, cumulative evidence indicates that analysis, although it would lead to expressional under-
HERVs could be reactivated in diverse cancers due to estimation of some HERVs.10,25,26 Therefore, the findings
dysregulated epigenetic modifications and play an related to HERVs should be interpreted in the context of
important role in modulating carcinogenesis.8,9 The the techniques and methods used.
flanking LTRs could lead to expressional alteration of Given the essential role of HERVs in other cancers
adjacent genes by acting as promoters or enhancers. For and the fact that the contribution of HERVs to GBC is
example, the expression of Krüppel-associated box currently unknown, we would like to investigate
(KRAB) domain–containing zinc-finger protein (KZFP) whether HERVs are also reactivated and implicated in
genes has been demonstrated to be induced by adjacent various cellular processes in GBC, aiming to find new
HERVs in tumors.10 HERV expression may trigger a insights into GBC management. In this study, we
’viral mimicry’ pathway, which can increase antigen applied scRNA-seq to reveal the cellular composition of
presentation, expression of interferon-related genes, and the tumor microenvironment (TME) and report the
expression of immune checkpoint genes, potentially comprehensive single-cell transcriptome profile of
sensitizing tumor cells to recognition.11,12 Furthermore, HERVs in GBC tumor tissues and adjacent normal tis-
some HERVs may be translated to tumor-specific sues. In addition, we investigated the enhancer activity

2 [Link] Vol 85 November, 2022

 Articles

of HERV-derived sequences by dual-luciferase reporter Next, the unique molecular identifier (UMI) and cell
assay and explored their potential functional impacts in barcode were used to tag the synthetic cDNA at the 5′
each cell type. These results would help us further the end, that is equivalent to the 3′ end of the mRNA
understating of the functional role of HERVs and pro- molecule. The final single cell library was generated
vide clues to the discovery of better treatment strategies through several procedures including random priming
for GBC. and extension (RPE), RPE amplification PCR and WTA
index PCR. Library quantification was performed by
Methods Agilent Bioanalyzer 2200 High Sensitivity DNA Chip
Ethics approval and consent to participate and Qubit High Sensitivity DNA assay (Thermo Fisher
4 patients with gallbladder adenocarcinoma were Scientific). All generated libraries were sequenced by an
recruited and they all signed the consent forms. Tumor illumina sequencer (Illumina, San Diego, CA) with
tissue and adjacent normal tissue were collected from PE150 strategy (paired-end 150bp).
each GBC patient. The matched adjacent normal tissue
was taken from the mucosal tissue of the gallbladder at
least 2 cm from the edge of the tumor. This present Cell culture
study was approved by the Ethics Committee of Eastern The GBC-SD (RRID: CVCL_6903), human gallbladder
Hepatobiliary Surgery Hospital (EHBHKY2021-K-006). carcinoma cell line was purchased from Zrbiorise
(Shanghai, China). GBC-SD and HEK293T (ATCC,
RRID: CVCL_0063) cells were maintained in high-
Single-cell dissociation glucose DMEM (Gibco) supplemented with 10% fetal
Fresh tumor tissues and adjacent normal tissue bovine serum (FBS, Gibco), penicillin (100 IU/ml,
collected were stored in MACS Tissue Storage Solution Gibco) and streptomycin (100 μg/ml, Gibco) in a hu-
(Miltenyi Biotec) before processing. The single-cell midified atmosphere containing 5% CO2 at 37 ◦ C. The
suspension was generated as described below. First, GBC-SD and HEK293T cell lines were characterized by
the tissue sample was minced into small pieces that Azenta Life Sciences (Jiangsu, China) using short tan-
were ∼1 mm3 in size on ice after washing with dem repeat (STR) markers (Supplementary File 1).
phosphate-buffered saline (PBS, Gibco). Collagenase IV
(Worthington) and DNase I (Worthington) were subse-
quently used to enzymatically dissociated these pieces Plasmids
and this step lasted for 30 min at 37 ◦ C. After dissoci- PGL3-promoter, PGL3-control and pRL-TK were pur-
ation, the sample was passed through a 70 μm cell chased from HedgehogBio Science and Technology Ltd
strainer (Falcon) and the mixture was then centrifuged (Shanghai, China). All potential enhancer sequences of
for 5 min at 300×g to remove the supernatant. Next, the HERV (Supplementary Table S4) were amplified from
red blood cells were lysed with red blood cell lysis buffer the HEK293T genome using nested PCR with
(Miltenyi Biotec) and the sample was washed with PBS 2 × Phanta Max Master Mix (Vazyme, P525). The PGL3-
containing 0.04% BSA (Thermo Fisher Scientific). A promoter was enzymatically cut using KpnI-HF (NEB,
35 μm cell strainer (Falcon) was added to re-filtered cell R3142S). All potential enhancer sequences of HERV
pellets after re-suspending in PBS containing 0.04% were inserted into the PGL3-promoter separately using
BSA. Cell viability was assessed by staining dissociated NEBuilder HiFi DNA Assembly Master Mix (NEB,
cells through Calcein-AM (Thermo Fisher Scientific) E2621). The sequences of all plasmids were confirmed
and Draq7 (BD Biosciences). Finally, the dead cells in by Sanger sequencing.
the single cell suspension were removed using a MACS
dead cell removal kit from Miltenyi Biotec.
Dual-luciferase reporter assay
A total of 5 × 104 GBC-SD cells were plated in 24-well
Single-cell library preparation and sequencing cell culture plates (corning, 3524) and transfected with
Whole transcriptomic information of each sample was 500 ng plasmids per well by using Lipofectamine 2000
captured with the BD Rhapsody single-cell system. To Transfection Reagent (Invitrogen, 11668019). The ratio
realize single-cell capture, the suspension was randomly of the experimental vector to the co-reporter vector pRL-
distributed across more than 200,000 microwells by TK was 24:1. After 48 h, the cells were collected for
limited dilution method. Beads containing oligonucle- luciferase activity evaluation by using a dual-luciferase
otide barcodes were required to be added excessively to reporter system (HANBIO, China) according to the
make sure that almost each microwell has one bead. manufacturer’s protocol.
Trapped cells were lysed, allowing the released mRNA
molecules to hybridize with barcoded capture oligos on
the beads in the microwell. Beads were transferred to a Quantification of gene and HERV expressions
single tube and then reverse transcription and ExoI Raw sequencing data was analyzed using whole tran-
(Thermo Fisher Scientific) digestion were carried out. scriptome analysis (WTA) pipeline of BD Rhapsody™

[Link] Vol 85 November, 2022 3

 Articles

on a local installation (BD® Single-Cell Multiomics significant; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001;
Analysis Setup User Guide, Doc ID: 47383). The steps ****, p ≤ 0.0001.
in the WTA analysis mainly include removing reads
with low quality, annotating R1 and R2 reads, collapsing
Differential expression analysis
reads into raw molecules, determining putative cells and
The specific HERV families (adjusted p-value < 0.05 and
generating expression matrices. For details, please refer
log2FC > 0.5) of each major cell type from GBC tissues
to BD® Single-Cell Multiomics Bioinformatics Hand-
were selected by the function FindAllMarkers(). Differ-
book (Doc ID: 54169). During the annotation process,
entially expressed HERV families and HERV loci be-
STAR (v2.5.2b)27 was used to align reads to the human
tween tumor- and normal-derived cells were identified
reference genome (GRCh38.p12, UCSC). Gene anno-
using the function FindMarkers(). HERV families or
tation file was downloaded from GENCODE (https://
HERV loci were considered statistically significant if
[Link]/, Release 31) and HERV anno-
their adjusted p-values by Bonferroni were less than
tation file which was compiled by RepeatMasker was
0.05. The contribution of significantly upregulated
obtained from UCSC Table Browser for GRCh38
HERV loci to the total increment in HERV expression in
([Link] Only the
each cell type was calculated as the ratio of the sum of
uniquely mapped reads were calculated to estimate gene
expression increments from significantly upregulated
and HERV expression levels. Besides, many HERV loci
HERV loci over the sum of expression increments from
were found to overlap with exons of host genes, so the
all upregulated HERV loci (log2FC > 0). The expression
overlapping regions were removed for HERV loci to
increment of a HERV locus was defined as the average
avoid quantification bias.
change in expression between cells from the tumor and
adjacent normal tissue.
Quality control and cell type determination
The gene-cell count matrix was imported into R package
Seurat (version 4.0.6) for subsequent analyses.28 Low HERV-derived enhancer prediction in epithelial cells
quality cells and cell doublets or multiplets were To search for HERV loci that may have enhancer activity
removed based on the gene count and mitochondrial to improve the expression of adjacent differentially
contamination (Supplementary Table S1) and genes expressed genes (DEGs) (genes with adjusted p-value
expressed in less than 3 cells were filtered out for each < 0.05), pairs between HERV loci and DEGs within
sample. After applying these filtering criteria, the 500 kb in the genome were selected as the base set.29–31
filtered gene-cell count matrix was normalized by the The initial screening was based on whether HERVs and
function NormalizeData() and 2000 most variable genes DEGs were both upregulated (gene.log2FC > 0 &
of each sample were selected to correct the batch effect HERV.log2FC > 0) in GBC-derived epithelial cells. The
derived from individual samples. Significant principal further filtering was carried out according to these 3
components identified by function ElbowPlot() were aspects: co-expression between HERV locus and DEG,
used for graph-based clustering and t-distributed sto- enhancer-gene link predicted by GeneHancer, and
chastic neighbor embedding (tSNE) visualization. Sub- DNase and histone modification signal from ENCODE
clustering of cell types of interest was done with the project. The co-expression was defined if there was a
same method. The identity of each cluster was charac- positive correlation (Spearman correlation coefficient
terized by the expression of the canonical marker genes. > 0.3) between HERV and DEG expression. Gene-
For HERV(locus)-cell count matrix, only cells that Hancer interactions were downloaded from the GeneLoc
passed through the previous filter and HERV loci that database ([Link]
were expressed in at least 3 cells were retained and [Link]) and ENCODE Candidate Cis-Regulatory
HERV(family)-cell count matrix was calculated by Elements (cCREs) data was obtained from the UCSC
aggregating counts of each HERV locus belonging to the Table Browser ([Link]
same family. The filtered HERV-cell count matrices
were also normalized by the function NormalizeData().
Identification of malignant cells
There are two key considerations when separating ma-
Comparison of HERV expression lignant cells from non-malignant epithelial cells. On the
To describe the difference of HERV expression between one hand, copy number variation (CNV) is thought to be
GBC and normal samples, the number of active HERV a characteristic of malignant cells and inferCNV
loci and total expression level of HERV loci were ([Link] was used
calculated and displayed. We defined a HERV locus to detect somatic chromosomal copy number alter-
active if it was expressed in at least 3 cells for each ations. This algorithm was implemented for each pa-
sample. The statistical method used for each compari- tient and fibroblasts and endothelial were considered as
son was two-sided Wilcoxon rank sum test and signifi- the reference. The CNV signal of each epithelial cell was
cance levels were indicated by these symbols: ns, not summarized as CNV score, which was the mean square

4 [Link] Vol 85 November, 2022

 Articles

of the CNV estimates across all genomic locations. On T cells, natural killer (NK) or natural killer T (NKT) cells,
the other hand, phenotypic similarities of malignant monocytes, macrophages, dendritic cells (DCs),
cells will drive them to cluster together, so we per- neutrophils, mast cells, fibroblasts and endothelial cells
formed sub-clustering for epithelial cells with a high (Fig. 1b and c, Supplementary Figure S1a and b). T cells
resolution. A cell cluster composed mainly of cells and epithelial cells were the most abundant cells in
with high CNV scores was designated as malignant. both tumor and adjacent normal tissues (Fig. 1d,
Here, subcluster 11 and 12 showing highest CNV Supplementary Figure S1c). Moreover, epithelial cells,
scores were classified as malignant cells (Supplementary monocytes and macrophages were predominantly
Figure S3c). enriched in tumor tissues.
We characterized the HERV expression patterns in
patients with GBC and found that samples from GBC
Single-cell trajectory construction
tissues had more active HERV loci and higher HERV
The epithelial cell trajectory was generated by Monocle
expression levels than samples from adjacent normal
(Version 2.20.0) algorithm.32 Differentially expressed
tissues (Fig. 1e). The derepression of HERVs in tumor
genes selected by Seurat were used to define each cell’s
tissues implies that HERVs may play a role in the
progress and “DDRTree” method was applied to reduce
initiation and progression of GBC. Furthermore,
data dimensionality. HERV families that changed along
dimension reduction analyses based on HERV expres-
with the developmental trajectory were calculated by
sion alone showed that the cells with the same cell labels
function differentialGeneTest().
tended to cluster together, indicating that HERVs were
actively transcribed in tumors in a cell type-specific
Identification of genes and putative pathways manner (Fig. 1f). To identify which HERVs contribute
associated with HERV families to tumor heterogeneity in terms of cellular composition,
To investigate how dramatically upregulated HERV the differential expression of all HERV families and
families (adjusted p-value < 0.05, log2FC > 0.5) poten- HERV loci was calculated in each cell type (Fig. 1g,
tially influence the cellular function, we performed Supplementary File 3). Interestingly, HERVE, HERVK
correlation analysis to find the genes whose expression and HERVH, the most reported families in cancer
pattern was similar to that of HERV family (Spearman research, were found to be mainly expressed in epithe-
correlation coefficient > 0.3) and gene ontology (GO) lial cells (Supplementary Figure S1d).
enrichment analysis was used to identify pathways
enriched by top 100 correlated genes.
HERVs as enhancers potentially regulate adjacent
DEGs
HERV-gene interaction prediction For all cell types except NK/NKT cells, HERVs were
To predict the interactions between upregulated HERV extensively activated in tumors compared to adjacent
loci and neighboring DEGs (within 500 kb) in GBC- normal tissues and HERVs were expressed at the
derived cell types, we calculated the expression correla- highest level in epithelial cells (Fig. 2a, Supplementary
tions between them. HERV locus and DEG were Figure S2a and b). Based on the criteria of log2FC
thought to be associated if HERV expression was posi- > 0.5 and adjusted p-value < 0.05, only 4 significantly
tively correlated with DEG expression (Spearman increased HERV loci were detected in T cells derived
correlation coefficient > 0.3). from tumors compared to T cells derived from adjacent
normal tissues and only 5 HERV loci were signifi-
Role of funders cantly increased in GBC-derived DCs (Supplementary
The funders played no role in study design, data Figure S2c). Actually, for each cell type, the increment
collection, data analyses, interpretation, or writing of caused by significantly elevated HERVs loci accounted
report. for only a small part of the total increment caused by all
HERVs (Fig. 2b), suggesting that the HERV loci that
have not passed the cut-off should also be considered by
Results researchers. This also reminds us to explore the
Aberrant activation of HERVs in GBC abnormal expression of HERVs in GBC at both the lo-
To investigate the cellular heterogeneity and charac- cus level and family level.
terize the molecular signature in GBC, tumor tissues The exaptation of HERVs as regulatory elements
and matched normal tissues were collected from 4 pa- influencing the transcription of host genes has been
tients with gallbladder adenocarcinoma to perform established by many studies.25,34 To find HERVs which
scRNA-seq (Fig. 1a, Supplementary Table S2). After act as enhancers to drive the expression of the neigh-
quality filtering, a total of 28,301 cells were obtained and boring DEGs, scRNA-seq data for epithelial cells com-
catalogued into 11 main cell types annotated with ca- bined with other informative data were used to predict
nonical marker genes, including epithelial cells, B cells, the possible links (Fig. 2c). First, we picked out the pairs

[Link] Vol 85 November, 2022 5

 Articles

b Epithelial cellsc d 100 Epithelial cells

Cell type percentage(%)

B cells EPCAM 5 B cells
25 T cells CD79A 5 T cells
NK/NKT cells 75
CD19 5
tSNE_2

NK/NKT cells
Monocytes
0 CD3D 5 Monocytes
Macrophages 50
Dendritic cells CD3E 5 Macrophages
Neutrophils CD3G 5 Dendritic cells
−25
Mast cells CD2 6 25 Neutrophils
Fibroblasts KLRD1 5 Mast cells
Endothelial cells
−25 0 25 GNLY 7 0
Fibroblasts
tSNE_1 FCN1 5 Normal Tumor Endothelial cells
e VCAN 5
10000 10285 CD68 5 Epithelial cells Depletion
LGMN 5 B cells Enrichment
Number of active loci

T cells
7500 CD1C 5
6612
6224 6213 6332 6236 NK/NKT cells
CD1E 4 RO/E
5260 Monocytes
5000 4738 CLEC10A 4 Macrophages
0.75
FCGR3B 6 Dendritic cells 1.00
CXCR2 5 Neutrophils 1.25
2500
CMTM2 5 Mast cells 1.50
CSF3R 6 Fibroblasts
0 Endothelial cells
TPSAB1 7
N1 T1 N2 T2 N3 T3 N4 T4
MS4A2 5 Normal Tumor

**** COL1A1 8 f 50 Epithelial cells

THY1 5
B cells
HERV expression level

750 PECAM1 5 T cells

25
VWF 5 NK/NKT cells
tSNE_2

Monocytes
Epithelial cells

B cells

T cells

NK/NKT cells

Monocytes

Macrophages

Mast cells
Dendritic cells

Neutrophils

Fibroblasts

Endothelial cells

500 **** **** **** 0 Macrophages

Dendritic cells
−25 Neutrophils
250
Mast cells
Fibroblasts
−50 Endothelial cells
0 −50 −25 0 25 50
N1 T1 N2 T2 N3 T3 N4 T4
tSNE_1
g
Epithelial cells_Tumor

B cells_Tumor
Percent
T cells_Tumor Expressed
0
NK/NKT cells_Tumor 20
40
Monocytes_Tumor
60
Macrophages_Tumor Average
Expression
Dendritic cells_Tumor 2
1
Neutrophils_Tumor
0
Mast cells_Tumor −1

Fibroblasts_Tumor

Endothelial cells_Tumor
H VE nt
VH int
M ER nt
65 D
PA LT int
− 8

V2 M R5 t
4B ST 2D

LT t
rim int

LT 3
M 1G3
LT 1F1
ER F
LT 1A
ER B
M 1B

M G
M T1I
LT T1K

M A
M ER B
57 C
LT int
ER 2
ER R D
P1 s
M −int

M TR 8
41 A
M int
M T1B
E 4
ER ER D0

LT E− t
49 t
M −int
R 0
ER 2
L 7D
ST A

t
E n

M −in

− n
R in

in
S− 5−H
BL R5

M 2B

M R7

L R4

M T2B

LT R5
M 47B
M 10

M R2

16
M ST

ER 37

M TR6
ER 11

ER 1 1

U T 1
1
i

M −i

M A−i

VL 4−i
ER −
ER −

−P B−

A−
2

H L 4

M R4

5
LT

L
H 1D

LT
R

L
R

E
L
1
VK
ER
H

ER

Fig. 1: Single-cell landscape and HERV expression patterns in GBC and adjacent normal tissues. a Overview of the study design. b tSNE projection of 28,301 single
cells from tumor and adjacent normal tissues, color-coded by cell type. c Violin plot displaying expression of canonical marker genes for major cell types. d Relative
proportion of major cell types in different tissue types (top). Dot plot showing the distribution of major cell types in different tissue types estimated by Ro/e,33 the ratio of
observed to expected cell numbers in each cell type (bottom). e Comparison of the number of active HERV loci (top) and the overall expression of HERVs (bottom) in the
matched samples. Statistical significance was evaluated by two-sided Wilcoxon rank sum test. f tSNE plot of all single cells from tumor tissues only based on the HERV
expression. The top 1000 variable HERV loci were used. g Dot plot displaying expression of cell type-specific HERV families for major cell types.

6 [Link] Vol 85 November, 2022

 Articles

a b
1000 Epithelial cells 23.32%
B cells 32.22%
HERV Expression Level

**** **** **** ns **** **** *** **** **** **** ****
750 T cells 7.50%
Monocytes 9.91%
Macrophages 11.43%
500 Normal Dendritic cells 15.14%
Tumor Neutrophils 12.40%

250 Mast cells 16.28%

Fibroblasts 16.70%
Endothelial cells 5.79%
0
ls

ls

es

es

ils

ls

ts

ls
ls
lls

ls
el

el

el
l

as
ce

ce

ph
ce

ce

yt

ag

tc

lc
lc

oc

bl
tro
ph
KT

ic
B

ia
ia

as

ro
on

it

eu

el
el

ro

b
dr
N

M
M

Fi

th
ith

ac
K/

N
en

do
Ep

En

c Co-expression
Gene expression

Co-upregulation HERV expression

GeneHancer Candidate
HERV enhancer
HERV
Gene HERV Gene

Gene

Normal Tumor
ENCODE cCREs
Epigenomic signal

Relative luciferase activity

20
*** **
d e
15

PGL3-promoter SV40 promoter Fluc Poly A 10

***
PGL3-control SV40 promoter Fluc Poly A SV40 enhancer 5
**** ***
*** ** ***
1
0
PGL3-promoter-potential
Candidate enhancer SV40 promoter Fluc Poly A
1

1F du 8-c 5

-c 9
L t-d -c 1
R - d 6- r

20
- d 3- 6
M C-d p7 r16
R 2- 5- 16

enhancer derived from HERV

LT 1B up cto

LT A p2 r

87 r1
LT 3- p hr

8 1

D p1 hr1
13 u ch

hr
-in p3 r
LT TR up hr

u p ch
h

10 du ch
M R3 -du 3-c
1 -d e

ST u -c
ER 1 l v
M R4A tro

pRL-TK HSV TK promoter Rluc Poly A

n
M Co
E

Fig. 2: HERV-derived enhancers prediction. a Comparison of overall expression of HERVs between GBC and normal tissue cells for each major cell type. Statistical
significance was evaluated by two-sided Wilcoxon rank sum test. b Contribution of significantly upregulated HERV loci to the increase in HERV expression. c Workflow of
HERV-derived enhancers prediction. d Schematic of dual-luciferase reporter assay. The pGL3-Promoter Vector contained an SV40 promoter upstream of the luciferase
gene. HERV containing putative enhancer elements can be inserted upstream of the promoter-luc+ transcriptional unit. PGL3-control containing SV40 enhancer se-
quences was used as a positive control. e Dual-luciferase reporter assay to test enhancer activity of the genomic fragments derived from HERV in GBC-SD. Data are shown
as the means ± SD. **, p < 0.01; ***, p < 0.001; ****, p < 0.0001. Significance values for each comparison were calculated by Student’s t-test. Use Welch’s t-test when two
groups don’t have the same variance. At least three biological repeats were carried out.

of HERVs and neighboring DEGs in which both HERV expression between HERV and neighboring DEG in
loci and neighboring DEGs were upregulated in GBC-derived epithelial cells; 2) predefined enhancer-
epithelial cells from GBC. Then, further filtering was gene link by GeneHancer database which predicted
performed based on the following criteria: 1) co- the associations between regulatory elements

[Link] Vol 85 November, 2022 7

 Articles

a 50 b 0.020 c EDN1 d 100 Epi_EDN1(low)

Depletion

Cell Type Percentage(%)

Enrichment
4 Epi_EDN1(high)

CNV Score
0.015

Expression Level
25 75 Malignant
tSNE_2

3 Epi_EDN1
0.010 (low)
50 RO/E
0 2
0.005 Epi_EDN1 0.8
1 25 1.0
(high)
−25
0.000 1.2
0 0
Malignant

nt
w

gh
−40 −20 0 20 Epi_EDN1 Epi_EDN1

na
Normal Tumor

lo

hi
1(

ig
tSNE_1

1(
(low) (high)

al
N
ED

M
ED
Epi_EDN1(low) Normal Tumor

i_

i_
Ep
Epi_EDN1(high)

Ep
Epi_EDN1(low) Epi_EDN1(high) Malignant
Malignant g
e 1000 f Malignant vs Non−malignant Average
1

Component 2
*
HERV Expression Level

Expression 0
****
**** 1.0 −1
750 Malignant
0.5 −2
0.0
−3
500
−0.5 0 5 10
Epi_EDN1 Component 1
(high) Percent
250
1

Component 2
Expressed
20
Epi_EDN1 0
0 40
(low)
60 −1
)

nt
w

gh

na
lo

hi
1(

ig
1(

ER 7Y

nt

M 2
LO 4D

ST t

B
−2
N

in

in
al

1N

41
N

−i
ED

1−

A−
M

LT

ER
ED

VH

R
LT
LT

R
i_

LT
i_

−3
Ep

M
Ep

5 0 5 10 0 5 10 0 5 10
HHLA2
h Epi_EDN1(low) Epi_EDN1(high) Malignant j CD274 3 Component 1

HERVH−int PDCD1LG2 2 k Malignant

100.0 CD80 3 LTR7 LTR7Y

10.0 4 R = 0.4, p < 2.2e−16 R = 0.38, p < 2.2e−16

CD86 7.5
1.0
VTCN1 4
0.1
Relative Expression

5.0
LTR7 TNFRSF14 4
100.0 2.5
ICOSLG 4
10.0
PVR 4 0.0
1.0
0.1 4 HERVH−int LTR7B
CD40
LTR7Y 3 R = 0.35, p = 5.9e−15 R = 0.32, p = 3.8e−12
CD70 7.5
100.0
CD200 3
10.0
5.0
1.0 CD276 4
0.1 4 2.5
IDO1
0 5 10 15
Pseudo−time 0.0
)
)

gh
w

0 1 2 0 1 2
i
lo

hi

nt
1(

1(

na

HHLA2
N

ig
ED

ED

al
M
i_

i_
Ep

Ep
protein targeting

Epi_EDN1(low) Epi_EDN1(high) Malignant

SEC61G 0.72 0.68 0.59 0.59 0.23 0.21 −0.14 0.18
RPS6 0.7 0.66 0.55 0.58 0.2 0.2 −0.15 0.14
SEC61G
to ER

RPS28 0.7 0.68 0.53 0.58 0.23 0.21 −0.2 0.19

RPS6
RPS12 0.65 0.61 0.49 0.53 0.19 0.19 −0.15 0.11
0.1 RPS28
RPS24 0.63 0.59 0.46 0.54 0.16 0.19 −0.14
0.1 RPS12
RPS29 0.61 0.58 0.46 0.49 0.16 0.2 −0.2
UBA52 0.59 0.55 0.5 0.45 0.15 0.16 −0.11 0.14 RPS24
CCL20 0.68 0.68 0.59 0.62 0.21 0.21 −0.21 0.24 RPS29
CSF1 0.63 0.62 0.54 0.56 0.24 0.25 −0.14 0.25 UBA52
chemotaxis cell proliferation cell differentiation
granulocyte mononuclear mononuclear

CCL3L1 0.62 0.6 0.54 0.58 0.22 0.22 −0.14 0.2 CCL20
CXCL1 0.61 0.59 0.51 0.49 0.22 0.13 −0.15 0.19 Corr CSF1
CXCL11 0.61 0.59 0.51 0.55 0.22 0.19 −0.14 0.17 1.0 CCL3L1 Expression
PPIA 0.59 0.55 0.46 0.45 0.15 0.13 −0.03 0.1 CXCL1
2
CXCL3 0.58 0.58 0.49 0.55 0.21 0.2 −0.14 0.18 0.5 CXCL11
CXCL10 0.52 0.57 0.47 0.55 0.22 0.15 −0.11 0.25 PPIA 1
CCL3 0.54 0.54 0.48 0.54 0.23 0.23 −0.14 0.14 0.0
CXCL3 0
CD70 0.72 0.72 0.61 0.57 0.25 0.22 −0.2 0.28
−0.5 CXCL10 −1
CD24 0.69 0.71 0.59 0.58 0.21 0.17 −0.21 0.19 CCL3
TNFSF9 0.67 0.67 0.59 0.59 0.18 0.19 −0.18 0.21
−1.0 CD70
−2
TWSG1 0.64 0.61 0.55 0.52 0.25 0.2 −0.12 0.23 CD24
RIPK2 0.63 0.67 0.56 0.57 0.24 0.18 −0.12 0.26 TNFSF9
CD320 0.6 0.6 0.51 0.5 0.24 0.21 −0.17 0.23
TWSG1
GAL 0.6 0.58 0.49 0.56 0.22 0.18 −0.12 0.18
RIPK2
TMEM176A 0.65 0.64 0.61 0.51 0.22 0.2 −0.17 0.22
CD320
NFKBIZ 0.58 0.56 0.54 0.49 0.24 0.14 −0.14 0.23
GAL
SOX13 0.66 0.64 0.53 0.56 0.28 0.18 −0.16 0.25
TMEM176A
DUSP10 0.57 0.6 0.52 0.5 0.26 0.24 −0.19 0.21
TMEM176B 0.49 0.51 0.5 0.42 0.19 0.16 −0.15 0.21 NFKBIZ
IFI16 0.53 0.52 0.5 0.48 0.25 0.17 −0.13 0.2 SOX13
DUSP10
nt

B
ER 7 Y

LO 4D

ST t

TMEM176B
in

in
R

1N

41
−i

1−

A−
R

LT

ER

IFI16
VH

R
LT
LT

LT
M
M

M
H

Fig. 3: Identification of malignant cells and HERVH derepression in the process of malignant transformation. a tSNE projection of
epithelial cells, color-coded by cell type. b Boxplot showing the CNV score for each epithelial cell subtype. c Violin plot displaying expression of
EDN1 in two non-malignant subtypes. d Relative proportion of epithelial cell subtypes in different tissue types (left). Dot plot showing the

8 [Link] Vol 85 November, 2022

 Articles

(enhancers and promoters) and target genes based on transformation confirmed these assumptions (Fig. 3g,
multiple information sources35; 3) HERV region Supplementary Figure S3h). Tumor and adjacent
showing enhancer-like signature supported by high normal tissues contained all subtypes at different frac-
DNase and H3K27ac with low H3K4me3 signal from tions (Fig. 3d, Supplementary Figure S3e). As expected,
ENCODE cCREs.36 As long as pairs of HERVs and malignant cells were dominantly identified in tumors,
neighboring DEGs met any two of the above criteria, whereas Epi_EDN1(low) cells were more enriched in
HERVs were considered to be candidate enhancers adjacent normal tissues.
modeling the expression of the neighboring DEGs. The We compared the total expression level of HERVs
predicted HERV-DEG pairs were provided in among these subtypes and found that malignant cells
Supplementary File 4. exhibited the highest level (Fig. 3e). HERVs showed a
We selected some HERVs predicted by our pipeline higher signal of transcription in Epi_EDN1(high) than
to validate that the genomic fragments derived from in Epi_EDN1(low), indicating that the activation of
HERVs had enhancer activity using dual-luciferase re- HERVs has already occurred in the epithelial cells at
porter assay (Fig. 2d and e, Supplementary Table S3). intermediate state before transforming into malignant
Interestingly, non-LTR sequence derived from HERV, cells. Next, we identified HERV families that showed
such as MLT1F-int-dup5-chr16 (chr16:29915600- differential expression in these subtypes. Strikingly,
29916253), also showed significant enhancer activity. HERVH-associated elements (HERVH-int, LTR7Y)
The neighboring genes, such as UCA1, SEPHS2 and upregulated in Epi_EDN1(high) (vs. Epi_EDN1(low))
DCTPP1, that are potentially regulated by HERV- were further elevated in malignant cells (Fig. 3f,
derived enhancers, have been reported to play onco- Supplementary Figure S3f, Supplementary File 5). At
genic roles in tumor proliferation and metastasis the locus level, 37/55 of the upregulated HERV loci in
(Supplementary Table S4).37–40 However, some neigh- malignant cells (vs. non-malignant cells) belong to the
boring genes, TMEM219 and YPEL3, mediate anti- HERVH family (Supplementary Figure S3g,
tumor activities.41,42 Supplementary File 5). These results imply the func-
tional importance of HERVH in GBC formation. To
HERVH expression is increased with the malignant further explore the aberrant activation of HERVH, we
transformation of epithelial cells constructed a transcriptional trajectory with the defini-
GBC originates from epithelial cells, which has been tive malignant and non-malignant epithelial cells
demonstrated by previous studies.43,44 To distinguish (Fig. 3g, Supplementary Figure S3h). Among the
malignant and normal epithelial cells resident in GBC HERVs displaying transcriptional alterations with the
tissues, large-scale CNVs were inferred with stromal tumor progression, HERVH-int, LTR7Y and LTR7 of
cells as references. The epithelial cells were divided into the HERVH family showed the highest significance
subclusters and the subclusters with markedly higher (Fig. 3h, Supplementary File 6), suggesting that
CNV scores were identified as malignant cells (Fig. 3a HERVH may be a promising biomarker for early GBC
and b, Supplementary Figure S3a–c). The remaining diagnosis. To elucidate the impacts of HERVH dere-
non-malignant subclusters were further annotated as pression on this biological process, we identified the
Epi_EDN1(low) and Epi_EDN1(high) according to the genes whose expression was correlated with the
expression level of EDN1 (Fig. 3a and c). The Epi_ expression of these HERVH elements in malignant
EDN1(high) subtype, which simultaneously exhibited cells. These genes were highly enriched for GO terms
high expression of a mesenchymal marker (MMP7) and related to protein targeting to ER, granulocyte chemo-
a cancer stem cell marker (CD44), was thought to taxis, cell proliferation and differentiation (Fig. 3i,
represent an intermediate cell state (ICS) between Supplementary Figure S3i). It is possible that the alter-
normal and malignant cells (Supplementary Figure S3d) ation of these genes’ expression was caused by the
and the Epi_EDN1 (low) subtype was expected to derepression of HERVH. HERVH has been reported to
represent normal epithelial cells. The Epi_EDN1 (high) play a crucial role in pluripotency maintenance in hu-
cells at the middle stage and Epi_EDN1 (low) cells at the man pluripotent stem cells25,45 and functions of HERVH
early stage in the transcriptional trajectory of malignant in malignant cells revealed by us are consistent with

distribution of epithelial cell subtypes in different tissue types estimated by Ro/e (right). e Comparison of overall expression of HERVs among
epithelial cell subtypes. Statistical significance was evaluated by two-sided Wilcoxon rank sum test. f HERV families upregulated in malignant
cells compared to non-malignant cells (adjusted p-value < 0.05 and log2FC > 0.5). g Differentiation trajectory of malignant and normal
epithelial cells inferred by Monocle2, color-coded by cell type. h Dynamic expression changes of HERVH elements along the differentiation
trajectory. i Heatmaps depicting expression correlation between gene and HERV family in malignant cells (left, Spearman correlation coefficient
was also shown) and relative expression of these correlating genes across epithelial cell subtypes (right). j Violin plot displaying expression of the
collected immune checkpoints in epithelial cell subtypes. k Scatter plot showing expression correlation between HHLA2 and HERVH element.
Spearman correlation coefficient and two-tailed p-value were shown.

[Link] Vol 85 November, 2022 9

 Articles

a c 100
Naïve CD4+ T

Cell type percentage(%)

Depletion
Naïve CD4+ T Naïve CD4+ T CD4+ Th
25 75
Enrichment
CD4+ Th

tSNE_2
CD4+ Th
Regulatory/Exhausted CD4+ T Regulatory/Exhausted CD4+ T
Regulatory/Exhausted CD4+ T RO/E
0 CD8+ T GZMK
50 CD8+ T GZMK CD8+ T GZMK
CD8+ T GZMB 0.50
Proliferative T CD8+ T GZMB
CD8+ T GZMB 0.75
−25 CD4−/CD8− T 25 Proliferative T
Proliferative T 1.00
CD4−/CD8− T
1.25
CD4−/CD8− T
−20 0 20 40 0 1.50
tSNE_1 Normal Tumor Normal Tumor

b CD3D
5

CD3E
5
d Normal
5
CD3G 600 Tumor
6
CD2
5
ns ns ** **** ** * ** ** **

HERV Expression Level

CD4
4
LEF1
5
400
CCR7
5
SELL
4
CD40LG
5
KLRB1 200
6
IL7R
5
FOXP3
5
TIGIT
5 0
CTLA4

Th

T
5

ZM

ZM
CD8A

4+

4+

ive

8−
4+
D

D
G

t
ra
D
C

/C
5

ife
C
GZMK
ve

ed

4−
8+

8+

ol
aï

st

D
Pr
D

D
6
N

C
au

C
GZMB

xh
/E
4

y
MKI67

or
at
ul
+T

Th

eT

−T

eg
CD ory/

ZM

ZM
T

R
ativ
4+
D4

D8
4+
ted lat

TG

TG
CD
C

/C
fer
us egu
ïve

4−
8+

8+

oli
R

CD
Na

Pr
CD

CD
ha
Ex

e Proliferative T
Regulatory/ Percent
Exhausted CD4+ T Expressed CD8+ T GZMB Percent
Expressed
CD4−/CD8− T
Percent
Expressed
Percent 5 7 4
Expressed 6 8 8
20 7 9 12
Tumor 40 Tumor Tumor Tumor
8 10 16
60 9 11
Average 10 12 Average
Expression Average Average Expression
Normal 0.4 Normal Expression Normal Expression Normal 0.4
0.0 0.4 0.4 0.0
−0.4 0.0 0.0 −0.4
HERVH−int −0.4 MER65−int

nt
−0.4
H TR Y
VH A
t
R B
M 2C
1D

7Y
ER i n

−i
L 7
ER 13

LT 41
M −
R

1
LT

in
LT

LT

qu
le
ar
H

al r
orm mo
_N _Tu
f tive
T
tive
T
era era
Pr olif Pr olif
Expression
CD24 0.62 0.2 0.54 0.33 0.28 0.01 Corr CD24 2
1.0
IL1B 0.57 0.39 0.52 0.4 0.42 0.13 IL1B
T cell activation

VSIG4 0.49 0.35 0.5 0.42 0.37 0.25 0.5 VSIG4 1

CD86 0.47 0.35 0.49 0.42 0.4 0.25 CD86
0.0 0
TNFSF13B 0.46 0.29 0.49 0.3 0.35 0.21 TNFSF13B
0.2 −0.5 CD4
CD4 0.48 0.29 0.49 0.34 0.33
FGL2 0.48 0.36 0.49 0.36 0.33 0.21 −1.0 FGL2
GPR183 0.56 0.34 0.56 0.4 0.39 0.13 GPR183
CLEC7A 0.49 0.32 0.49 0.35 0.43 0.22 CLEC7A
IL1A 0.51 0.37 0.42 0.38 0.36 0.11 IL1A
LILRB4 0.5 0.33 0.48 0.44 0.37 0.3 LILRB4
ICAM1 0.5 0.23 0.42 0.38 0.31 0.19 ICAM1
MMP9 0.53 0.31 0.54 0.41 0.45 0.13 MMP9
nuetrophil activation

OLR1 0.52 0.32 0.49 0.47 0.41 0.27 OLR1

CYFIP1 0.5 0.32 0.46 0.4 0.39 0.28 CYFIP1
PLAUR 0.49 0.32 0.5 0.44 0.38 0.2 PLAUR
RAB31 0.45 0.32 0.49 0.4 0.43 0.28 RAB31
PLAU 0.39 0.37 0.41 0.41 0.41 0.26 PLAU
CD14 0.43 0.33 0.48 0.38 0.43 0.3 CD14
FCGR2B 0.46 0.32 0.46 0.37 0.42 0.16 FCGR2B
CTSB 0.43 0.29 0.42 0.37 0.4 0.34 CTSB
IL18 0.44 0.37 0.43 0.37 0.39 0.28 IL18
CREG1 0.44 0.3 0.44 0.38 0.39 0.31 CREG1
LYZ 0.36 0.23 0.37 0.31 0.38 0.33 LYZ
7Y

nt

D
13

41

12

1
−i
R

LT
VH
R

ER

R
LT

M
LT

LT
ER

M
H

Fig. 4: Proliferative T cells showing prominent expression levels of HERVs. a tSNE projection of T cells, color-coded by cell type. b Violin plot
displaying expression of canonical marker genes for T cell subtypes. c Relative proportion of T cell subtypes in different tissue types (left). Dot
plot showing the distribution of T cell subtypes in different tissue types estimated by Ro/e (right). d Comparison of overall expression of HERVs
between GBC and normal tissue cells for each T cell subtype. Statistical significance was evaluated by two-sided Wilcoxon rank sum test. e HERV

10 [Link] Vol 85 November, 2022

 Articles

that. In addition, we noticed that the gene UCA1 HERVs were expressed at different levels in different
showed a strong association with HERVH-int and was subtypes of T cells. When we compared cells from tu-
specifically increased in malignant cells (Supplementary mor tissues with those from adjacent normal tissues,
Figure S3j), which fit well with the fact that UCA1, significant HERV transcription changes were observed
which has been found to enhance the proliferation, in all T cell subtypes except naïve CD4+ T and CD4+ Th
migration and invasion of bladder cancer cell, is a cells (Fig. 4d). Proliferative T cells were characterized by
lncRNA generated by HERVH elements.46 the prominent expression levels of HERVs in our data.
HERVs are highly valued in immunotherapy, such as HERVH-associated elements (LTR7Y and HERVH-int),
immune checkpoint inhibitors, for their ability to sensi- LTR13A, MER41B, LTR12C and MLT1D were upregu-
tize tumor cells to immunological recognition.47 Human lated in proliferative T cells from tumor tissues (Fig. 4e).
endogenous retrovirus-H long terminal repeat- The top correlated genes with upregulated HERVH
associating 2 (HHLA2), a recently emerging immune were associated with T cell activation (Fig. 4f,
checkpoint, is considered to be derived by HERVH.48 Supplementary Figure S4d). Strikingly, we observed that
HHLA2 was thought to play a dual role in carcinogen- many genes involved in neutrophil behaviors, such as
esis: one for immunostimulation and another for neutrophil activation and degranulation, were tran-
immunosuppression.49 We examined the transcription of scriptionally increased in tumor-derived proliferative T
HHLA2 in epithelial cell subtypes, as well as some im- cells and displayed association with HERVH, MER41B
mune checkpoints potentially expressed in tumor cells and LTR12C. The association between neutrophil ac-
(Fig. 3j). We found that HHLA2 expression was detect- tivity and one HERVK locus has been reported in pe-
able in Epi_EDN1(high) (represent ICS) and malignant ripheral blood mononuclear cells (PBMCs) from elderly
cells and was higher than that of PD-L1 (CD274), which people.50 Here, we hypothesized that the neutrophil
is the most studied immune checkpoint. The weak signal state is regulated by proliferative T cells in GBC and that
for PD-L1 transcription may explain, at least in part, the HERV may play an important role in this process.
limited efficacy of immunotherapy targeting PD-L1 in Besides, high expression of HERVH-int in regulato-
some patients with GBC. The elevated expression level of ry/exhausted CD4+ T cells, MER65-int in CD8+ T
HHLA2 in intermediate and malignant cells and its dual GZMB and LTR7Y and Harlequin-int in CD4−/CD8− T
role (one for immunostimulation and another for cells was also detected in tumor tissues (Fig. 4e).
immunosuppression) in carcinogenesis suggest that
HHLA2 may serve as a favorable candidate for GBC
treatment. Further, we found an expressional correlation
The association of IgG genes with the MLT1C locus
in plasma B cells
of HHLA2 and HERVH elements (Fig. 3k), implying that
B cells that infiltrate the TME play a multifaceted role in
combining DNA demethylating agents that aim to induce
modulating the tumor immunity.51 Three major B cell
HERVH activation with targeting HHLA2 therapy may
subtypes, including follicular B cells (MS4A1), plasma B
be a promising therapeutic strategy.
cells (IGHG1) and granzyme B-secreting B cells
(GrB+ B cells, GZMB), were identified from our data
Proliferative T cells show abnormally high based on the signature genes (Fig. 5a and b,
expression of HERVs Supplementary Figure S5a and b). The composition of
T cells display heterogeneity in cellular composition and these B cell subtypes was distinct between GBC and
functional states in the TME. Here, T cells were further normal samples (Fig. 5c, Supplementary Figure S5c).
partitioned into 7 distinct subpopulations annotated by The relative proportion of plasma B cells was observed
marker genes, including naïve CD4+ T, CD4+ T helper to be increased in GBC, whereas follicular B cells were
(CD4+ Th), regulatory/exhausted CD4+ T, CD8+ T enriched in normal samples and comprised the majority
GZMK, CD8+ T GZMB, proliferative T and CD4−/ of B cells in the adjacent normal tissues.
CD8− T cells (Fig. 4a and b, Supplementary Figure S4a Direct comparison of the expression level of HERVs
and b). Compared with the adjacent normal tissues, in GBC and normal tissues for each B cell subtype
tumor tissues showed reduced proportions of naïve revealed that tumor-derived follicular B and plasma B
CD4+ T and CD4−/CD8− T cells and increased pro- cells exhibited the increased abundance of HERV tran-
portions of proliferative T and regulatory/exhausted scription (Fig. 5d). Compared with non-tumor-derived
CD4+ T cells (Fig. 4c, Supplementary Figure S4c). More follicular B cell, cells in tumors showed higher expres-
regulatory/exhausted CD4+ T cells accumulated in tu- sion levels of HERVH-int and MSTC (Fig. 5e). Plasma B
mor tissues, indicating a change from immune activa- cells, which are terminally differentiated B cells, can
tion to immune suppression during tumor progression. secrete antibodies that are an essential component of

families upregulated in GBC-derived T cell subtypes (adjusted p-value < 0.05 and log2FC > 0.5). f Heatmaps depicting expression correlation
between gene and HERV family in tumor-derived proliferative T cells (left, Spearman correlation coefficient was also shown) and relative
expression of these correlating genes in all proliferative T cells (right).

[Link] Vol 85 November, 2022 11

 Articles

a Follicular B b 6 d Normal
MS4A1

HERV expression level

Plasma B
400 **** *** ns Tumor
20 GrB+ B
tSNE_2

300
0 9
IGHG1
200
−20
100
6
GZMB
−40 0
−20 −10 0 10 20

B
B

B
ar

ar

+
a
tSNE_1

sm

rB
+
sm
ul

ul
rB
lic

G
lic

a
a

G
l

Pl
Pl

l
Fo

Fo
c e Follicular B Average
100 Follicular B Depletion Expression
Cell type percentage(%)

Plasma B Follicular B Enrichment

0.4
75 GrB+ B Tumor
RO/E 0.0
1 −0.4
50 Plasma B
2
Percent
25
3 Normal Expressed
5.0
GrB+ B
7.5
0 10.0
Normal Tumor Normal Tumor HERVH−int MSTC

Plasma B_Tumor Plasma B_Tumor Plasma B_Tumor

f 8
R = 0.63 p < 2.2e−16 R = 0.55 p < 2.2e−16 R = 0.43 p = 4.5e−11
7.5 7.5
6
IGHG1

IGHG4

IGHG3
5.0 5.0
4

2.5 2 2.5

0.0 0 0.0
0 2 4 6 0 2 4 6 0 2 4 6
MLT1C−dup586−chr14 MLT1C−dup586−chr14 MLT1C−dup586−chr14

g Plasma B Percent
h 9

IGHG1
Expressed
0
8
5
Tumor 10 IGHG4
15
20 8
Average
IGHG3
Expression
Normal 0.4
al

or
m

0.0
m
or

Tu
N

−0.4
B_
B_

a
a

m
m

MLT1C−dup586−chr14
as
as

Pl
Pl

Fig. 5: IgG genes are associated with the neighboring MLT1C locus in plasma B cells. a tSNE projection of B cells, color-coded by cell type. b
Violin plot displaying expression of canonical marker genes for B cell subtypes. c Relative proportion of B cell subtypes in different tissue types
(left). Dot plot showing the distribution of B cell subtypes in different tissue types estimated by Ro/e (right). d Comparison of overall expression
of HERVs between GBC and normal tissue cells for each B cell subtype. Statistical significance was evaluated by two-sided Wilcoxon rank sum
test. e HERV families upregulated in GBC-derived Follicular B cells (adjusted p-value < 0.05 and log2FC > 0.5). f Scatter plot showing expression
correlation between the indicated gene and HERV locus. Spearman correlation coefficient and two-tailed p-value were shown. g Dot plot
displaying expression level of the MLT1C locus in plasma B cells. h Violin plot displaying expression of IgG genes in plasma B cells.

humoral immunity.52 Immunoglobulin G (IgG), a type IgG genes (IGHG1, IGHG4 and IGHG3) (Fig. 5f). The
of antibody, is produced and released by plasma B cells. expression association and upregulation of both them
Through our analysis, we discovered a transcriptional led us to speculate that MLT1C-dup586-chr14 may be a
association between one MLT1C locus (MLT1C-dup586- regulator that can promote the expression of these
chr14, chr14:105750174-105750614) and neighboring neighboring IgG genes (Fig. 5g and h).

12 [Link] Vol 85 November, 2022

 Articles

Overexpression of HERVH in all myeloid cell types tumor migration and invasion through various mecha-
Emerging evidence emphasizes the key role of myeloid nisms, such as remodelling the extracellular matrix
cells in modulating cancer progression.53 In our study, (ECM) and modulating the tumor immune system.57
the myeloid cells were separated into monocytes, mac- Sub-clustering of fibroblasts revealed 2 main distinct
rophages, DCs and neutrophils (Fig. 1b). DCs were subtypes, including inflammatory CAFs (iCAFs;
further categorized into conventional DC (cDC1 and FDGFRA) and myo-cancer-associated fibroblasts (myo-
cDC2), monocyte-derived DC (mo–DC) and mature DC CAFs; RGS5) (Fig. 7a and b, Supplementary Figure S7a
based on the prominent expression markers (Fig. 6a and and b). In both GBC tissues and adjacent normal tis-
b, Supplementary Figure S6a and b). For both tumor sues, iCAFs were the predominant cell type. Compared
and adjacent tissues, cDC2 and mo-DC accounted for to that in normal tissues, the proportion of myoCAFs in
the majority of DCs. Notably, mo-DCs and mature DCs GBC tissues was increased, whereas the level of iCAFs
were enriched in GBC tissues, while cDC1s and cDC2s was slightly decreased (Fig. 7c, Supplementary
were depleted (Fig. 6c, Supplementary Figure S6c). Figure S7c).
HERV expression patterns were different among DC Regarding the expression level of HERVs, both
subtypes (Fig. 6d). cDC1s and mo-DCs in GBC showed iCAFs and myoCAFs in tumors presented significant
significant elevations in HERV expression, but no sig- elevation (Fig. 7d). HERVs that were abnormally acti-
nificant differences were detected in cDC2s and mature vated in GBC were screened (Fig. 7e, Supplementary
DCs. Mo-DCs are believed to arise from monocytes in File 5). For iCAFs, MER65-int, MER4D0 and MER4-
the context of inflammation or infection and this sub- int are the top upregulated HERVs. At the locus level,
population can induce T cell activation in various tumor MER65-int-dup5-chr19 (chr19:41729515-41730138) was
models.54,55 We found that LTR7Y were overexpressed in associated with its neighboring gene CEACAM5 in
mo-DCs derived from GBC tissues (Fig. 6e). Genes iCAFs (Fig. 7f). CEACAM5, encoding carcinoembryonic
presenting similar expression patterns with LTR7Y in antigen (CEA), has been used as a tumor biomarker in
mo-DCs were mainly implicated in antiviral response and clinical detection.58 After checking their locations, we
T cell activation (Fig. 6f, Supplementary Figure S6d), found that part of the MER65-int-dup5-chr19 sequence
suggesting that the enhancement of mo-DCs’ antiviral serves as the exon of CEACAM5. Besides, this associa-
and inflammatory abilities may be regulated by LTR7Y. tion was also observed in GBC-derived macrophages
What’s more, LTR7Y exhibited strong association with (Fig. 6i). Whether the transcription of this MER65-int
gene MMP7 (Fig. 6g), which functions as an oncogenic element induces the high expression of CEACAM5 in
factor to mediate occurrence and progression of several iCAFs and macrophages from GBC tissues deserves
types of cancers.56 This association was also seen in ma- further investigation. For myoCAFs, MER65-int specif-
lignant cells, proliferative T cells and macrophages, ically activated in GBC-derived cells exhibited an asso-
which further supports the hypothesis that elevated ciation with genes related to ECM organization (Fig. 7g,
MMP7 level in these cells may be associated with the Supplementary Figure S7d). This suggests the possi-
activation of LTR7Y (Supplementary Figure S6e). bility that active MER65-int may be involved in dysre-
The significant increases in expression of HERVH gulation of ECM genes.
elements were also observed in other myeloid cells, In endothelial cells from GBC tissues, MER50,
monocytes (HERVH-int, LTR7C and LTR7), macro- LTR7C, LTR6A and MLT1A0 were upregulated
phages (HERVH-int and LTR7Y) and neutrophils (Fig. 7h). However, how these HERV families affect
(LTR7C) (Fig. 6h). Genes associated with HERVH ele- cellular functions is not clear.
ments in tumor-derived monocytes and macrophages
were mostly involved in immune cell activation and
differentiation (Supplementary Figure S6d). These bio- Discussion
logical processes play key roles in TME. Here, we GBC might not be detected until it’s advanced and
speculated that these profound perturbations occurred catching the cancer early will benefit patient survival.
in myeloid cells might be related to the abnormal tran- The prognosis of GBC is poor and only some patients
scription of HERVH elements. Taken together, we exhibited desired outcome. Therefore, the pursuit for
discovered that HERVH elements were upregulated in early diagnostic biomarkers and more effective treat-
almost all myeloid cell types and they may have similar ment strategies for GBC is ongoing. HERVs have
impacts on myeloid cell behavior and function. received a lot of attention because of their strong asso-
ciation with cancer development and immunotherapy.
MER65-int is activated in fibroblasts and has an In this study, we determined the cellular composition
expressional correlation with ECM genes and presented a comprehensive single-cell transcrip-
Cancer-associated fibroblasts (CAFs) are one of the most tional profile of HERVs for GBC tissues, highlighting
dominant components of the tumor stroma and have that HERVs were transcribed in a cell type-specific
heterogeneous phenotypes and functions. A large manner. We found that HERVs, as enhancers, have
number of studies support that CAFs can promote the potential to alter host gene expression and further

[Link] Vol 85 November, 2022 13

 Articles

a b CLEC9A
3 c
20 cDC1 3
XCR1 100
cDC2 cDC1 Depletion

Cell type percentage(%)

mo−DC 5 cDC2 cDC1
CD1C Enrichment
10 75
Mature DC mo−DC
tSNE_2

5 cDC2 RO/E
CD14 Mature DC
0 50 0.8
6
C1QB mo−DC 1.0
25
C1QC
6 1.2
−10 Mature DC 1.4
5 0
CCR7
Normal Tumor Normal Tumor
−20 −10 0 10 20 30 LAMP3
4

tSNE_1
1

C
C

D
cD

cD

o−

e
ur
m

at
mo−DC_Normal

M
mo−DC_Tumor
d 500 * ns ** ns Normal
f
HERV expression level

Tumor PTPN22 0.57 PTPN22

400
IFI6 0.51 IFI6
300 defense response STAT1 0.49 STAT1
to virus

200 PARP9 0.49 PARP9

AIM2 0.48 AIM2
100
IFI27 0.47 IFI27 Expression
0 Corr
BECN1 0.46 BECN1 2
1.0
1

C
C

TBK1 0.45 TBK1 1

cD

cD

o−

0.5
e
ur
m

at

TNFSF9 0.62 TNFSF9 0

0.0
M

regulation of T cell activation

e mo−DC Average
Expression
CD24 0.57
−0.5 CD24 −1

IL1A 0.53 −1.0 IL1A

0.4
PRDM1 0.52 PRDM1
0.0
Tumor
−0.4 NCSTN 0.52 NCSTN

Percent LAG3 0.48 LAG3

Expressed PIK3R1 0.48 PIK3R1
15
Normal 20 CD47 0.47 CD47
25
30 CASP3 0.47 CASP3
35 LTR7Y
LTR7Y 40
g 2.5
mo−DC_Tumor MMP7
i Macrophages_Tumor Macrophages
Percent Expressed
2.5
R = 0.71, p < 2.2e−16 R = 0.38, p < 2.2e−16 2.5
2.0 2.0 3 5.0
Expression Level

Tumor 7.5
CEACAM5

1.5 1.5 10.0

MMP7

2
1.0 1.0 Average Expression

0.5 0.5 1 0.4

Normal
0.0
0.0 0.0
0 −0.4
0 2 4 6 mo-DC mo-DC
Tumor
0 2 4 6
LTR7Y Normal CEACAM5
MER65−int−dup5−chr19

h Monocytes Average Macrophages Average Neutrophils Average

Expression Expression Expression

0.4 0.4 0.4

Tumor Tumor 0.0 Tumor 0.0
0.0
−0.4 −0.4 −0.4

Percent Percent Percent

Expressed Expressed Expressed
Normal Normal 10 Normal 5
20
40 20 10
60 30 15
40 20
VL 9B

TH int

M C
LT C
7C

A
VH int
R t
LT 7C
TH 9B

LT A
TH R7

R D
C

H R4B int
V int
R nt
M R65 1A

M 6A
LT TB
LT 7Y
TH 9B
LT E1A

LT B
2
LT −in

50
ST
E1
ST
E1

LT E1
12

1E
M 48

−
LT H−i
ER 5−

ER TR

R
R

E −
ER −
E E

R
R
S

−E

M
R
M T1
H R6

L
L
E

M
M

Fig. 6: The activation of HERVH in myeloid cells. a tSNE projection of DCs, color-coded by cell type. b Violin plot displaying expression of canonical marker genes for DC
subtypes. c Relative proportion of DC subtypes in different tissue types (left). Dot plot showing the distribution of DC subtypes in different tissue types estimated by Ro/e
(right). d Comparison of overall expression of HERVs between GBC and normal tissue cells for each DC subtype. Statistical significance was evaluated by two-sided
Wilcoxon rank sum test. e HERV families upregulated in GBC-derived mo-DCs (adjusted p-value < 0.05 and log2FC > 0.5). f Heatmaps depicting expression correla-
tion between gene and HERV family in tumor-derived mo-DCs (left, Spearman correlation coefficient was also shown) and relative expression of these correlating genes in
all mo-DCs (right). g Scatter plot showing expression correlation between the indicated gene and HERV family (left). Spearman correlation coefficient and two-tailed p-
value were shown. Violin plot displaying expression of gene MMP7 in mo-DCs (right). h HERV families upregulated in GBC-derived myeloid cell subtypes (adjusted p-value
< 0.05 and log2FC > 0.5). i Scatter plot showing expression correlation between the indicated gene and HERV locus (left). Spearman correlation coefficient and two-tailed
p-value were shown. Dot plot displaying expression of gene CEACAM5 in macrophages (right).

14 [Link] Vol 85 November, 2022

 Articles

a iCAF
b c
PDGFRA 5
20 myoCAF 100 iCAF Depletion

Cell type percentage(%)

myoCAF Enrichment
75 iCAF
tSNE_2

0 RO/E

RGS5 6 50 0.9
1.0
−20 25 myoCAF 1.1
1.2
0 1.3
−20 0 20 iCAF myoCAF Normal Tumor Normal Tumor
tSNE_1
d e iCAF Percent myoCAF Percent
Normal Expressed Expressed
10 10
Tumor 20 20
HERV expression level

400 **** **** 30 Tumor 30

Tumor
40
Average
50
Expression
200 Average
Expression 0.4
Normal Normal
0.0
0.4
−0.4
0 0.0
iCAF myoCAF −0.4

t
ER nt

10 int
M E1B int
H ER t
ER 0

R C
VL B

39

ER int

nt

2C

0
in
LT −in
M −in

in
M 4D

1A
ST
ER 4
M 5−i

LT R12

−i
1−
ER −
TH 37−

B−
ER

1
VH

65
4

LT
R
M
6

34
M

LT

M
ER

ER

ER
M
H
M

M
f g myoCAF_Normal myoCAF_Tumor
iCAF_Tumor
SFRP2 0.07 −0.05 −0.05 0.47 0.12 0.22 0.1
ITGA11 0.21 −0.02 0.01 0.41 0.16 0.17 0.1 SFRP2
R = 0.3, p < 2.2e−16
extracellular matrix organization

MMP2 0.04 −0.07 −0.16 0.39 0.02 0.12 0.14

MFAP5 −0.05 0.04 −0.04 0.38 0.08 0.22 0.18 ITGA11
2
Corr
COL6A3 −0.02 −0.02 −0.03 0.38 0.12 0.14 0.21 1.0
POSTN 0.04 −0.02 0.1 0.37 0.13 0.15 0.21 MMP2
CEACAM5

0.5
COL1A1 0.1 0.05 0.21 0.37 0.12 0.08 0.14
MFAP5
COL11A1 0.04 0.04 0.08 0.35 0.13 0.19 0.17 0.0
1 COL3A1 0.08 0.07 0.18 0.34 0.03 0.08 0.19 COL6A3
−0.5 Expression
MMP9 0.09 0.04 0.05 0.34 0.05 0.14 0.08
COL1A2 0.14 0.05 0.21 0.33 0 0.11 0.15 −1.0 POSTN 2

COL6A2 −0.02 −0.02 −0.02 0.33 −0.06 0.01 0.09 1

COL1A1
0 FN1 0 −0.02 0.1 0.31 0.02 0.16 0.07 0

0 2 4 6 COL6A1 −0.02 −0.03 0.08 0.3 −0.05 0 0.01 −1

COL11A1
MER65−int−dup5−chr19 −2
39

nt

t
C

0
n

in

1A
ST

12
−i

−i

B−
ER

COL3A1
VH

65

LT
R
M

34
M

LT

CEACAM5
ER

M
ER

ER

h
M
H

MMP9
M

3 Endothelial cells Percent COL1A2

Expressed
10 COL6A2
20
Expression Level

2
30 FN1
Tumor
40
COL6A1
1
Average
Expression

0 0.4
Normal
iCAF_Normal iCAF_Tumor 0.0

−0.4
6A
7C
50

0
1A
ER

R
R

LT

LT
LT
M

Fig. 7: The activation of MER65-int in fibroblasts. a tSNE projection of fibroblasts, color-coded by cell type. b Violin plot displaying expression of canonical marker genes
for fibroblast subtypes. c Relative proportion of fibroblast subtypes in different tissue types (left). Dot plot showing the distribution of fibroblast subtypes in different
tissue types estimated by Ro/e (right). d Comparison of overall expression of HERVs between GBC and normal tissue cells for each fibroblast subtype. Statistical sig-
nificance was evaluated by two-sided Wilcoxon rank sum test. e HERV families upregulated in GBC-derived fibroblast subtypes (adjusted p-value < 0.05 and log2FC > 0.5).
f Scatter plot showing expression correlation between the indicated gene and HERV locus (top). Spearman correlation coefficient and two-tailed p-value were shown.
Violin plot displaying expression of the indicated gene in iCAFs (bottom). g Heatmaps depicting expression correlation between gene and HERV family in tumor-derived
myoCAFs (left, Spearman correlation coefficient was also shown) and relative expression of these correlating genes in all myoCAFs (right). h HERV families upregulated in
GBC-derived endothelial cells (adjusted p-value < 0.05 and log2FC > 0.5).

change the characteristics of tumors. We also suggested strategy for GBC treatment. Furthermore, for each cell
that HERVH may be a candidate for early GBC diag- type, we explored which biological functions are poten-
nosis and targeting HHLA2 might be an appealing tially affected by various HERV families. Our study

[Link] Vol 85 November, 2022 15

 Articles

could provide a framework for future discoveries of transcriptional information of each HERV locus is
functional HERVs as molecular and cellular therapeutic lost.66,67 A number of computational tools, such as
targets for GBC. RepEnrich, TEtranscripts, REdiscoverTE and so on, are
Here, dual-luciferase reporter assay was applied to produced to quantify transposable element (TE, HERV
demonstrate that some HERVs can act as enhancers. belongs to TE) expression for bulk RNA-seq data and the
HERVs also work as alternative promoters, resulting in best choice of the strategy should be guided by the spe-
transcriptional initiation of host genes, including cific biological question.65,68–70 Nevertheless, current re-
cancer-related genes.59 In addition to acting as tran- sources available for TE quantification at single-cell
scriptional regulators, HERV generate noncoding RNA resolution are relatively limited, although scRNA-seq
(ncRNA) or protein product to modulate the gene reg- technology opens the possibility to investigate TE tran-
ulatory network, as a perpetrator or protector in carci- scription variability among different cell populations,
nogenesis.5,60 These functions were not explained in our factors driving such diversity and cellular phenotypes
study, but deserve to be further explored. Within each influenced by TE derepression.24 Pioneering pipelines,
cell population, we identified biological functions that including scTE71 and a framework with assembled tran-
may be biased by the aberrantly expressed HERVs, but scripts,72 have been generated to report the TE expression
whether and by which mechanisms distinct HERVs at single-cell resolution, but they both have certain limi-
affect cellular behavior remain to be investigated. tations to application, namely counting only at the family
The widespread reactivation of HERVs has been level and counting only assembled transcripts respec-
linked to some malignancies, suggesting their potential tively. In this study, we considered using uniquely map-
to help cancer screen.61 In our study, the expression of ping reads to achieve TE quantification to ensure
HERVH was progressively elevated with malignant mapping accuracy and obtain locus information. As the
transformation of epithelial cells and displayed the biological significance of HERVs is becoming under-
highest significance compared to other HERV families, stood, we believe that there will be a rapid advance in
leading us to speculate that HERVH may be an early scRNA-seq computational pipelines tailored for HERV to
indicator of GBC. The role of HERVH as biomarker has be compatible with a wide variety of scRNA-seq protocols.
been discussed in several cancers, including colorectal In summary, our work highlights the functional role
carcinoma, prostate cancer and lung cancer.62–64 Further of HERVs in GBC and provides a new resource for
experiments, at the nucleic acid or protein level, are cancer diagnosis and management.
needed to verify whether HERVH is an effective diag-
nostic biomarker. HHLA2, a newly discovered immune Contributors
checkpoint belonging to B7 family, is thought to be J.W., X.J. and J.C. conceived the research. J.C. directed the analysis. J.W.
and M.R. performed the analysis. J.W. and M.R. contributed equally.
derived by HERVH.48 The expression of HHLA2 was J.Y., M.H., X.W. and W.M. conducted the experiments. All the authors
detected in a subset of intermediate and malignant cells, discussed the results and wrote the paper. J.C., J.W. and M.R. have
which is higher than that of some hot immune check- verified the underlying data. All authors read and approved the final
points such as PD-L1. Considering the fact that there are version of the manuscript.
substantial differences in expression of immune
Data sharing statement
checkpoints among diverse cancers and even in the
Raw single-cell RNA sequencing data reported in this study are available
same type of cancer, the expression profile of immune in Genome Sequence Archive for Human (GSA-Human) under acces-
checkpoints varies from patient to patient, although the sion number HRA001917.
expression of HHLA2 was modest in our data, it might
be an ideal immunotherapeutic target for a specific
Declaration of interests
group of GBC patients. The family diversity of HERVs The authors declare that they have no competing interests.
contributes to their functional complexity in cells. The
discovery of the association of HHLA2 with HERVH
indicates that epigenetic agents specially tailored to Acknowledgements
target HERVH should be a research priority in com- This study was supported by the National Natural Science Foundation of
China (31970176, 81972256) and the research grants from the Innova-
bined epigenetic and immune therapy. tion Capacity Building Project of Jiangsu province (BM2020019).
The fact that HERVs are present in multiple copies
within the human genome poses a challenge for quan- Appendix A. Supplementary data
tification.65 A common strategy used in routine analysis is Supplementary data related to this article can be found at [Link]
to keep only reads that uniquely map to the HERV loci. org/10.1016/[Link].2022.104319.
The advantage of this approach is that the expressional
References
signal for each locus can be obtained, but the disadvan- 1 Canale M, Monti M, Rapposelli IG, et al. Molecular targets and
tage is that it tends to underestimate the transcript level emerging therapies for advanced gallbladder cancer. Cancers
of HERVs, especially evolutionarily young HERVs. (Basel). 2021;13(22):5671.
2 Baiu I, Visser B. Gallbladder cancer. JAMA. 2018;320(12):1294.
Another approach, using the “multi-mapper” strategy, 3 De Lorenzo S, Garajova I, Stefanini B, Tovoli F. Targeted therapies
largely preserves the HERV-derived reads, but the for gallbladder cancer: an overview of agents in preclinical and

16 [Link] Vol 85 November, 2022

 Articles

clinical development. Expert Opin Investig Drugs. 2021;30(7): 28 Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating
759–772. single-cell transcriptomic data across different conditions, tech-
4 Song X, Hu Y, Li Y, Shao R, Liu F, Liu Y. Overview of current nologies, and species. Nat Biotechnol. 2018;36(5):411–420.
targeted therapy in gallbladder cancer. Signal Transduct Target Ther. 29 van Arensbergen J, van Steensel B, Bussemaker HJ. In search of
2020;5(1):230. the determinants of enhancer-promoter interaction specificity.
5 Zhang M, Liang JQ, Zheng S. Expressional activation and func- Trends Cell Biol. 2014;24(11):695–702.
tional roles of human endogenous retroviruses in cancers. Rev Med 30 Marbach D, Lamparter D, Quon G, Kellis M, Kutalik Z,
Virol. 2019;29(2):e2025. Bergmann S. Tissue-specific regulatory circuits reveal variable
6 Lander ES, Linton LM, Birren B, et al. Initial sequencing and modular perturbations across complex diseases. Nat Methods.
analysis of the human genome. Nature. 2001;409(6822):860–921. 2016;13(4):366–370.
7 Alcazer V, Bonaventura P, Depil S. Human endogenous retrovi- 31 Corces MR, Granja JM, Shams S, et al. The chromatin accessibility
ruses (HERVs): shaping the innate immune response in cancers. landscape of primary human cancers. Science.
Cancers. 2020;12(3):610. 2018;362(6413):eaav1898.
8 Petrizzo A, Ragone C, Cavalluzzo B, et al. Human endogenous 32 Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves
retrovirus reactivation: implications for cancer immunotherapy. complex single-cell trajectories. Nat Methods. 2017;14(10):979–982.
Cancers. 2021;13(9):1999. 33 Guo X, Zhang Y, Zheng L, et al. Global characterization of T cells in
9 Steiner MC, Marston JL, Iniguez LP, et al. Locus-specific charac- non-small-cell lung cancer by single-cell sequencing. Nat Med.
terization of human endogenous retrovirus expression in prostate, 2018;24(7):978–985.
breast, and colon cancers. Cancer Res. 2021;81(13):3449–3460. 34 Chuong EB, Elde NC, Feschotte C. Regulatory evolution of innate
10 Ito J, Kimura I, Soper A, et al. Endogenous retroviruses drive KRAB immunity through co-option of endogenous retroviruses. Science.
zinc-finger protein family expression for tumor suppression. Sci 2016;351(6277):1083–1087.
Adv. 2020;6(43):eabc3020. 35 Fishilevich S, Nudel R, Rappaport N, et al. GeneHancer: genome-
11 Brocks D, Schmidt CR, Daskalakis M, et al. Erratum: DNMT and wide integration of enhancers and target genes in GeneCards.
HDAC inhibitors induce cryptic transcription start sites encoded in Database (Oxford). 2017;2017:bax028.
long terminal repeats. Nat Genet. 2017;49(11):1661. 36 Consortium EP, Moore JE, Purcaro MJ, et al. Expanded encyclo-
12 Wolff F, Leisch M, Greil R, Risch A, Pleyer L. The double-edged paedias of DNA elements in the human and mouse genomes.
sword of (re)expression of genes by hypomethylating agents: from Nature. 2020;583(7818):699–710.
viral mimicry to exploitation as priming agents for targeted immune 37 Xuan W, Yu H, Zhang X, Song D. Crosstalk between the lncRNA
checkpoint modulation. Cell Commun Signal. 2017;15(1):13. UCA1 and microRNAs in cancer. FEBS Lett. 2019;593(15):1901–
13 Wang-Johanning F, Radvanyi L, Rycaj K, et al. Human endogenous 1914.
retrovirus K triggers an antigen-specific immune response in breast 38 Carlisle AE, Lee N, Matthew-Onabanjo AN, et al. Selenium detox-
cancer patients. Cancer Res. 2008;68(14):5869–5877. ification is required for cancer-cell survival. Nat Metab.
14 Schiavetti F, Thonnard J, Colau D, Boon T, Coulie PG. A human 2020;2(7):603–611.
endogenous retroviral sequence encoding an antigen recognized on 39 Nunziata C, Polo A, Sorice A, et al. Structural analysis of human
melanoma by cytolytic T lymphocytes. Cancer Res. 2002;62(19): SEPHS2 protein, a selenocysteine machinery component, over-
5510–5516. expressed in triple negative breast cancer. Sci Rep. 2019;9(1):16131.
15 Mullins CS, Linnebacher M. Endogenous retrovirus sequences as a 40 Lu J, Dong W, He H, et al. Autophagy induced by overexpression of
novel class of tumor-specific antigens: an example of HERV-H env DCTPP1 promotes tumor progression and predicts poor clinical
encoding strong CTL epitopes. Cancer Immunol Immunother. outcome in prostate cancer. Int J Biol Macromol. 2018;118(Pt
2012;61(7):1093–1100. A):599–609.
16 Smith CC, Beckermann KE, Bortone DS, et al. Endogenous retro- 41 Cai Q, Dozmorov M, Oh Y. IGFBP-3/IGFBP-3 receptor system as
viral signatures predict immunotherapy response in clear cell renal an anti-tumor and anti-metastatic signaling in cancer. Cells.
cell carcinoma. J Clin Invest. 2018;128(11):4804–4820. 2020;9(5):1261.
17 Saini SK, Orskov AD, Bjerregaard AM, et al. Human endogenous 42 Zhang J, Wen X, Ren XY, et al. Correction to: YPEL3 suppresses
retroviruses form a reservoir of T cell targets in hematological epithelial-mesenchymal transition and metastasis of nasopharyn-
cancers. Nat Commun. 2020;11(1):5660. geal carcinoma cells through the Wnt/beta-catenin signaling
18 Attermann AS, Bjerregaard AM, Saini SK, Gronbaek K, pathway. J Exp Clin Cancer Res. 2021;40(1):400.
Hadrup SR. Human endogenous retroviruses and their implication 43 Hundal R, Shaffer EA. Gallbladder cancer: epidemiology and
for immunotherapeutics of cancer. Ann Oncol. 2018;29(11):2183– outcome. Clin Epidemiol. 2014;6:99–109.
2191. 44 Misra S, Chaturvedi A, Misra NC, Sharma ID. Carcinoma of the
19 Chiappinelli KB, Strissel PL, Desrichard A, et al. Inhibiting DNA gallbladder. Lancet Oncol. 2003;4(3):167–176.
methylation causes an interferon response in cancer via dsRNA 45 Fort A, Hashimoto K, Yamada D, et al. Deep transcriptome
including endogenous retroviruses. Cell. 2015;162(5):974–986. profiling of mammalian stem cells supports a regulatory role for
20 Goel S, DeCristo MJ, Watt AC, et al. CDK4/6 inhibition triggers retrotransposons in pluripotency maintenance. Nat Genet.
anti-tumour immunity. Nature. 2017;548(7668):471–475. 2014;46(6):558–566.
21 Perez RK, Gordon MG, Subramaniam M, et al. Single-cell RNA-seq 46 Wang F, Li X, Xie X, Zhao L, Chen W. UCA1, a non-protein-coding
reveals cell type-specific molecular and genetic associations to RNA up-regulated in bladder carcinoma and embryo, influencing
lupus. Science. 2022;376(6589):eabf1970. cell growth and promoting invasion. FEBS Lett. 2008;582(13):1919–
22 Wang X, Miao J, Wang S, et al. Single-cell RNA-seq reveals the 1927.
genesis and heterogeneity of tumor microenvironment in pancre- 47 Wang L, Amoozgar Z, Huang J, et al. Decitabine enhances
atic undifferentiated carcinoma with osteoclast-like giant-cells. Mol lymphocyte migration and function and synergizes with CTLA-4
Cancer. 2022;21(1):133. blockade in a murine ovarian cancer model. Cancer Immunol Res.
23 Hornburg M, Desbois M, Lu S, et al. Single-cell dissection of 2015;3(9):1030–1041.
cellular components and interactions shaping the tumor immune 48 Mager DL, Hunter DG, Schertzer M, Freeman JD. Endogenous
phenotypes in ovarian cancer. Cancer Cell. 2021;39(7):928–944.e6. retroviruses provide the primary polyadenylation signal for two
24 O’Neill K, Brocks D, Hammell MG. Mobile genomics: tools and new human genes (HHLA2 and HHLA3). Genomics.
techniques for tackling transposons. Philos Trans R Soc Lond B Biol 1999;59(3):255–263.
Sci. 2020;375(1795):20190345. 49 Ying H, Xu J, Zhang X, Liang T, Bai X. Human endogenous
25 Wang J, Xie G, Singh M, et al. Primate-specific endogenous retrovirus-H long terminal repeat-associating 2: the next
retrovirus-driven transcription defines naive-like stem cells. Nature. immune checkpoint for antitumour therapy. EBioMedicine.
2014;516(7531):405–409. 2022;79:103987.
26 Goke J, Lu X, Chan YS, et al. Dynamic transcription of distinct 50 Autio A, Nevalainen T, Mishra BH, Jylha M, Flinck H, Hurme M.
classes of endogenous retroviral elements marks specific pop- Effect of aging on the transcriptomic changes associated with the
ulations of early human embryonic cells. Cell Stem Cell. expression of the HERV-K (HML-2) provirus at 1q22. Immun
2015;16(2):135–141. Ageing. 2020;17:11.
27 Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal 51 Sarvaria A, Madrigal JA, Saudemont A. B cell regulation in cancer
RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. and anti-tumor immunity. Cell Mol Immunol. 2017;14(8):662–674.

[Link] Vol 85 November, 2022 17

 Articles

52 Nutt SL, Hodgkin PD, Tarlinton DM, Corcoran LM. The generation 63 Manca MA, Solinas T, Simula ER, et al. HERV-K and HERV-H env
of antibody-secreting plasma cells. Nat Rev Immunol. proteins induce a humoral response in prostate cancer patients.
2015;15(3):160–171. Pathogens. 2022;11(1):95.
53 Engblom C, Pfirschke C, Pittet MJ. The role of myeloid cells in 64 Zare M, Mostafaei S, Ahmadi A, et al. Human endogenous retro-
cancer therapies. Nat Rev Cancer. 2016;16(7):447–462. virus env genes: potential blood biomarkers in lung cancer. Microb
54 Kuhn S, Yang J, Ronchese F. Monocyte-derived dendritic cells are Pathog. 2018;115:189–193.
essential for CD8(+) T cell activation and antitumor responses after 65 Lanciano S, Cristofari G. Measuring and interpreting transposable
local immunotherapy. Front Immunol. 2015;6:584. element expression. Nat Rev Genet. 2020;21(12):721–736.
55 Chow KV, Lew AM, Sutherland RM, Zhan Y. Monocyte-derived 66 Treangen TJ, Salzberg SL. Repetitive DNA and next-generation
dendritic cells promote Th polarization, whereas conventional sequencing: computational challenges and solutions. Nat Rev
dendritic cells promote Th proliferation. J Immunol. Genet. 2011;13(1):36–46.
2016;196(2):624–636. 67 Goerner-Potvin P, Bourque G. Computational tools to un-
56 Liao HY, Da CM, Liao B, Zhang HH. Roles of matrix mask transposable elements. Nat Rev Genet.
metalloproteinase-7 (MMP-7) in cancer. Clin Biochem. 2021;92:9–18. 2018;19(11):688–704.
57 Sahai E, Astsaturov I, Cukierman E, et al. A framework for 68 Criscione SW, Zhang Y, Thompson W, Sedivy JM, Neretti N.
advancing our understanding of cancer-associated fibroblasts. Nat Transcriptional landscape of repetitive elements in normal and
Rev Cancer. 2020;20(3):174–186. cancer human cells. BMC Genomics. 2014;15:583.
58 Berg KCG, Eide PW, Eilertsen IA, et al. Multi-omics of 34 colorectal 69 Jin Y, Tam OH, Paniagua E, Hammell M. TEtranscripts: a package
cancer cell lines - a resource for biomedical studies. Mol Cancer. for including transposable elements in differential expression
2017;16(1):116. analysis of RNA-seq datasets. Bioinformatics. 2015;31(22):3593–
59 Lamprecht B, Walter K, Kreher S, et al. Derepression of an endoge- 3599.
nous long terminal repeat activates the CSF1R proto-oncogene in 70 Kong Y, Rose CM, Cass AA, et al. Transposable element expression
human lymphoma. Nat Med. 2010;16(5):571–579, 1p following 9. in tumors is associated with immune infiltration and increased
60 Bannert N, Hofmann H, Block A, Hohn O. HERVs new role in antigenicity. Nat Commun. 2019;10(1):5228.
cancer: from accused perpetrators to cheerful protectors. Front 71 He J, Babarinde IA, Sun L, et al. Identifying transposable element
Microbiol. 2018;9:178. expression dynamics and heterogeneity during development at the
61 Gao Y, Yu XF, Chen T. Human endogenous retroviruses in cancer: single-cell level with a processing pipeline scTE. Nat Commun.
expression, regulation and function. Oncol Lett. 2021;21(2):121. 2021;12(1):1456.
62 Perot P, Mullins CS, Naville M, et al. Expression of young HERV-H 72 Shao W, Wang T. Transcript assembly improves expression quan-
loci in the course of colorectal carcinoma and correlation with tification of transposable elements in single-cell RNA-seq data.
molecular subtypes. Oncotarget. 2015;6(37):40095–40111. Genome Res. 2021;31(1):88–100.

18 [Link] Vol 85 November, 2022