# Local Raman Spectroscopy of Chiral Majorana Edge Modes in Kitaev Spin
Liquids and Topological Superconductors

James Jun He1, Naoto Nagaosa1,2 1 RIKEN Center for Emergent Matter Science
(CEMS), Wako, Saitama 351-0198, Japan
2 Department of Applied Physics, The University of Tokyo, Tokyo 113-8656,
Japan

###### Abstract

The Raman scattering with local optical excitation from the Majorana edge
modes of Kitaev spin liquids and topological superconductors is studied
theoretically. Although the effective one-dimensional model is common between
these two cases, the coupling to the electromagnetic field is different. It is
found that the Raman spectrum at low energy scales with $\omega^{3}$ in Kitaev
spin liquids while it shows the gap in topological superconductors. This is in
sharp contrast to the infrared absorption, where the spectrum shows the gap in
Kitaev spin liquids, while it behaves as $\sim\omega^{2}$ in topological
superconductors. This indicates that the electrodynamics of Majorana edge
modes depends on their higher-dimensional origins. The realistic estimate of
the Raman scattering intensity is given for $\alpha$-RuCl3 as the candidate
for Kitaev spin liquid.

Majorana fermions in condensed matter systems Read ; Kitaev2001 ; Fu2008 ;
Wilczek ; Vic2009  attract intensive interests recently. They are considered
as a promising candidate for the robust quantum computation with the non-
Abelian statistics Ivanov ; Alicea2011 ; Fujimoto ; SCZhangPNAS ; Sato ;
Fujimoto ; Nomura . Robustness comes from the no or weak coupling with the
surrounding systems. Majorana fermions are neutral particles, and hence is
intuitively expected to be decoupled from the electromagnetic field, which
makes its observation difficult. For example, propagating Majorana modes are
believed to exist in very different systems including quantum spin liquids
(QSLs) Kitaev ; Matsuda  and topological superconductors (TSCs) Fu2009 ;
Beenakker2009 ; Balents ; Qi2010 ; JJH_CP2019 ; Qinglin2017 ; ZWang ; XGWen .
In QSLs, the Majorana fermions are neutral since they come from spin operators
while their neutrality in TSCs is manifested by their equal weight of
electrons and holes at the zero energy. The origin of neutrality is that the
Majorana field operator $c(x)$ satisfies the relation
$c^{\dagger}(x)=e^{i\theta}c(x)$, with $\theta$ being some phase, and hence
the density operator $c^{\dagger}(x)c(x)=e^{i\theta}c(x)^{2}$ is a numeric
constant.

However, it has been shown in a previous study JJH2021  that the optical
response of Majorana edge modes in TSCs is nonzero. This is because of the
existence of Cooper pair condensate. Then, it raises the question whether
optical response also exists for the Majorana edge modes in QSLs or not. An
analytically solvable model of QSLs is given by Kitaev Kitaev  which supports
Majorana edge modes. There are two types of excitations in this model:
vortices and Majorana fermions. The Raman scattering due to both of them in
the bulk has been investigated Knolle1 ; Nasu ; Sandilands  while the infrared
absorption is associated with both the vortices and Majorana fermions Bolens .
But the optical response of the edge states remains an open question.
Particularly, optical microscopy and spectroscopy methods Jarzembski ;
Verma2017 ; Zhang2017  have been successful in detecting local optical
responses with high spatial resolution (up to a few nanometers Lee ;
Klingsporn ), which is ideal for the investigation of topological edge states
ZXShen ; Lai .

In this paper, we theoretically investigate the Raman scattering in a Kitaev
QSL with a magnetic field and in a spinless $p$-wave SC, both of which host
gapless chiral Majorana edge channels while the bulk energy spectrum is
gapped. In the QSL, we found that the edge channel induces Raman scattering
with an intensity proportional to $\omega^{3}$, $\omega$ being the photon
frequency shift. We estimate the realistic value of the intensity and find it
detectable with current local optical techniques. In contrast, the edge
channel in the $p$-wave SC does not induce Raman scattering unless
$\hbar\omega$ is larger than the bulk gap. Such essential differences between
our results for QSLs and TSCs indicate that there does not exist a general
effective edge theory which is sufficient to describe the electrodynamics of
Majorana edge modes of various origins.

## Quantum spin liquid —

Consider the following Kitaev model Hamiltonian Kitaev

$\displaystyle
H=\sum_{i,\nu}J_{\nu}\sigma_{\bm{r}_{i}}^{\nu}\sigma_{\bm{r}_{i}+\bm{d}_{\nu}}^{\nu}+h_{\nu}\sigma_{\bm{r}_{i}}^{\nu},$
(1)

where $\bm{r}_{i}$ is the displacement vector at site $i$ on a honeycomb
lattice as shown in FIG. 1(a). Each site has three nearest neighbors connected
by vectors $\bm{d}_{\nu=1,2,3}$ respectively and $\sigma^{\nu}_{\bm{r}_{i}}$
are the spin operators at site $i$. When the external magnetic field $h_{\nu}$
is small, the Hamiltonian for the low-energy sector can be written in the
Majorana representation Kitaev ; XT

$\displaystyle H=\frac{iJ_{\nu}}{2}\sum_{\langle
ij\rangle}c_{i}c_{j}+\frac{iK}{2}\sum_{\langle ij\rangle^{\prime}}c_{i}c_{j}.$
(2)

Here $\langle ij\rangle$ denotes nearest neighbours connected by
$\bm{d}_{1,2,3}$ while $\langle ij\rangle^{\prime}$ are next-nearest
neighbours as indicated by the dashed arrows in FIG. 1(a). $K$ is determined
by $h_{\nu}$. Note that we have excluded the vortex excitations in Eq. (2)
because they are gapped and we consider the energy regime inside the gap.
Thus, our results are applicable when the temperature is small compared to the
vortex gap and the photon frequency $\omega$ is smaller than this gap.

When $h_{\nu}=0$ and thus $K=0$, there are two types of QSLs given by Eq. (1)
— a gapless phase when
$|J_{1}|\leqslant|J_{2}|+|J_{3}|,|J_{2}|\leqslant|J_{1}|+|J_{3}|,|J_{3}|\leqslant|J_{2}|+|J_{1}|$,
and gapped phases otherwise. Here we consider the case $J_{1}=J_{2}=J_{3}=J$
and $h_{1}h_{2}h_{3}\neq 0$, then $K\approx h_{1}h_{2}h_{3}/J$ which makes the
system gapped except on the edges where chiral Majorana modes exist. Kitaev
The energy spectra with open zigzag and armchair edge modes are shown in
FIG.1(b) and (c) respectively.

Figure 1: (a) The honeycomb lattice for the Kitaev model. The solid/dashed
arrows denotes the nearest/next-nearest neighbours. The dashed rectangle
denotes a unit cell used for numerical calculations. (b)/(c) The energy
spectrum with open zigzag/armchair edges and $J=2,K=0.1$. The in-gap edge
states are pointed out.

Spins can induce photon scattering through super-exchange coupling Fleury  and
the Raman operator for the bond of two spins on the sites $i$ and $j$ is (up
to an undetermined constant) Fleury ; Shastry ; Nasu ; Knolle1

$\displaystyle
R_{ij}\sim(\hat{\bm{e}}_{\text{in}}\cdot\bm{d}_{ij})(\hat{\bm{e}}_{\text{out}}\cdot\bm{d}_{ij})J^{\nu}_{ij}\sigma_{i}^{\nu}\sigma_{j}^{\nu},$
(3)

where $J^{\nu}_{ij}$ is the coupling between spins and $\bm{d}_{ij}$ is the
vector connecting them. The unit vector $\hat{\bm{e}}_{\text{in/out}}$ is the
polarization direction of the incident/outgoing light. In our case, it can be
written in Majorana representation as

$\displaystyle
R_{ij}\sim\frac{iJ}{2}(\hat{\bm{e}}_{\text{in}}\cdot\bm{d}_{ij})(\hat{\bm{e}}_{\text{out}}\cdot\bm{d}_{ij})c_{i}c_{j}.$
(4)

The Raman scattering intensity for a given photon frequency shift $\omega$ is

$\displaystyle I(q,\omega)=\int_{-\infty}^{\infty}dte^{i\omega
t}\langle\hat{R}(q,t)\hat{R}(-q,0)\rangle.$ (5)

where $\hat{R}(q,t)$ is the Fourier component along the $x$-direction of
Eq.(4) in the Heisenberg picture.

We are interested in the quantity $I(q,\omega)$ induced by the one-dimensional
Majorana edge modes. Before a calculation in the full two-dimensional model,
we study tentatively a one-dimensional effective edge theory with the Majorana
Hamiltonian

$\displaystyle H_{\rm eff}=\sum_{k>0}vkc_{k}^{\dagger}c_{k}.$ (6)

And the one-dimensional version of Eq.(4) is (in the reciprocal space)
$R(q)\sim v\sum_{k}\sin kc_{-k-q}c_{k}$ which becomes
$R(q)=\sum_{k}vkc_{k+q}^{\dagger}c_{k}$ in the continuum limit. Substitution
of this operator into Eq.(5) and application of the effective Hamiltonian
$H_{\rm eff}$ lead to (using Wick’s theorem and assuming the temperature
$T=0$)

$\displaystyle\langle\hat{R}(q,t)\hat{R}(-q,0)\rangle\sim\sum_{k,k^{\prime}}kk^{\prime}e^{-ivqt}[\delta_{k,-k^{\prime}}+\delta_{k+q,k^{\prime}}],$
(7)

when $0<k<-q$ or $0<k^{\prime}<q$, and zero otherwise. $\delta_{k,k^{\prime}}$
is the Kronecker delta. Assuming infinite size and changing the wave vector
summation to an integral, we obtain
$\langle\hat{R}(q,t)\hat{R}(-q,0)\rangle\sim v^{2}|q|^{3}e^{i|vq|t}.$ Thus,
using Eq. (5),

$\displaystyle I(q,\omega)\sim v^{2}|q|^{3}\delta(\omega+|vq|).$ (8)

Assuming that the incident light has a Gaussian distributed intensity
proportional to $e^{-x^{2}/\Delta_{x}^{2}}$, the total Raman intensity is

$\displaystyle I(\omega)=\int
dqe^{-(q\Delta_{x})^{2}/2}I(q,\omega)\sim\frac{|\omega|^{3}}{v^{2}}e^{-(\omega\Delta_{x})^{2}/2v^{2}}$
(9)

for $\omega<0$ and it vanishes when $\omega>0$. The constant $\Delta_{x}$
denotes the width of the light distribution. At small $\omega$, we get
$I\sim|\omega|^{3}$. The factor $1/v^{2}$ is related to the density of states,
which is expected to affect the scattering strength. Eq. (9) has a maximum of

$\displaystyle I(\omega_{0})\sim v/\Delta_{x}^{3}$ (10)

at the frequency shift

$\displaystyle\omega_{0}=\sqrt{3}v/\Delta_{x}.$ (11)

When $\omega\ll\omega_{0}$, $I(\omega)$ is proportional to $\omega^{3}$.
$I(\omega_{0})$ decreases as $\Delta_{x}$ increases because the Raman
scattering process of the chiral Majorana modes requires translational
symmetry breaking.

Figure 2:  The zero-temperature local Raman spectra for various choices of
light-polarizations on (a) zigzag edges and (b) armchair edges. The
solid/dashed curves are obtained with the detection spot on the edges/in the
bulk. The size of the detection spot is $n_{x}\times n_{y}$, with $n_{x}$ and
$n_{y}$ being the number of rectangular unit cells along $x$\- and
$y$-directions respectively. (a)$n_{x}\times n_{y}=1\times 2$. The number of
sites along $y$-direction is $L_{y}=20$. (b), $n_{x}\times n_{y}=2\times 1$.
The number of sites along $x$-direction is $L_{x}=30$. $J=2,K=0.1$ in both (a)
and (b). The symbols XX, YY and XY denote the polarizations of the incident
(the first letter) and outgoing (the second letter) lights.

Back to the two-dimensional Kitaev model, there are two typical kinds of edges
when open boundaries are considered — zigzag and armchair. Let us assume that
the light shines on a small rectangular area of size $n_{x}\times n_{y}$
($n_{x}$ and $n_{y}$ denote the numbers of rectangular unit cells, as shown in
FIG. 1(a), along $x$\- and $y$-directions respectively.). Then the zero-
temperature Raman intensity with various polarizations and at different
positions (bulk or edge) are shown in FIG.2.

The solid (dashed) curves are obtained with the light-shining area on the edge
(in the bulk). In the low-frequency regime, the Raman intensity is finite on
the edges but vanishes in the bulk, manifesting the existence of the edge
states. A closer look at this regime shows that the intensity is proportional
to $\omega^{3}$ (see the insets in FIG.2), as indicated by previous one-
dimensional effective model analysis. Thus, a local Raman spectroscopy
measurement provides a valid method to detect the Majorana edge modes in the
QSL model. The topological gap suggested by the bulk response is around
$2E_{g}=0.5$, consistent with the energy spectrum in FIG. 1(b) and (c). For
$\hbar\omega>E_{g}$, the bulk states is involved and the $\omega^{3}$ feature
is gone. The Raman scattering in this regime is contributed by the optical
transition between the edge and bulk states. In the case of zigzag edges, the
band edge is at momenta (about $\pi\pm 0.35\pi$) away from the Majorana-edge-
mode-crossing point $\pi$, as seen in FIG. 1(b). Since the Raman scattering
with $\hbar\omega\gtrsim E_{g}$ comes from the transitions between edge modes
with $k_{x}\approx\pi$ and the bulk states on the band edge, a momentum
mismatch exists, which suppresses the Raman intensity. This is the reason that
the Raman intensity decreases just after $\hbar\omega\approx E_{g}$, resulting
in a bump-like feature of the curve in FIG. 2(a). On the other hand, there is
no momentum mismatch when the armchair edges are considered. Thus the Raman
intensity in FIG. 2(b) is monotonically increasing. The symbols XX, YY and XY
denote the polarizations of the incident (the first letter) and outgoing (the
second letter) lights. The signals are the strongest when both the incident
and outgoing lights are polarized along the edge, i.e. XX for zigzag edges and
YY for armchair edges.

An exact calculation of finite temperature results needs to consider the
complete Hilbert space including higher-energy sectors with vortex
excitations. However, we limit ourselves to finite but small (compared to the
vortex excitation gap) temperatures. Then the excitation of vortices is
exponentially suppressed and the effect of temperature is reflected on the
Fermi distribution functions of the low-energy Majorana states. In this way,
we obtain the low-temperature behaviour as shown in FIG. 3. For $\omega=0$,
the intensity is proportional to $T^{3}$.

Figure 3: Temperature dependences of the local Raman intensities (XX
polarizations) for various frequencies on the zigzag edge (solid curves) and
in the bulk (dashed). The inset is the intensity on the edge at $\omega=0$
drew with the horizontal axis changed to $(10k_{B}T)^{3}$.

## Topological superconductors —

One may expect similar Raman scattering features to appear with chiral
Majorana edge modes in other systems such as topological superconductors
(TSCs). To compare the Raman scattering by Majorana edge modes in QSLs and
TSCs, let us consider a spinless $p$-wave superconductor whose Hamiltonian is

$\displaystyle
H_{p}=\sum_{\bm{k}}{\xi(\bm{k})}\psi_{\bm{k}}^{\dagger}\psi_{\bm{k}}+\Delta[(\sin
k_{x}+i\sin k_{y})\psi_{\bm{k}}^{\dagger}\psi_{-\bm{k}}^{\dagger}+h.c.]$ (12)

where the operator $\psi_{\bm{k}}$ annihilates an electron of energy
$\xi(\bm{k})=-2t[\cos(k_{x})+\cos(k_{y})]-\mu$.

The Raman operator in models of itinerant electrons is essentially given by
the density operator. In tight-binding models, it becomes Abrikosov  (assuming
$\hat{\bm{e}}_{\text{in}}=\hat{\bm{e}}_{\text{out}}=\hat{x}$ without loss of
generality),

$\displaystyle
R_{p}(\bm{q})\sim\sum_{\bm{k}}\psi^{\dagger}_{{\bm{k}}+\bm{q}}\frac{\partial^{2}\xi(\bm{k})}{\partial
k_{x}^{2}}\psi_{\bm{k}},$ (13)

By rewriting Eqs. (12-13) in the real-space representation along $y$-direction
and applying open boundary conditions there, we obtain the Raman intensity for
this system using Eq. (5) as shown in FIG. 4. Both the spectrum in the inset
and the bulk response indicate that the bulk energy gap is around
$2E_{g}=0.54$. The results for the two edges coincide due to the obvious
spatial inversion symmetry. Most importantly, no Raman scattering happens for
frequencies $\omega<E_{g}$. This is a remarkable difference from the results
for the QSL model.

The vanishing of the edge Raman scattering in TSCs may be understood by noting
that the density operator,
$\rho(x)=\psi(x)^{\dagger}\psi(x)=\frac{1}{2}[1+ic(x)\tilde{c}(x)]$ (the
normal fermion operator is expressed in terms of Majorana operators by
$\psi(x)=c(x)+i\tilde{c}(x)$), which also gives the Raman operator here, must
involve two Majorana fields $c$ and $\tilde{c}$. Since only one Majorana field
exists below energy $E_{g}$, Raman scattering cannot happen. When
$\omega>E_{g}$, the bulk states get involved, with which another species of
Majorana modes can be defined although they are not energy eigenstates. Thus,
the Raman intensity becomes finite. The Raman response in the range
$E_{g}<\omega<2E_{g}$ is a result of the collaboration of edge and bulk
states, similar to the QSL case.

Figure 4: The Raman intensities in the bulk and on the edge of a spin-less
$p$-wave superconductor. The inset shows the energy spectrum. The size along
the $y$-direction is $80$ sites. The parameters $t=1,\mu=1,\Delta=0.3$.

Comparing the Raman scattering in the QSL and that in the TSC makes it clear
that the electrodynamics of Majorana edge modes depends on their environments,
i.e. on the bulk properties. This is related to the fact that Majorana modes
do not preserve charges themselves. Charge conservation is recovered only when
the whole system, including the bulk and the edges, is taken into account.
Thus, is no general theories of the electrodynamics of Majorana edge modes.

Although Majorana edge states alone in TSCs do not contribute to Raman
scattering, they can absorb light and induce an optical conductivity (the real
part) $\Re[\sigma_{xx}]$ that is proportional to $\omega^{2}$ JJH2021 . In a
Kitaev spin liquid, photon absorption involves the gapped vortex excitations
and thus is negligible at the low-energy limit. Bolens  We obtain TABLE 1 as a
comparison between QSLs and TSCs.

| $\Re[\sigma_{xx}]$ | Raman
---|---|---
QSL | 0 | $\sim\omega^{3}$
TSC | $\sim\omega^{2}$ | 0
Table 1: Raman scattering and real-part optical conductivity due to Majorana
edge modes in QSLs and TSCs.

## Conclusion and discussion—

We have shown that the chiral Majorana edge modes in the Kitaev model of QSL
induce Raman scattering whose intensity is proportional to the frequency cube
$\omega^{3}$ in the small-$\omega$ limit. For a finite temperature $T$ that is
small compared to the bulk gap, the Raman intensity is proportional to $T^{3}$
when $\omega\rightarrow 0$. The bulk system must be taken into account when
electrodynamics of Majorana edge modes is considered. Majorana edge modes from
different origins may lead to different electrodynamic properties. For
instance, the Majorana edge modes alone in TSCs cannot induce Raman scattering
while they can absorb light.

A candidate material that shows the properties of the Kitaev model is
$\alpha$-RuCl3 Sandilands ; Jansa . The spin-excitation gap of this material
is about $\Delta_{0}=1$meV Jansa . The topological gap $\Delta_{1}$ opened by
an in-plane magnetic field $B$ is about $\Delta_{1}\approx 4\times
10^{-3}B^{3}$meV$\cdot$ T-3 Jansa . If $B=5$T, then $\Delta_{1}=0.5$meV and
the Raman scattering purely due to the chiral Majorana edge modes should be
detected with the photon energy shift $\hbar\omega<0.25$meV (or the frequency
shift $\omega<0.38$THz) and the temperature $T\ll 10$K.

Now we estimate the Raman intensity based on our numerical results. Figure
1(b) indicates $v\sim 1/\Delta_{1}$. Thus, a stronger field and a smaller size
results in larger Raman intensity according to Eq. (10). Since $J\approx
30$meV Jansa , the gap in our QSL model calculation as shown in FIG. 2 is
about $0.25J\approx 8$meV. The real gap $\Delta_{1}=0.5meV$ is $16$ times
smaller and thus the real edge state velocity should be
$v_{real}=v_{model}/16$. The peak position of FIG. 2(a) is about half of the
gap, i.e. $\hbar\omega_{0,model}\approx 4$meV. The $\Delta_{x}$ we used in our
calculation is $\Delta_{x,model}\approx 1$nm. According to Eq. (11), a
realistic peak position with $\Delta_{x,real}=5nm$ is
$\hbar\omega_{0,real}\approx\frac{1}{16\times 5}\hbar\omega_{0,model}=0.05$meV
or $\omega_{0,real}=75$GHz. Since $I(\omega_{0,model})\approx 4\times
10^{-3}I_{bulk}$ according to FIG. 2(a), we obtain
$I(\omega_{0,real})=\frac{1}{16\times 5^{3}}I(\omega_{0,model})\approx 2\times
10^{-6}I_{bulk}$. As a reference, the bulk Raman intensity $I_{bulk}$ in our
model calculation is of the order unity, which corresponds to the value
measured by Sandilands et al.Sandilands . It may seem too small a signal to be
observed if a factor of $10^{-6}$ is there. However, the gigantic local field
generated by a scanning tunneling microscope (STM) tip can enhance the Raman
scattering by a factor of $10^{7}$ or even larger Zhang2017 ; Verma2017  and
thus the signal can be well detectable.

Real materials may be a stack of many layers instead of rather than exactly
two-dimensional. In principle, the edge states of each layer may couple with
each other and form two-dimensional surface states, instead of the one-
dimensional modes we have discussed here. However, a Heisenberg coupling
between layers involves vortex excitations which are gapped by $\Delta_{0}$.
So the inter-layer coupling is irrelevant in the low-energy limit.

###### Acknowledgements.

## Acknowledgments —

N.N. was supported by Ministry of Education, Culture, Sports, Science, and
Technology Nos. JP24224009 and JP26103006, the Impulsing Paradigm Change
through Disruptive Technologies Program of Council for Science, Technology and
Innovation (Cabinet Office, Government of Japan), and Core Research for
Evolutionary Science and Technology (CREST) No. JPMJCR16F1 and No. JPMJCR1874,
Japan. JJH was supported by the RIKEN Incentive Research Projects.

## References

  * (1) N. Read and D. Green, Phys. Rev. B 61, 10 267 (2000).
  * (2) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * (3) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * (4) F. Wilczek, Nat. Phys. 5, 614 (2009).
  * (5) K. T. Law, P. A. Lee and T. K. Ng, Phys. Rev. Lett. 103, 237001 (2009).
  * (6) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * (7) S. Fujimoto, Phys. Rev. B 77, 220501(R) (2008).
  * (8) K. Nomura, S. Ryu, A. Furusaki, and N. Nagaosa, Phys. Rev. Lett. 108, 026802 (2012).
  * (9) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
  * (10) J. Alicea, Y. Oreg, G. Refael, F. von Oppen and M. P. A. Fisher Nat. Phys. 7,412 (2011).
  * (11) B. Lian, X. -Q. Sun, A. Vaezi, X. -L. Qi and S. -C. Zhang, Proc. Nat. Acad. Sci. USA 115, 10938 (2018).
  * (12) A. Kitaev, Annals of Physics 321, 2–111 (2006).
  * (13) Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Nature 559, 227 (2018).
  * (14) L. Fu and C. L. Kane, Phys. Rev. Lett. 102, 216403 (2009).
  * (15) A. R. Akhmerov, J. Nilsson, and C. W. J. Beenakker, Phys. Rev. Lett. 102, 216404 (2009).
  * (16) X. L. Qi, T. L. Hughes and S. C. Zhang, Phys. Rev. B 82, 184516 (2010).
  * (17) T. Meng and L. Balents, Phys. Rev. B 86, 054504 (2012).
  * (18) W. Ji and X. G. Wen, Phys. Rev. Lett. 120, 107002 (2018).
  * (19) J. J. He, T. Liang, Y. Tanaka and N. Nagaosa, Commun. Phys. 2, 149 (2019).
  * (20) Q. L. He, L. Pan, A. L. Stern, E. C. Burks, X. Che, G. Yin, J. Wang, B. Lian, Q. Zhou, E. S. Choi, K. Murata, X. Kou, Z. Chen, T. Nie, Q. Shao, Y. Fan, S. C. Zhang, K. Liu, J. Xia, and K. L. Wang, Science 357, 294 (2017).
  * (21) Z. Wang, J. O. Rodriguez, L. Jiao, S. Howard, M. Graham, G. D. Gu, T. L. Hughes, D. K. Morr, and V. Madhavan, Science 367, 104 (2020).
  * (22) J. J. He, Y. Tanaka, and N. Nagaosa, to be published.
  * (23) J. Knolle, G.-W. Chern, D. L. Kovrizhin, R. Moessner, and N. B. Perkins, Phys. Rev. Lett. 113, 187201 (2014).
  * (24) J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, and R. Moessner, Nature Phys 12, 912–915 (2016).
  * (25) L. J. Sandilands, Y. Tian, K. W. Plumb, Y.-J. Kim, and K. S. Burch, Phys. Rev. Lett. 114, 147201 (2015).
  * (26) W. Zhang, Z. Fang, and X. Zhu, Chem. Rev. 117, 5095 (2017).
  * (27) P. Verma, Chem. Rev. 117, 6447 (2017).
  * (28) A. Jarzembski, C. Shaskey, and K. Park, Front. Energy 12, 43 (2018).
  * (29) Jordan M. Klingsporn, Nan Jiang, Eric A. Pozzi, Matthew D. Sonntag, Dhabih Chulhai, Tamar Seideman, Lasse Jensen, Mark C. Hersam, and Richard P. Van Duyne, J. Am. Chem. Soc. 136, 3881 (2014).
  * (30) K. Lee, M. I. B. Utama, S. Kahn, A. Samudrala, N. Leconte, B. Yang, S. Wang, K. Watanabe, T. Taniguchi, M. V. P. Altoé, G. Zhang, A. Weber-Bargioni, M. Crommie, P. D. Ashby, J. Jung, F. Wang, and A. Zettl, Science Advances 6, eabd1919 (2020).
  * (31) K. Lai, W. Kundhikanjana, M. A. Kelly, Z.-X. Shen, J. Shabani, and M. Shayegan, Phys. Rev. Lett. 107, 176809 (2011).
  * (32) M. Allen, Y.-T. Cui, E. Y. Ma, M. Mogi, M. Kawamura, I. C. Fulga, D. Goldhaber-Gordon, Y. Tokura, and Z.-X. Shen, Proceedings of the National Academy of Sciences, 116(29), 14511–14515 (2019).
  * (33) A. Bolens, H. Katsura, M. Ogata, and S. Miyashita, Phys. Rev. B 97, 161108 (2018).
  * (34) X. Feng, G. Zhang and T. Xiang, Phys. Rev. Lett. 98, 087204 (2007).
  * (35) P. A. Fleury and R. Loudon, Phys. Rev. 166, 514 (1968).
  * (36) B. S. Shastry and B. I. Shraiman, Phys. Rev. Lett. 65, 1068 (1990).
  * (37) A. Abrikosov and V. Genkin, Sov. J. Exp. Theor. Phys. 38, 842 (1974).
  * (38) N. Jansa, A. Zorko, M. Gomilšek, M. Pregelj, K. W. Krämer, D. Biner, A. Biffin, C. Rüegg, and M. Klanjšek, Nature Physics 14, 786-790 (2018).


# Integrative analysis of ATAC-seq and RNA-seq for cells
infected by human T-cell leukemia virus type 1

Azusa Tanaka  Department of Human Genetics, Graduate School of Medicine, The
University of Tokyo  Yasuhiro Ishitsuka  Institute of Mathematics for
Industry, Kyushu University  Hiroki Ohta  Department of Human Sciences,
Obihiro University of Agriculture and Veterinary Medicine
Norihiro Takenouchi  Department of Microbiology, Kansai Medical University
Masanori Nakagawa  Kyoto Prefectural University of Medicine  Ki-Ryang Koh
Department of Hematology, Osaka General Hospital of West Japan Railway Company
Chiho Onishi  Laboratory of Ultrastructural Virology, Institute for Life and
Medical Sciences, Kyoto University
Hiromitsu Tanaka  Department of Biophysics, Graduate school of Science, Kyoto
University  Department of Developmental Biology, Graduate School of Medicine,
Chiba University  Akihiro Fujimoto  Department of Human Genetics, Graduate
School of Medicine, The University of Tokyo  Jun-ichirou Yasunaga  Department
of Hematology, Rheumatology and Infectious Disease, Faculty of Life Sciences,
Kumamoto University  Masao Matsuoka  Department of Hematology, Rheumatology
and Infectious Disease, Faculty of Life Sciences, Kumamoto University

###### Abstract

Human T-cell leukemia virus type 1 (HTLV-1) causes adult T-cell leukemia (ATL)
and HTLV-1-associated myelopathy (HAM) after a long latent period in a
fraction of infected individuals. These HTLV-1-infected cells typically have
phenotypes similar to that of CD4+ T cells, but the cell status is not well
understood. To extract the inherent information of HTLV-1-infected CD4+ cells,
we integratively analyzed the ATAC-seq and RNA-seq data of infected cells.
Compared to CD4+ T cells from healthy donors, we found anomalous chromatin
accessibility in HTLV-1-infected CD4+ cells derived from ATL cases in terms of
location and sample-to-sample fluctuations in open chromatin regions. Further,
by focusing on systematically selected genes near the open chromatin regions,
all the gene expressions in ATL cases were found to be distinct from those of
healthy CD4+ T cells. Based on a further analysis of chromatin accessibility,
we detected TLL1 (Tolloid Like 1) as one of the key genes that exhibit unique
gene expressions in ATL cases. A luciferase assay indicated that TLL1 has a
strong regulatory effect on TGF-$\beta$. Overall, this study provides results
about the status of HTLV-1 infected cells, which are qualitatively consistent
across the different scales of chromatin accessibility, transcription, and
immunophenotype.

##  I Introduction

It has been statistically estimated that there are more than $300,000$ types
of mammalian host viruses doi:10.1128/mBio.00598-13 . Among the many viruses
that have been discovered, only a few have been reported to cause cancers,
such as the DNA virus human papillomavirus (HPV) and the RNA virus hepatitis C
virus (HCV) KrumpYou2018 . One, human T-cell leukemia virus (HTLV-1), is a
retrovirus, oncogenic virus, and estimated to infect approximately $10$
million people worldwide (Gessain2012, ).

Adult T-cell leukemia (ATL) and HTLV-1-associated myelopathy (HAM) are both
associated with prior infection with HTLV-1. However, these two diseases have
different clinical and pathological presentations (Matsuoka2007, ;
Bangham2023, ). The genes encoded by HTLV-1, such as HBZ (HTLV-1 basic leucine
zipper factor) and Tax, have been reported to affect important signaling
pathways involved in cell proliferation, apoptosis, and infectivity
(Yamagishi2018, ). In particular, HBZ is maintained in all ATL cases and
functions as both a protein and RNA (Satou2006, ; Tanaka-Nakanishi2014, ;
Mitobe2015, ; Ma2016, ). Recent studies have elucidated that in ATL cells,
genomic mutations are highly enriched in T cell-related pathways, such as
NF-${\kappa}$B, and typically activate the pathways (Kataoka2015, ;
Kataoka2016_2, ). Furthermore, it has been frequently observed in ATL cases
that the aberrant expression of programmed cell death 1-ligand 1 (PD-L1) is
caused by disruption of the PD-L1 3’-untranslated region (UTR) (Kataoka2016_1,
).

Several questions about these diseases at the genomic scale remain, including
how the chromatin structure of ATL cells differs from that of CD4+ T cells
derived from healthy donors and how this difference influences transcription
and translation to finally cause symptoms. In general, cellular phenotypes are
largely affected by gene expressions that are strongly correlated with the
epigenetic mechanisms occurring in chromatin. To understand the epigenetic
mechanisms, it is important to understand how human DNA is packed and
chemically modified in the nucleus, which can be quantified by measuring
chromatin accessibility.

In this paper, we study the relationship between chromatin accessibility and
transcription in HTLV-1-infected cells at the whole genome level using Assay
for Transposase-Accessible Chromatin using sequencing (ATAC-seq) (ATAC1, ) and
RNA sequencing (RNA-seq) data. We performed a comparative analysis of HTLV-1
infected CD4+ cells from ATL cases, HAM cases, and CD4+ T cells from healthy
donors (healthy CD4+ T cells) based mainly on our previously developed
algorithm of systematic clustering (Tanaka2020, ).

Our analysis shows that CD4+ cells derived from ATL cases have anomalous
properties in terms of the locations and sample-to-sample fluctuations of open
chromatin regions compared with healthy CD4+ T cells. Additionally, genes
selected by our systematic clustering algorithm based on the immunophenotype
had distinct expressions between ATL cells and healthy CD4+ T cells. Using our
systematic clustering algorithm, we also found a relationship between
chromatin accessibility and immunophenotype that suggests some ATL cases
approach several types of myeloid cells.

Finally, we detected TLL1 (Tolloid Like 1) as one of a few genes having
anomalous expressions in ATL cases. A luciferase assay found TLL1 isoforms,
depending on the types of its isoforms, regulate differently the maturation of
TGF-$\beta$ (transforming growth factor $\beta$), which is known to play
important roles in cancer progression.

##  II Results

Figure 1: The genetic (intergenic/intronic/exonic) annotation of ATAC-seq
peaks, which quantify open chromatin regions. TSS(2kb $+/-$) corresponds to
$-2000$ to $2000$ base pairs from a transcription start site.

###  II.1 Chromatin accessibility: whole view of the genome

(a)

(b)

State No. | Description | $N_{\mathrm{CD4^{+}T}}$ | $N_{\mathrm{ATL}}$ | $N_{\mathrm{HAM}}$ |  $\dfrac{N_{\mathrm{ATL}}}{N_{\mathrm{CD4^{+}T}}}$ | $\dfrac{N_{\mathrm{HAM}}}{N_{\mathrm{CD4^{+}T}}}$
---|---|---|---|---|---|---
1 | Active TSS | 9939 | 12084 | 12691 | 1.22 | 1.28
2 | Flanking Active TSS | 12202 | 7015 | 6983 | 0.575 | 0.572
3 | Transcr. at gene 5’ and 3’ | 489 | 213 | 224 | 0.436 | 0.458
4 | Strong transcription | 5525 | 3278 | 893 | 0.593 | 0.162
5 | Weak transcription | 10659 | 18946 | 7269 | 1.78 | 0.682
6 | Genic enhancers | 4431 | 1374 | 1017 | 0.310 | 0.230
7 | Enhancers | 32260 | 26597 | 24283 | 0.824 | 0.753
8 | ZNF genes & repeats | 401 | 595 | 559 | 1.48 | 1.39
9 | Heterochromatin | 1340 | 794 | 565 | 0.593 | 0.422
10 | Bivalent/Poised TSS | 478 | 181 | 192 | 0.379 | 0.402
11 | Flanking Bivalent TSS/Enh | 1894 | 326 | 268 | 0.172 | 0.141
12 | Bivalent Enhancer | 2031 | 1074 | 682 | 0.529 | 0.336
13 | Repressed PolyComb | 1746 | 1077 | 382 | 0.617 | 0.219
14 | Weak Repressed PolyComb | 9086 | 12178 | 4445 | 1.34 | 0.489
15 | Quiescent/Low | 38566 | 92307 | 28746 | 2.39 | 0.745

Figure 2: (a) The number of ATAC-seq peaks (vertical axis) vs. indices
(horizontal axis) classified by each functional annotation from $1$ to $15$,
as shown in (b). (b) The number of ATAC-seq peaks classified into each
functional annotation for healthy CD4+T cells, ATL cells, and HAM cells.
$N_{\mathrm{CD4^{+}T}}$, $N_{\mathrm{ATL}}$, and $N_{\mathrm{HAM}}$: peak
number of healthy CD4+T, ATL, and HAM, respectively. Note that each peak
quantifies an open chromatin region.

The landscape of chromatin accessibility provides useful information for
understanding the mechanisms that govern cell-type-specific gene expressions.
Preliminarily, we overview the chromatin accessibility characterized by the
ATAC-seq of healthy CD4+ T cells, ATL cells, and HAM cells.

To obtain the chromatin accessibility landscape, we performed ATAC-seq on
HTLV-1-infected CD4+ cells obtained from the peripheral blood of $29$ ATL and
$6$ HAM cases. All samples selected for the ATAC-seq library preparation were
at least $98$% HTLV-1-infected cells. The ATAC-seq libraries were sequenced
with an average of $44$ million reads, resulting in a dataset comprising of
$1.3$ billion and $556$ million sequenced reads for ATL and HAM, respectively.
The data quality was high in all cases, with mitochondrial read rates of
$7.5$% for ATL and $7.3$% for HAM. For a comparison, we used ATAC-seq datasets
of CD4+ T cell samples from $5$ healthy donors and downloaded from GEO
accession GSE74912 (Howard, ).

To identify genome-wide accessible chromatin regions, for each of the three
groups, we concatenated ATAC-seq reads for the different samples, where the
sample number was $29$ for ATL, $6$ for HAM, and $5$ for healthy CD4+ T cells.
As explained in the Materials and Methods, we randomly selected $100$M reads
from the concatenated data of each group. We used the MACS2 algorithm to
select the locations of peaks to quantify the open chromatin regions from the
ATAC-seq datasets (Zhang2008, ), finding a total of $178811$, $89972$, and
$131609$ peaks in ATL, HAM, and healthy CD4+ T cells, respectively.

The ENCODE consortium shows that $10$% of peaks are localized in near
transcription start sites (TSSs), whereas the remaining $90$% of peaks are
mapped nearly equally to intronic and intergenic regions (Thurman2012, ).
Consistent with these data, as shown in Fig. 1, about $10$% of the ATAC-seq
peaks are overlapped with the TSSs and their surrounding regions, whereas the
majority of ATAC-seq peaks (about $85$%) of healthy CD4+ T cells, ATL cells,
and HAM cells reside in intergenic or intronic regions.

To determine the functional roles of the peaks in HTLV-1-infected cells and
healthy CD4+ T cells, we computed the overlapping ratio of these regions with
specific genomic features, such as active TSS, enhancers, heterochromatin,
etc. To assign a genomic feature to all genomic positions, we assumed a
chromHMM $15$-state model obtained from
(https://egg2.wustl.edu/roadmap/web_portal/chr_state_learning.html). We used
the data of E043 for healthy CD4+ T cells and E037 for HTLV-1-infected CD4+
cells. Note that E037 is the model of CD4+ memory T cells because a majority
of ATL cases has been reported to show CD45RO+, which is consistent with CD4+
memory T cells (Richardson1990, ).

As shown in Fig. 2, for HTLV-1-infected CD4+ cells, the number of peaks
compared with healthy CD4+ T cells was proportionally less in categories
related to Enhancer $(6,7,11,12)$ and Heterochromatin $(9)$, but it was higher
for the category of Quiescent/Low $(15)$ only for ATL cases.

This observation suggests that compared with healthy CD4+ T cells, distinct
enhancer mechanisms in HTLV1-infected cases are correlated with the distinct
chromatin structures.

(a)

(b)

Figure 3:  (a) The normalized number $\widetilde{\rho}_{\nu}(z)$ of reads at
position $z$ from the TSS for type $\nu$ as defined in (3). Note that position
$z$ takes the value $\\{-2000+20k\mid 0\leq k\leq 200\\}$, where $k$ is an
integer. (b) A histogram (heat map) $F^{\Delta,\xi}_{\nu}(z,\ell)$ of the
fragment length $\ell$ of the ATAC-seq data at position $z$ from the TSSs as
defined in (5). (left) CD4+T, (center) ATL, and (right) HAM. The graphs are
plotted for $z\in\\{n\Delta\mid-37\leq n\leq 37\\}$ and $\ell\in\\{n\xi\mid
0\leq n\leq 49\\}$, where $\Delta=1000/75$ and $\xi=10$.

###  II.2 Increased chromatin accessibility around transcriptional start sites
(TSSs) in ATL

(a) All regions

(b) Coding regions

(c) Non-coding regions

Figure 4: Histograms $F^{(1)}_{\mathbb{L},\mathbb{S}}(O;\Delta)$ of the
averaged length $O$ of peaks overlapping with the reference open region as
defined in (13) with bin width $\Delta$. $\mathbb{S}$ consists of $6$ HAM
samples, $29$ ATL samples, and $5$ healthy CD4+T samples. (a) All regions with
$\mathbb{L}=\mathbb{G}\cup\mathbb{G}^{c}$. (b) Coding regions with
$\mathbb{L}=\mathbb{G}$. (c) Non-coding regions with
$\mathbb{L}=\mathbb{G}^{c}$.

To determine how the chromatin structures observed in the HTLV1-infected cases
are statistically characterized depending on the positions in the genome, we
examined the positions of the reads from the ATAC-seq data.

We plotted a histogram $\widetilde{\rho}_{\nu}(z)$ of the reads as a function
of their positions $z$ relative to TSSs for cell type $\nu$, where $\nu$ is an
ATL, HAM, or healthy CD4+ T cell type. Note that the positions of the TSSs and
coding regions of all genes were obtained from the human genome (hg19). For
technical details of the histogram, see Materials and Methods. As shown in
Fig. 3a, the tail parts of the histogram for ATL cases take higher values
compared with HAM and healthy CD$4^{+}$ T cells, whereas the latter two cases
showed more similar forms as a whole.

To elucidate the statistics of the chromatin structures around the TSSs, we
also focused on the fragments, which can be reconstructed by the reads data;
both ends of a fragment correspond to a pair of two reads. Specifically, we
investigated the position-dependent accumulation of the fragments, which can
be an estimate of nucleosome positioning. We plotted a heat map
$F^{\Delta,\xi}_{\nu}(z,\ell)$ in which the mid-point of each fragment
relative to TSSs is placed on the horizontal axis as $z$ and the length of
each fragment is placed on the vertical axis as $\ell$ (Henikoff2011, ). For
technical details of the heat map, see Materials and Methods.

As shown in Fig. 3b, healthy CD4+ T cells and HTLV-1-infected CD4+ cells from
HAM samples show a pattern of enriched nucleosome-free fragments ($\ell<100$
bp) and mono-nucleosome fragments ($\ell=180\sim 247$ bp) surrounding the
TSSs, where the thresholds for the length of a fragment are $100$ bp, $180$
bp, and $247$ bp based on a previous study (Mazumdar2015, ). ATL samples
showed less enrichment of nucleosome-free and mono-nucleosome fragments.

These observations suggest that the statistics of open chromatin regions and
nucleosome positioning around the TSSs in ATL cases is distinct from HAM cases
and healthy CD4+ cells, both of which showed more similar forms again. We
continue to elaborate on characteristic behaviors of ATL cells distinct from
the other two cases.

###  II.3 Giant sample-to-sample fluctuations of chromatin accessibility in
ATL

(a) All regions

(b) Coding regions

(c) Non-coding regions

(d) All regions

Figure 5: Histograms $F^{(2)}_{\mathbb{L},\mathbb{S}}(O;\Delta)$ of the
variance in length $O$ of peaks overlapping the reference open region as
defined in (14) with bin width $\Delta$. $\mathbb{S}$ consists of $6$ HAM
samples, $29$ ATL samples, and $5$ CD4+T samples for (a), (b), and (c). (a)
All regions with $\mathbb{L}=\mathbb{G}\cup\mathbb{G}^{c}$. (b) Coding regions
with $\mathbb{L}=\mathbb{G}$. (c) Non-coding regions with
$\mathbb{L}=\mathbb{G}^{c}$. (d) All regions with
$\mathbb{L}=\mathbb{G}\cup\mathbb{G}^{c}$ for CD4+ T only, HAM only, and a
mixed set of CD4+ T and HAM with $5$ samples for each cell type.

To examine the distribution of open chromatin regions, we applied a systematic
clustering algorithm that robustly detects open chromatin regions relevant to
classifying the immunophenotype (Tanaka2020, ).

This algorithm characterizes open/closed chromatin regions in the following
way. First, MACS2 uses ATAC-seq reads data for a given sample $s$ as the input
and outputs a collection of peak data including the locations of peaks with
their $p$-values. These peaks are considered as open regions of chromatin for
sample $s$. The peaks are ordered in ascending $p$-values, and the first $M$
peaks are taken. The set of $M$ peaks are written as
$\hat{g}_{s}^{M}=((\gamma_{k},\alpha_{k},\beta_{k}),p_{k})_{k=1}^{M}$, where
the location of the $k$-th peak is the region $(\alpha_{k},\beta_{k})$ in the
$\gamma_{k}$-th chromosome, and $p_{k}$ is the $p$-value of the $k$-th peak.
Then, the set of the top $M=64000$ peaks is determined as the optimal set of
open chromatin regions, which effectively classifies the immunophenotypes of
the samples (Tanaka2020, ). In this procedure, peaks with high $p$-values
(unreliable peaks) can be treated as noise and ignored for later analysis.

First, we tried to capture the genomic positions where chromatins tend to be
open for at least one of the ATL, HAM, and healthy CD4+T types. To quantify
such chromatin regions, we constructed a new reference set $g_{0}$ of peaks as
follows: we concatenated all the reads from all the samples with the same cell
type. Then, we used MACS2 to obtain a set of peaks for a cell type of the
concatenated reads. Finally, we merged all the peaks obtained from the three
cell types. For the explicit construction of the reference set $g_{0}$, see
Materials and Methods. Next, we classified the reference set $g_{0}$ into the
set of all peaks overlapping gene-coding regions, which we denoted as
$\mathbb{G}$, and the set of all peaks overlapping with non-coding regions,
which we denoted as $\mathbb{G}^{c}$.

To quantify open chromatin regions characterized by the peaks
$\hat{g}_{s}^{M}$ of sample $s$ in each $k$-th peak from the reference set
$g_{0}$, we computed the length of overlapped peaks
$O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})$ between $g_{0}$ and
$\hat{g}_{s}^{M}$ in the peak location $(\alpha_{k},\beta_{k})$ picked up from
the reference region
$\mathbb{L}\in\\{\mathbb{G},\mathbb{G}^{c},\mathbb{G}\cup\mathbb{G}^{c}\\}$.
We set $M=64000$ as the provisionally optimal number for the immunophenotype
classification (Tanaka2020, ). Then, we focused on the average and variance of
$O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})$. For details of the calculations,
see Materials and Methods.

As shown in Fig. 4, healthy CD4+ T, HAM, and ATL cells showed similar
behaviors at large average lengths of the overlapped peaks. For small average
lengths, CD4+ T showed fewer peaks, compared with ATL cases. Additionally, as
shown in Fig. 5a, we found that healthy CD4+ T and HAM cases had similar
sample-to-sample fluctuations, and ATL cases had a higher frequency at
variances larger than $10^{5}$. As shown in Fig. 5b and Fig. 5c, the larger
sample-to-sample fluctuations in the ATL cases were found in both non-coding
regions and coding regions. In contrast, Fig. 5b shows apparent gaps between
ATL cases and the other two cases at intermediate variances only in coding
regions; ATL cases had a higher frequency at intermediate variances around
$10^{3}$ compared with the other two cases.

The above analysis does not clearly distinguish healthy CD4+ T and HAM cells,
in particular, with respect to the sample-to-sample fluctuations of the open
chromatin regions. Therefore, we mixed the datasets of 5 cases of HAM cells
and 5 cases of healthy CD4+ T cells to give set $\mathbb{S}$. As shown in Fig.
5d, the variance in the mixed dataset was larger than of HAM cells or healthy
CD4+ T cells alone. This finding indicates the distributions of open chromatin
regions are different between healthy CD4+ T and HAM cells, although the two
distributions showed similar sample-dependence. As a comparision, Fig. 5d and
Fig. 5a show that the sample-to-sample fluctuations of ATL cells are larger
than the fluctuations of the mixed data except for the tail, where the samples
are scarce.

Thus, ATL cases have a higher frequency at larger sample-to-sample
fluctuations at the whole genome level and a higher frequency at intermediate
sample-to-sample fluctuations only in coding regions. On the other hand, the
chromatin structures in HAM cases show less sample-dependence, which is
similar to CD4+ T cells, implying the existence of a certain trend common to
all the HAM samples.

(a)

(c)

(b)

Figure 6:  RNA-seq statistics: (a) A PCA of RNA-seq for healthy CD4+T cells
and ATL cells, where the percentages are the first and second contribution
ratios. (b) A GO analysis of RNA-seq for detecting the top 5 up-regulated gene
expressions (top) and the top 5 down-regulated expressions (bottom) in terms
of cell function for ATL. (c) A histogram of the log-fold change calculated in
(16) and (17) between the gene expressions (RNA-seq) of ATL cases and of
healthy CD4+ T cells for all genes (blue line) vs. 61 histone-related genes
(orange line). To select relevant genes, first we chose 69 genes whose names
start with HIST. Then, we removed 8 genes, including 7 genes for which healthy
CD4+ T cells had no peak and 1 gene for which ATL cells had no peak.

###  II.4 mRNA: distinctively expressed histone modifications in ATL

To analyze gene expressions, we examined the RNA-seq data of ATL cases. Unless
noted otherwise, the samples used were 8, 10, 21, 24 in Table 1, which had
common properties in terms of immunophenotypes and symptoms, as explained
below.

We analyzed the gene expression pattern of HTLV-1-infected CD4+ cells obtained
from the peripheral blood of 4 ATL cases and of 4 healthy CD4+ T cells. The
reads count of the RNA-seq data was normalized by TMM normalization (cite-key,
) and used as the input data. For technical details, see Materials and
Methods.

As shown in Fig. 6a, a principal component analysis (PCA) shows that the gene
expression patterns differ significantly between ATL and healthy CD4+ T cells.
In ATL cells, there were $1289$ genes up-regulated based on the condition of
$\log_{2}\rm{FC}>3$ and $p$-value $<0.01$ and $944$ genes down-regulated based
on the condition of $\log_{2}\rm{FC}<-3$ and $p$-value $<0.01$, where
$\mathrm{FC}$ is the fold change of the gene expression of ATL cells relative
to healthy CD4+ T cells.

In addition, as shown in Fig. 6b, a Gene Ontology (GO) analysis using enrichR
(enrichr, ) revealed that in ATL cases, the up-regulated genes are enriched in
histone modifications. Note that the combined scores of the genes with down-
regulated expression are lower than those of up-regulated genes enriched for
histone modification. Further, as shown in Fig. 6c, many histone-related
genes, such as HIST1H2AH, HIST1H3C, and HIST1H4C, are significantly up-
regulated in ATL cases, where all genes beginning with HIST in the first 4
letters are regarded as histone-related genes. For details of the analysis,
see equation (1).

This observation led us to consider a correlation between the anomalous
chromatin properties of ATL cells shown above and the gene expression levels
found here.

(a) CD4+T vs Mono

(b) CD4+T vs ATL

(c) Mono vs ATL

(d) CD4+T vs CD8+T

(e) HSC vs ATL

(f) HSC vs CD4+T

Figure 7:  Histograms $F_{t,t_{0}}^{M}(P;\Delta)$ of the expression of the
selected genes defined in (18), where bin width $\Delta$ is equal to $0.3$.
The histograms characterize the correlation between gene expressions
quantified by RNA-seq and ATAC-seq peaks selected systematically. (a)
$M=\infty$ (orange line), $\ 64000$ (blue line) for cell types with
$(t,t_{0})=(\rm{CD4}^{+}\rm{T},\rm{Mono})$. (b) $M=\infty,\ 64000$ for cell
types with $(t,t_{0})=(\rm{CD4}^{+}\rm{T},\rm{ATL})$. (c) $M=\infty,\ 64000$
for cell types with $(t,t_{0})=(\rm{Mono},\rm{ATL})$. (d) $M=\infty,\ 64000$
for cell types with $(t,t_{0})=(\rm{CD4}^{+}\rm{T},\rm{CD8}^{+}\rm{T})$. (e)
$M=\infty,\ 64000$ for cell types with $(t,t_{0})=(\rm{HSC},\rm{ATL})$. (f)
$M=\infty,\ 64000$ for cell types with
$(t,t_{0})=(\rm{HSC},\rm{CD4}^{+}\rm{T})$.

###  II.5 Correlation between chromatin accessibility and mRNA: exclusive mRNA
expressions in ATL

To gain a direct quantification of how chromatin structures are correlated
with gene expressions, we performed an integrated analysis of ATAC-seq and
RNA-seq data for ATL cells.

We used our algorithm to classify the top $M$ peaks into open chromatin
regions, where $M$ is originally determined as $M=64000$, such that the
clustering of the ATAC-seq samples is closest to the appropriate
immunophenotype (Tanaka2020, ). For a comparison, we also used all peaks
outputted from MACS2 as $M=\infty$. We tried to find genes for which at least
one peak from a given set of all peaks is assigned. Note that each peak can be
assigned to zero or more genes.

Here, the RNA-seq data analyzed by edgeR (Robinson2010, ; McCarthy2012, ) was
used to examine the expression of each gene in ATL cells and healthy CD4+T
cells. For details of the calculations, see Material and Methods.

We considered the set of cell types $\mathbb{T}$ as

$\displaystyle\mathbb{T}=\\{\rm{HSC},\rm{CD4}^{+}\rm{T},\rm{CD8}^{+}\rm{T},\rm{NK},\rm{Mono},\rm{ATL}\\},$

where HSC, NK, and Mono are hematopoietic stem cell, natural killer cell, and
monocyte, respectively. We computed the fold change ${\rm FC}_{i}(t,t_{0})$ of
the expression of gene $i$ between type $t,t_{0}\in\mathbb{T}$ as

$\displaystyle{\rm
FC}_{i}(t,t_{0}):=\dfrac{\overline{R}_{i}(t)}{\overline{R}_{i}(t_{0})},$ (1)

where $\overline{R}_{i}(t)$ means the average of the normalized reads count of
the RNA-seq data in gene $i$ over all samples with type $t$. Then, we focused
on the log fold change of the gene expression $P=\log_{2}{\rm
FC}_{i}(t,t_{0})$, where we only take into account the expression of gene $i$
in which at least one peak within the top $M$ peaks from $\hat{g}_{s}^{M}$ was
located for types $t,t_{0}$ of the samples. For details of the calculations,
see Material and Methods.

As shown in Fig. 7, when $M=\infty$, the distribution is close to unimodal for
almost all pairs of cell types. When $M=64000$, the distributions
significantly depend on the pairs of cell types. For example, as shown in Fig.
7b, 7c, and 7e, the distributions related to ATL cases are sharply bimodal-
like. In particular, the pair ATL and CD4+ T cells was found to be distinct
(Fig. 7b); the events around $P=0$ are completely undetected.

Thus, the above observations suggest that ATL cells have exceptionally
distinct structures in terms of the correlation between chromatin
accessibility and mRNA compared with CD4+ T cells.

Figure 8: Reads counts overlapping with the TREC region divided by reads
counts overlapping with the TRA gene region.

###  II.6 Immunophenotypes from chromatin accessibility and mRNA

Table 1: Clustering results of ATAC-seq from ATL cases in terms of immunophenotypes obtained using the method in (Tanaka2020, ). ATAC-seq information for ATL samples is DRR250714 for ATL8, DRR250710 for ATL10, DRR250711 for ATL21, DRR250712 for ATL24, DRR250713 for ATL4, DRR250715 for ATL2, and DRR250716 for ATL5. Other information about the sample labels will be added in next version. Sample | first | second | third | Clinical
---|---|---|---|---
labels |  |  |  | subtype
1 | Ery | Mono | CLP | Acute
2 | Mono | CD4+ T | B | Chronic
3 | CD4+ T | CD8+ T | NK | Chronic
4 | CD4+ T | CD8+ T | B | Acute
5 | Mono | Ery | CD4+ T | Chronic
6 | Ery | Mono | B | Chronic
7 | CD4+ T | CD8+ T | B | Acute
8 | CD4+ T | CD8+ T | NK | Acute
9 | CD4+ T | CD8+ T | B | Acute
10 | CD4+ T | CD8+ T | B | Acute
11 | CD4+ T | CD8+ T | B | Acute
12 | CD4+ T | CD8+ T | B | Acute
13 | CD4+ T | CD8+ T | B | Acute
14 | CD4+ T | CD8+ T | B | Acute
15 | CD4+ T | B | CD8+ T | Acute
16 | CD4+ T | CD8+ T | NK | Acute
17 | Ery | Mono | B | Chronic
18 | CD4+ T | CD8+ T | B | Acute
19 | CD4+ T | B | CD8+ T | Acute
20 | CD4+ T | CD8+ T | B | Chronic
21 | CD4+ T | B | CD8+ T | Acute
22 | CD4+ T | B | CD8+ T | Acute
23 | Ery | CD4+ T | Mono | Acute
24 | CD4+ T | CD8+ T | B | Acute
25 | CD4+ T | CD8+ T | B | Acute
26 | CD4+ T | CD8+ T | B | Acute
27 | CD4+ T | CD8+ T | B | Chronic
28 | CD4+ T | CD8+ T | B | Acute
29 | CD4+ T | CD8+ T | B | Acute

We next examined ATL cases by inferring the past cell status before infection
with HTLV-1 and the current cell status in terms of immunophenotypes compared
with normal hematopoietic cells.

It remains unclear why most ATL cells have an immunophenotype similar to CD4+
memory T cells even though HTLV-1 infects multiple hematopoietic cells
(Furuta2017, ). Specifically, by focusing on the T-cell Receptor Excision
Circles (TREC) region, we computed the number of reads mapped to the TREC
region using ATAC-seq data from ATL cells and $13$ human primary blood cell
types from healthy donors. As shown in Fig. 8, the number of mapped reads in
the TREC region in ATL/HAM cells is closer to that of CD4+ T and CD8+ T cells
compared with other hematopoietic cell types. This finding indicates that
cells before infection had already differentiated into T cells. Moreover, when
combined with our findings about chromatin levels above, the results suggest a
strong sample-dependent opening of the chromatin occurs in infected CD4+ T
cells during the development of ATL.

As mentioned above, the overall chromatin structures vary greatly among ATL
cases and are correlated with the gene expressions. This phenomenon presumably
affects the immunophenotype. To investigate the effect of large-scale
variations of chromatin structures on the immunophenotype, we used our
algorithm to quantitatively evaluate differences between ATL cells and
hematopoietic cells from healthy donors (Tanaka2020, ).

Using this algorithm, we calculated the Hamming distances between the peak-
based binarized genome of ATL cells and of hematopoietic cells from healthy
donors. For this purpose, we used ATL samples and $77$ ATAC-seq datasets from
$13$ human primary blood cell types. As summarized in Table 1, the majority of
ATL samples are close to CD4+ T cells, as expected by the above analysis about
the past cell status. We also found that the ATAC-seq patterns of some ATL
cases are close to myeloid cells such as erythroid cells and monocytes.

To ascertain whether the mRNA expression in the ATL cells reflects the
characteristics of myeloid cells, we analyzed the RNA-seq data from healthy
CD4+ T cells and HTLV-1-infected CD4+ cells from 4 ATL cases (samples 8, 10,
21, 24 in Table 1). We used the condition $\log_{2}\rm{FC}>1$ and $p$-value
$<0.01$ to identify up-regulated genes and found two candidates: CD71 (TFRC),
which is ubiquitously expressed by erythroid precursors (Dong2011, ;
Marsee2010, ), and KLF4, which is highly expressed in myeloid cells and
essential for monocyte differentiation (Feinberg2007, ).

###  II.7 Chromatin-based systematic selection of key genes in ATL cases

(a)

(b)

Figure 9:  (a) The number of selected genes over increments of $8000$ for $M$
from $8000$ to $80000$. (b) The width of the selected genes.

To identify genes that are specially expressed in ATL cells, we investigated
whether there are characteristic open chromatin regions in gene coding regions
that are common to all ATL cells but not to hematopoietic cells derived from
healthy individuals.

We therefore compared the chromatin accessibility between $29$ ATL samples and
$77$ ATAC-seq datasets from $13$ human primary blood cell types. We applied
our algorithm to detect such key genes, using $M$ as a parameter for
clustering the ATAC-seq data corresponding to immunophenotypes (Tanaka2020, ).
First, we defined a subset of regions of top $M$ peaks, where all ATL samples
had peaks but no $13$ human primary blood cell types from healthy donors did.
We call this subset ATL-specific open regions. Second, we detected the genes
in which at least one of ATL-specific open region was located for a given $M$.

(a)

(b)

(c)

Figure 10: (a,b) Histogram of ATAC-seq reads around the TLL1 region. (c)
Expression of TLL1 relative to GAPDH (Glyceraldehyde-3-phosphate
dehydrogenase) for ATL sample $s$. The relative expression of sample $s$ is
defined as $2^{\Delta C_{s}^{0}}$, where $\Delta C_{s}^{0}$ is computed using
the $\Delta\Delta$ Ct method; for details, see (19). N.D. stands for no
detection of TLL1.

As shown in Fig. 9, the number of key genes was $0$ at $M$=$2000$, $1$ at
$M=4000$, $3$ to $4$ between $M=8000$ and $56000$, and then dropped back to
$1$ for values larger than $M=64000$. Among these genes, we picked those that
stayed in intervals of M $\geq 8000$. Concretely, TLL1 (Tolloid-like-1) gene
appeared from $M=8000$-$16000$, EVC (Ellis-van Creveld) gene and CRMP1
(Collapsin response mediator protein 1) gene appeared from $M=16000$-$48000$,
TNRC6A (Trinucleotide Repeat Containing Adaptor 6A) gene appeared from
$M=32000$-$48000$, and UST (Uronyl-2-sulfotransferase) gene appeared from
$M=64000$. As a reference in Fig. 10a, we show the locations of the ATAC-seq
reads around TLL1 gene. Note that the detected genomic region for EVC and
CRMP1 is the same because the two genes overlap. Thus, by our selection, the
candidates determined as key genes are few.

Consistently, EVC has been reported to be overexpressed in ATL and plays an
important role in cellular Hedgehog activation (Takahashi2014, ). UST was also
highly expressed in ATL cases, though the relationship between the function of
UST and ATL has not been explicitly clarified (Sasaki2005, ). On the other
hand, to the best of our knowledge, TLL1 and TNRC6A have not been studied in
this context. While TLL1 is known to be necessary for normal septation and
positioning of the heart (Clark1999, ), a more recent report found that it is
associated with the development of hepatocellular carcinoma after the
eradication of HCV (Matsuura2017, ). Further, TLL1 is a member of the BMP1/TLD
(bone morphogenetic protein1/tolloid)-like proteinase family; BMP1 controls
latent TGF-${\beta}$ activation via the cleavage of LTBP1 (latent TGF-$\beta$
binding protein-1) (Ge2006, ), and TGF-$\beta$ plays important roles in cancer
progression.

Thus, we picked up TLL1 as a promising candidate among genes expressing ATL-
specific functions.

(a)

(b)

Figure 11:  (a) Schematic picture of pre-pro TGF-$\beta$ and mature
TGF-$\beta$. (b) TGF-$\beta$ activation was measured by 3TP-Lux protein
activation, which depended on the TLL1 isoform type. The $p$-value is 0.0028
($**$) for Case 2 vs. Case 3 and $0.019$ ($*$) for Case 2 vs. Case 4 (t-test).

###  II.8 TLL1 can strongly regulate TGF-$\beta$

We next considered the gene expression of TLL1 in ATL cells and TLL1’s effect
on the maturation process of TGF-$\beta$ in a HepG2 line, which is a
TGF-$\beta$ responsive cell line(Westerhausen1991, ).

First, as shown in Fig. 10c, we performed real-time PCR that showed TLL1 is
expressed in ATL cases but not in the peripheral blood mononuclear cells of
healthy donors. This result is consistent with the Human Protein Atlas, which
shows that TLL1 mRNA is not detected in most adult tissues including immune
cells or any hematopoietic lineage.

By analyzing the RNA-seq data, we also confirmed that TLL1 mRNA is not
expressed in normal hematopoietic cells, but it was expressed in all examined
ATL cases. The same was not true for TNRC6A gene, which, according to RNA-seq
data, did not show any systematic expression change between ATL cases and
normal hematopoietic cells.

Next, we considered the relationship between TLL1 and the maturation process
of TGF-$\beta$. TLL1 has two mRNA isoforms: TLL1 isoform 2 lacks many exons
from the 3’ end of TLL1 isoform 1. Thus, we asked if both isoforms regulate
TGF-${\beta}$ in a manner similar to BMP1.

To investigate this possibility, as shown in Fig. 11a, we performed a
luciferase assay using the pre-mature form of TGF-${\beta}$ co-expressed with
either TLL1 isoform 1 or 2 in a HepG2 cell line. As shown in Fig. 11b, we
found that compared to the sample without TLL1 (Case 2), TLL1 isoform 1 (Case
3) activates the pre-mature form of TGF-$\beta$ for maturation, but TLL1
isoform 2 (Case 4) represses the maturation. It should be noted that the
difference of luciferase activity between isoform 1 (Case 3) and isoform 2
(Case 4) with pre-mature TGF-$\beta$ approximated the difference between the
TLL1-less condition with pre-mature TGF-$\beta$ (Case 2) and without pre-
mature TGF-$\beta$ (Case 1). Thus, the results suggest that TLL1 is able to
strongly regulate TGF-${\beta}$ depending on the isoform expression ratio.

##  III Discussion

In this paper, we statistically characterized the anomalous chromatin
accessibility and gene expression of HTLV-1-infected cells and healthy CD4+ T
cells at the whole genome level.

Our analysis suggests that compared to healthy CD4+ T cells, ATL cells have
the following properties: the chromatin accessibility increases near TSSs, a
higher frequency of larger sample-to-sample fluctuations at the whole genome
level, and a higher frequency of intermediate sample-to-sample fluctuations in
gene coding regions. Consistently, histone-related genes were up-regulated.
The expression of the genes systematically selected by the chromatin
accessibility was found to be distinct from healthy CD4+ T cells but not other
hematopoietic cell types that we have studied. Further, whereas the
immunophenotype determined by the systematically selected open chromatin
regions was classified to be near CD4+ T cells for most samples, some samples
were classified as myeloid cells.

Based on the above integrative analysis of chromatin accessibility and gene
expression, we found that there are chromatin regions which are open in all
the ATL cases but closed in all the analyzed samples of the 13 hematopoietic
cell types derived from healthy donors. One of the genes overlapping with the
chromatin regions that satisfy such conditions is TLL1, which was
experimentally shown to have a large potential to regulate TGF-$\beta$.

Contrary to ATL cases, the statistics of the chromatins in HAM cells resembled
those of healthy CD4+ T cells, including sample-to-sample fluctuations. This
observation implies that there is a certain sample-independent trend in the
chromatin structure of HAM cases.

It should be noted that we were unable to analyze large number of samples due
to the difficulty in obtaining samples for given clinical conditions. Thus,
some quantities, such as frequencies shown in Fig. 6 and Fig. 7, were not
estimated with statistical validity. To validate the hypothesis about the
uniqueness of ATL cells across different scales from chromatin and
transcription to immunophenotypes, more samples should be used in future
studies.

Our finding about ATL samples might motivate us to consider a rather general
relationship between increased chromatin accessibility and the onset of
leukemia. Indeed, in a previous study of Acute Myeloid Leukemia (AML), it was
reported that mutations in cohesin genes increase chromatin accessibility,
which controls the activity of transcription factors leading to leukemogenesis
(Mazumdar2015, ). It was also reported that HMGN1 amplification is also
associated with increased chromatin accessibility; It confers a
transcriptional and chromatin phenotype associated with stem cells and
leukemia (Cabal-Hierro2020, ). It remains for future studies to check how the
relationship between the increased chromatin accessibility and the onset of
leukemia can be generalized beyond the above cases such as ATL and AML.

Let us discuss our additional findings in the following, though they are
preliminary results. First, to compare ATL cases with another type of leukemia
in terms of immunophenotypes, we analyzed the chromatin structure of cutaneous
T-cell lymphoma (CTCL) using the ATAC-seq data of CTCL (GSE85853) (Tanaka2020,
), which is reported to have a clinical and histopathological phenotype
similar to that of ATL (Qu2017, ). As shown in Table 2, we found that the
chromatin structures in some CTCL cases are closer to myeloid cells rather
than CD4+ T cells, though CTCL is conventionally classified as T-cell
leukemia. Especially, in the cases of Patient 11 and 60, both of whom are
romidepsin responders, the chromatin structure changed from myeloid cell-like
to CD4+ T-like after drug treatment on the 0th day. Although the molecular
mechanisms underlying this process are not understood, these results suggest
that the change in immunophenotype reflects a molecular response in treatment.
This observation could shed light on finding better therapeutic targets and
prediction of drug response.

Table 2: Clustering results of ATAC-seq from 9 CTCL cases in terms of immunophenotype as a function of time using the method in (Tanaka2020, ). The Histone deacetylase inhibitor (HDACi) was romidepsin. $+$ and $-$, positive response and negative response, respectively (Qu2017, ). SRR | first | second | third | Patient tag | HDACi
---|---|---|---|---|---
Number |  |  |  | & Time tag | responder
4044872 | Ery | Mono | CLP | Patient-11 on 0-th day |
4044873 | CD4+ T | CD8+ T | Ery | on 0-th day |
4044874 | CD4+ T | CD8+ T | NK | on 7-th day | $+$
4044875 | CD4+ T | CD8+ T | NK | on 7-th day |
4044876 | CD4+ T | CD8+ T | B | on 14-th day |
4044877 | Ery | CD4+ T | Mono | on 14-th day |
4044878 | Ery | Mono | CLP | Patient-20 on 7-th day | $-$
4044879 | Ery | Mono | CLP | Patient-39 on 0-th day |
4044880 | Ery | Mono | CLP | on 0-th day | $-$
4044881 | Ery | Mono | CLP | on 7-th day |
4044882 | Ery | Mono | CLP | on 7-th day |
4044885 | CD4+ T | CD8+ T | B | Patient-59 on 7-th day | $+$
4044886 | CD4+ T | CD8+ T | B | on 7-th day |
4044887 | Ery | Mono | CD4+ T | Patient-60 on 0-th day | $+$
4044888 | CD4+ T | B | CD8+ T | on 7-th day |
4044889 | Ery | Mono | CLP | Patient-61 on 0-th day | $+$
4044890 | CD4+ T | CD8+ T | B | Patient-62 on 0-th day | $+$
4044891 | CD4+ T | Ery | CD8+ T | on 0-th day |
4044892 | CD4+ T | CD8+ T | NK | Patient-1366 on 0-th day | $+$

Second, looking at the ATAC-seq data further, a footprint analysis for the
identification of differential motifs revealed that ETS1, IRF2, and RUNX2 had
deeper footprints and higher DNA accessibility at the flanking locations of
their motifs in healthy CD4+ T cells, while NRF1, KLF4, and KLF9 had deeper
footprints in ATL cells (Fig. 12) (Li2019, ). These observations suggest that
transcription factors such as NRF1, KLF4, and KLF9 play an important role in
ATL.

Third, to understand TLL1’s function, we conducted experimental studies on
MT-2, an ATL cell line that does not express TLL1. We prepared three samples
of MT-2: one transduced with an empty vector, another with TLL1 isoform 1, and
the third with TLL1 isoform 2. We analyzed the RNA-seq data from the three
samples by computing the gene expression ratio between them.

Table 3: Relative expression (RE) of genes in TLL1-transduced MT-2 cells over MT-2 cells transduced with an empty vector. Type 1(2) corresponds to MT-2 cells transduced with TLL1-isoform 1 (isoform 2). Genes were picked up if the RE of type 1 was larger than $11$ and the RE of type 2 was smaller than $3$ or if the RE of type 1 was smaller than $1$ and the RE of type 2 was larger than $5$. Gene symbol | RE of type 1 | RE of type 2
---|---|---
CCL3 | 14.53 | 2.94
CCR6 | 11.70 | 2.85
MIR155 | 11.30 | 2.44
POSTN | 11.06 | 2.49
HBG2 | 0.058 | 15.00
HBB | 0.15 | 12.29
HBA2 | 0.14 | 5.79
HBA1 | 0.09 | 5.04

Among the $22963$ genes, we picked up the genes which have nonzero reads count
in MT-2 transduced with an empty vector. A part of such genes showed
significantly altered expressions depending on the type of isoform. As shown
in Table 3, four genes were significantly up-regulated when TLL1-isoform 1 was
transduced: CCR6 is related to the regulation of Treg migration (Yamazaki2008,
), microRNA-155 (MIR155) modulates Treg cell differentiation and its
expression is up-regulated in HTLV-1 transformed cells (Yao2012, ;
Pichler2008, ; Watanabe2017, ), the chemokine CCL3 regulates myeloid
differentiation (Staversky2018, ), and POSTN has been reported to be involved
in TGF-${\beta}$ activation (Sidhu2010, ). As for the four genes significantly
up-regulated when TLL2-isoform2 was transduced, all were globin genes. These
findings reiterate the dependence of TLL’s function on its isoforms in ATL
cells.

As an additional note, Fig. 10b suggests that the chromatin regions around the
TLL1 loci tend to be open also for the HAM cases. This observation indicates
that open chromatin regions around TLL1 are not the only cause of leukemia
onset. Rather, it suggests that open chromatin regions are potentially caused
by the infection itself and related to the latent period or expansion of the
virus.

This study is the first comprehensive analysis of open chromatin structures in
ATL samples. The findings will deepen understanding of the ATL pathogenesis.

(a) ETS1

(b) IRF2

(c) RUNX2

(d) NRF1

(e) KLF4

(f) KLF9

Figure 12: A footprint analysis of transcription factors. (a)-(f) Distances
from motifs (horizontal axis) vs. the averaged number of reads over all parts
of a given motif (vertical axis) outputted from HINT-ATAC (Li2019, ).

##  IV Methods and Materials

Figure 13: Pipeline of the data processing for ATAC-seq.

###  IV.1 Sequencing sample preparation

Peripheral blood mononuclear cells from ATL patients, HAM patients, and HTLV-1
carriers were thawed and washed with PBS containing 0.1% BSA. To discriminate
dead cells, we used a LIVE/DEAD Fixable Dead Cell Stain Kit (Invitrogen). For
cell surface staining, cells were stained with APC anti-human CD4 (clone:
RPA-T4) (BioLegend) and anti-SynCAM (TSLC1/CADM1) mAb-FITC (MBL) antibodies
for 30 minutes at 4 ℃ followed by washing with PBS. HTLV-1-infected cells
(CAMD1+ and CD4+) were purified by sorting with a FACS Aria (Beckman Coulter)
to reach 98–99% purity. Data were analyzed using FlowJo software (Treestar).
Soon after the sorting, 10000-50000 HTLV-1-infected cells were centrifuged and
used for ATAC-seq. Total RNA was isolated from the remaining cells using an
RNeasy Mini Kit (Qiagen). Library preparation and high-throughput sequencing
were performed at Macrogen Inc. (Seoul, Korea). The diagnostic criteria and
classification of the clinical subtypes of ATL were performed as previously
described (Shimoyama1991, ). 77 ATAC-seq datasets from 13 human primary blood
cell types were obtained from the Gene Expression Omnibus (GEO) with accession
number GSE74912 (Corces2016, ), and RNA-seq datasets of CD4+ T and CD8+ T
cells from healthy donors were obtained from GSE74246 (Corces2016, ). The RNA-
seq data for ATL samples can be downloaded from DBBJ (DNA Data Bank of Japan)
with accession number DRR250721 for ATL8, DRR250717 for ATL10, DRR250718 for
ATL21, and DRR250719 for ATL24. The ATAC-seq data for ATL samples can be
downloaded from DBBJ with accession number DRR250714 for ATL8, DRR250710 for
ATL10, DRR250711 FOr ATL21, and DRR250712 for ATL 24.

###  IV.2 Pre-processing of ATAC-seq

High-throughput sequencing provides a set of reads as output. ATAC-seq reads
were aligned using BWA version 0.7.16a with default parameters. SAMtools was
used to convert SAM files into compressed BAM files and sort the BAM files by
chromosome coordinates. PICARD software (v1.119)
(http://broadinstitute.github.io/picard/) was used to remove PCR duplicates
using the MarkDuplicates options. Reads with mapping quality scores less than
30 were removed from the BAM files. For peak calling, MACS2 (v2.1.2) software
was used with option \--nomodel --nolambda --keep-dup all -p 0.01. ATAC-seq
tracks were visualized using Integrative Genomics Viewer (IGV), and
footprinting analysis was performed using HINT-ATAC (Li2019, ). Note that the
paired-end output of the sequence was used to reconstruct the fragments, where
paired two reads correspond to both ends of a fragment.

###  IV.3 Pre-processing of RNA-seq

RNA-seq data were aligned to human reference genome hg19 using STAR 2.6.0c
with the \--quantMode GeneCounts function (Dobin2013, ). RNA-seq data analysis
was performed using edgeR, where the reads counts of the RNA-seq data were
normalized using TMM normalization (cite-key, ) to be converted into pseudo
reads counts. Let $n_{i}^{0}(s)$ be the reads count of the RNA-seq data for
each gene $i$ for a given cell sample $s$ and $N(s):=\sum_{i}n_{i}^{0}(s)$ be
the total reads count over all genes. Using the TMM normalization with
$n_{i}^{0}(s)$ for gene $i$ and sample $s$ as the input data, one can obtain
the normalization factor $r(s)$ for a given sample $s$ using the command
calcNormFactors of edgeR. After acquiring $r(s)$, the pseudo reads count
$n_{i}(s)$ is calculated using the command estimateCommonDisp of edgeR, which
we used as the starting point of the RNA-seq data analysis in the main text.

An additional analysis was done to evaluate robustness. We computed the
geometric mean of $N(s)r(s)$ as
$N_{0}=\left(\prod_{s\in\mathbb{S}}r(s)N(s)\right)^{1/|\mathbb{S}|}$, where
$\mathbb{S}$ is the set of samples and $|\mathbb{S}|$ is the number of
elements of $\mathbb{S}$. We then checked whether the pseudo reads count
$n_{i}(s)$ is close to a normalized reads count
$n_{i}^{\prime}(s):=n_{i}^{0}(s)\frac{N_{0}}{r(s)N(s)}$ for sample $s$. We
found that the maximum deviation between $n_{i}^{\prime}(s)$ and $n_{i}(s)$ is
smaller than $5$ over all genes in the case of CD4+ T vs. Mono in Fig. 7a. In
this case, the effects from the differences between $n_{i}^{\prime}(s)$ and
$n_{i}(s)$ are quite small except for quantities related to genes with almost
zero reads. Therefore, even if we use $n_{i}^{\prime}(s)$ as the starting
point of the analysis instead of the pseudo reads count $n_{i}(s)$,
qualitatively the same conclusion as that with $n_{i}(s)$ is expected to be
obtained.

A PCA was done using the covariance matrix of $\log_{10}(n_{i}(s)+1)$, where
the first and second principal components were calculated using the prcomp
command with option scale=FALSE in R.

###  IV.4 Cell lines and clinical samples

All ATL cell lines were cultured in RPMI 1640 medium supplemented with 10% FBS
and antibiotics. HepG2 was cultured in DMEM. To construct MT-2 cells stably
expressing TLL1 isoform 1 or 2, the coding sequence of each isoform was
transduced using a lentivirus vector constructed as described in subsection
IV.6.

###  IV.5 Real-time PCR

cDNA products were analyzed by real-time PCR using PowerUp SYBR Green Master
Mix and StepOnePlus Real-Time PCR System (Applied Biosystems) according to the
manufacturer’s instructions. Primer sequences for the GAPDH gene have been
described previously (Ponchel2003, ), and primer sequences for the TLL1 gene
were 5- TTGTTTTCTACGGGGAGCTATGG-3 and 5- ATATCGCCCCAAAATACAGCG-3. The relative
quantification was calculated according to the method described in Applied
Biosystems ABI prism 7700 SDS User Bulletin #2. Note that the ATL sample used
for this experiment is not listed  in Table 1.

###  IV.6 Lentiviral vector construction and transfection of recombinant
lentivirus

The coding region of TLL1 isoform 2 was synthesized using gBlocks Gene
Fragment (Integrated DNA Technolofies), which was used as the template for
synthesizing TLL1 isoform 1 by PCR amplification. TLL1 isoform 1 and 2
fragments were subcloned into pCS2-EF-MCS (gift from H. Miyoshi, RIKEN
Bioresource Center). An empty vector that expresses only hrGFP was used as the
control for the lentiviral transduction. 293T cells at 80% confluence in a
10-cm dish were co-transfected with 10 ${\mu}$g lentivirus vector, 10 ${\mu}$g
psPAX2, 5 ${\mu}$g pMD2.G, and PEI (Polyethylenimine). 48 hours after the
transfection, supernatant containing the virus was collected and concentrated
by ultracentrifugation. MT-2 cells were transfected with the lentivirus, and
two weeks after the transduction, GFP-positive cells were purified by sorting
with FACS Canto. RNA was isolated using the Qiagen RNeasy Mini Kit and then
used for the RNA-seq analysis.

###  IV.7 Luciferase assay

The coding region of human TGF-${\beta}$, whose length is 1173 bp, was
generated by PCR amplification and subcloned into a pFUSE-hIgG1-Fc2 vector.
HepG2 cells were plated on 12-well plates at ${1\times 10^{5}}$ cells per
well. After 24 hours, the cells were transfected with 50 ng/well of luciferase
reporter plasmid (p3TP-Lux) (Zhao2011, ), 5 ng/well of Renilla luciferase
control vector (phRL-TK) together with 35 ng/well of TLL1 expression plasmid,
and 35 ng/well of TGF-${\beta}$-expressing plasmid or empty vector. Plasmids
were transfected using TransIT-LT1 (Mirus) according to the manufacturer’s
instructions. After 48 hours, the cells were collected, and luciferase
activities were measured using the Dual-Luciferase Reporter Assay Kit
(Promega). Relative luciferase activity was calculated as the ratio of firefly
to Renilla luciferase activity. Three independent experiments, each with
triplicate transfections, were performed, and typical results are shown.

###  IV.8 Explicit definitions of the computed quantities

We explicitly define quantities discussed in the main text. First, we assume
that the set of reads from DNA of an ATAC-seq sample and the reads count for
each gene from an RNA-seq sample are given. Also, a set of fragments for the
ATAC-seq sample is also given by using a pair of reads; both ends of a
fragment correspond to a pair of two reads. The positions of TSSs and coding
regions of all genes were obtained from the human genome (hg19) as a set of
intervals on the genome. Therefore, a read from the ATAC-seq data is interval
$[x_{1},x_{2}]$ on the genome, which corresponds to a region including an edge
of a fragment. A fragment has length $\ell$ and location $x$ as the mid-point
of the two edges on the genome (Tanaka2020, ). The reads from the RNA-seq data
provide the read count $n_{i}(s)$ for each gene $i$, where $s\in\mathbb{S}$ is
the sample index. We denote by $\mathbb{S}_{\nu}$ the set of all analyzed
samples with type $\nu$.

####  IV.8.1 The normalized number of reads in Fig. 3a:

Let us consider the number $\rho_{s}(z)$ of ATAC-seq reads from a sample
$s\in\mathbb{S}$ located on position $z$ from the nearest TSS. Then, we take
the sample average among type $\nu$ as

$\displaystyle\overline{\rho}_{\nu}(z):=\dfrac{1}{|\mathbb{S}_{\nu}|}\sum_{s\in\mathbb{S}_{\nu}}\rho_{s}(z),$
(2)

where $\nu\in\\{\mathrm{CD4^{+}T},\mathrm{HAM},\mathrm{ATL}\\}$. In Fig. 3a,
we plot the normalized quantity $\widetilde{\rho}_{\nu}(z)$ obtained after
dividing $\overline{\rho}_{\nu}(z)$ by the value at the TSS ($z=0$) such that

$\displaystyle\widetilde{\rho}_{\nu}(z):=\dfrac{\overline{\rho}_{\nu}(z)}{\overline{\rho}_{\nu}(0)}.$
(3)

####  IV.8.2 The averaged number of fragments in Fig. 3b:

Let us consider the number of fragments $\phi_{s}(z,\ell)$ from sample
$s\in\mathbb{S}$ satisfying the following two conditions: (1) their centers
are located at $z$, and (2) they have length $\ell$.
$\overline{\phi}_{\nu}(z,\ell)$ describes the sample average among type $\nu$
as

$\displaystyle\overline{\phi}_{\nu}(z,\ell):=\frac{1}{|\mathbb{S}_{\nu}|}\sum_{s\in\mathbb{S}_{\nu}}\phi_{s}(z,\ell),$
(4)

where $\nu\in\\{\mathrm{CD4^{+}T},\mathrm{HAM},\mathrm{ATL}\\}$. In Fig. 3b,
we plot the histogram $F^{\Delta,\xi}_{\nu}(z,\ell)$ for bin width $\Delta$
and $\xi$ for $z$ and $\ell$, respectively, as

$\displaystyle F^{\Delta,\xi}_{\nu}(z,\ell):=\sum_{z-\Delta/2\leq
z^{\prime}<z+\Delta/2}\sum_{\ell\leq\ell^{\prime}<\ell+\xi}\overline{\phi}_{\nu}(z^{\prime},\ell^{\prime}).$
(5)

####  IV.8.3 The reference set of peaks:

To analyze open chromatin regions, we used MACS2 with the input of reads from
the ATAC-seq data. Concretely, we used MACS2 with the option \--nomodel
--nolambda --keep-dup all -p $0.01$, which corresponds to $p_{G}=10^{-2}$
(Tanaka2020, ). This algorithm outputs the collection of peaks $\hat{g}_{s}$
for a given sample $s$ as candidates of open chromatin regions, which can be
described as

$\displaystyle\hat{g}_{s}:=((\gamma_{k},\alpha_{k},\beta_{k}),p_{k})_{k\geq
1},$ (6)

where $\gamma_{k}$ is the chromosome number, $\alpha_{k}$ is the starting
point, and $\beta_{k}$ is the ending point in terms of genome position with
$p_{k}$ as the $p$-value of the $k$-th peak. As in (Tanaka2020, ), $p_{k}\leq
p_{k^{\prime}}$ for $k<k^{\prime}$. In particular, we consider the set of the
top $M$ peaks and denote it by

$\displaystyle\hat{g}_{s}^{M}:=\begin{cases}((\gamma_{k},\alpha_{k},\beta_{k}),p_{k})_{k=1}^{M}&(\text{if
$|\hat{g}_{s}|\geq M$}),\\\ \hat{g}_{s}&(\text{otherwise}).\end{cases}$ (7)

Next, we concatenate the data of all reads from all ATAC-seq samples with cell
type $\nu\in\\{\mathrm{CD4^{+}T},\mathrm{HAM},\mathrm{ATL}\\}$. Then, we
randomly extract 100 million reads from the concatenated data for type $\nu$
as the input for the MACS2 algorithm to obtain the collection of peaks,

$\displaystyle g_{\nu}:=((\gamma_{k},\alpha_{k},\beta_{k}),p_{k})_{k\geq 1}.$
(8)

Using a coalescing process of $g_{ATL},g_{HAM},g_{CD4^{+}T}$, we construct a
new reference set of peaks $g_{0}$ as follows. Operationally, the coalescing
of two peaks is done as follows. If two peaks $(\gamma,\alpha,\beta)$ and
$(\gamma^{\prime},\alpha^{\prime},\beta^{\prime})$ for
$\gamma=\gamma^{\prime}$ satisfy
$\alpha^{\prime}\leq\alpha\leq\beta^{\prime}$, the two peaks become one peak
as $(\gamma,\alpha^{\prime},\max\\{\beta,\beta^{\prime}\\})$. This operation
is repeated for a newly obtained set of peaks until no more coalescing
processes occur.

####  IV.8.4 The length of overlapped peaks in Fig. 4, 5::

To quantify the similarity between two collection of peaks, $g_{0}$ and
$\hat{g}_{s}^{M}$, first, we fix a set of peaks $\mathbb{L}$ as any of (1) the
set $\mathbb{G}$ of all peaks in $g_{0}$ overlapping gene coding regions, (2)
the set $\mathbb{G}^{c}$ of all peaks in $g_{0}$ corresponding to non-coding
regions, and (3) the union $\mathbb{G}\cup\mathbb{G}^{c}$ of two sets
$\mathbb{G},\mathbb{G}^{c}$. Note that for a given peak, its center was
calculated by the command of annotatePeaks.pl in the HOMER algorithm and used
to judge whether the peak joins $\mathbb{G}$ or $\mathbb{G}^{c}$ Heinz2015 .

Then, we focus on the length of the overlapped peaks
$O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})$, which is the number of base pairs
in a peak of $\hat{g}_{s}^{M}$ inside the $k$-th peak $(\alpha_{k},\beta_{k})$
of $\mathbb{L}\subset g_{0}$. Then, we compute the average and variance of
$O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})$ as follows:

$\displaystyle\overline{O}_{k}(\mathbb{L},\mathbb{S}):=\dfrac{1}{|\mathbb{S}|}\sum_{s\in\mathbb{S}}O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M}),$
(9) $\displaystyle
V_{k}(\mathbb{L},\mathbb{S}):=\dfrac{1}{|\mathbb{S}|}\sum_{s\in\mathbb{S}}\big{(}O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})-\overline{O}_{k}(\mathbb{S},\mathbb{L})\big{)}^{2}.$
(10)

We set $M=64000$ as the provisionally optimal number for immunophenotype
classification (Tanaka2020, ).

The following functions describe the frequency of the average and variance of
$O_{k}^{\mathbb{L}}(g_{0},\hat{g}_{s}^{M})$:

$\displaystyle\rho^{(1)}_{\mathbb{L},\mathbb{S}}(O):=\sum_{k\in\mathbb{L}}\delta(O,\overline{O}_{k}(\mathbb{L},\mathbb{S})),$
(11)
$\displaystyle\rho^{(2)}_{\mathbb{L},\mathbb{S}}(V):=\sum_{k\in\mathbb{L}}\delta(V,V_{k}(\mathbb{L},\mathbb{S})),$
(12)

where $\delta(a,b)=1$ for $a=b$, otherwise $0$.

Lastly, in Fig. 4 and Fig. 5, the histograms are defined as

$\displaystyle F^{(1)}_{\mathbb{L},\mathbb{S}}(O;\Delta):=\sum_{O-\Delta/2\leq
O^{\prime}<O+\Delta/2}\rho^{(1)}_{\mathbb{L},\mathbb{S}}(O^{\prime}),$ (13)
$\displaystyle F^{(2)}_{\mathbb{L},\mathbb{S}}(V;\Delta):=\sum_{V-\Delta/2\leq
V^{\prime}<V+\Delta/2}\rho^{(2)}_{\mathbb{L},\mathbb{S}}(V^{\prime}).$ (14)

####  IV.8.5 Fold change of selected gene expression in Fig. 7:

For the RNA-seq data, we consider the set of cell types $\mathbb{T}$ as

$\displaystyle\mathbb{T}=\\{\rm{HSC},\rm{CD4}^{+}\rm{T},\rm{CD8}^{+}\rm{T},\rm{NK},\rm{Mono},\rm{ATL}\\}.$
(15)

We compute the fold change in gene $i$, ${\rm FC}_{i}(t,t_{0})$ between type
$t,t_{0}\in\mathbb{T}$ as

$\displaystyle{\rm
FC}_{i}(t,t_{0}):=\dfrac{\overline{R}_{i}(t)}{\overline{R}_{i}(t_{0})},$ (16)

where $\overline{R}_{i}(t)$ is the average normalized expression $n_{i}(s)$ of
the RNA-seq data for gene $i$ over all samples with type $t$. Here, we
consider only the genes where $\log_{2}{\rm FC}_{i}(t,t_{0})$ is well-defined;
in other words, we consider the genes $i$ satisfying the following conditions:
(1) a peak of sample type $t$ that intersects the coding region of the gene,
and (2) the same condition for type $t_{0}$, where the corresponding $k$-th
peak satisfies $k\leq M$. We denote $\mathbb{G}_{t,t_{0}}(M)$ as the set of
all the genes satisfying these conditions. Then, we focus on the following
function quantifying the frequency of the log-fold change

$\displaystyle\rho_{t,t_{0}}^{M}(P):=\dfrac{1}{|\mathbb{G}_{t,t_{0}}(M)|}\sum_{i\in\mathbb{G}_{t,t_{0}}(M)}\delta(P,\log_{2}{\rm
FC}_{i}(t,t_{0})).$ (17)

In Fig. 7, we plot the histogram of the frequency of the log-fold change

$\displaystyle F_{t,t_{0}}^{M}(P;\Delta):=\sum_{P-\Delta/2\leq
P^{\prime}<P+\Delta/2}\rho_{t,t_{0}}^{M}(P^{\prime}).$ (18)

####  IV.8.6 $\Delta\Delta C_{t}$ method in Fig. 10c:

To obtain the threshold cycle for real-time PCR of an mRNA sample (see ”Real-
time PCR” for details), for sample $s\in\mathbb{S}$, we denote $C_{s}^{TLL1}$
as the threshold cycle for gene TLL1 and $C_{s}^{GAPDH}$ for the gene GADPH.
Then, we define the difference $\Delta C_{s}:=C_{s}^{TLL1}-C^{GAPDH}_{s}$ and
consider the normalized difference as

$\displaystyle\Delta C_{s}^{0}:=\Delta C_{s}-\min_{s\in\mathbb{S}}\Delta
C_{s}.$ (19)

### Data availability

All ATAC-seq and RNA-seq data needed to reproduce this study will be deposited
at the DNA Data Bank of Japan (DDBJ) under the accession number XXXXX.

### Ethics approval and consent to participate

Experiments using clinical samples were conducted according to the principles
expressed in the Declaration of Helsinki and approved by the Institutional
Review Board of Kyoto University (permit numbers G310 and G204). ATL patients
provided written informed consent for the collection of samples and subsequent
analysis.

### Author’s contributions

A. Tanaka: Conceptualization, NGS sample preparation, NGS data analysis,
investigation, performing experiments, project administration, generating
figures and tables, funding acquisition, and writing original draft. J.I.
Yasunaga: Collecting clinical samples, data investigation, funding acquisition
and experimental advices. H. Ohta and Y. Ishitsuka: Data investigation,
methodology, generating figures, and writing original draft. C. Onishi and H.
Tanaka: Assisting plasmid preparation, experiments and analyses. N.
Takenouchi, M. Nakagawa and K. Koh: Collecting clinical samples. A. Fujimoto:
Assisting NGS data analysis. M. Matsuoka: Collecting clinical samples,
supervision, funding acquisition, project administration, and writing original
draft. All authors participated in discussions and interpretation of the data
and results.

###### Acknowledgements.

We thank P. Karagiannis for proofreading the manuscript and many valuable
comments. This research was supported by JSPS KAKENHI Grant Numbers JP19K16740
(AT), JP18J40119 (AT), XXXX (MM), and XXXX (JiY) and by a grant from the Naito
Foundation (AT).

## References

  * (1) Anthony SJ, Epstein JH, Murray KA, Navarrete-Macias I, Zambrana-Torrelio CM, Solovyov A, et al.  A Strategy To Estimate Unknown Viral Diversity in Mammals.  mBio. 2013;4(5):e00598–13.  doi:10.1128/mBio.00598-13.
  * (2) Krump NA, You J.  Molecular mechanisms of viral oncogenesis in humans.  Nature Reviews Microbiology. 2018;16(11):684–698.  doi:10.1038/s41579-018-0064-6.
  * (3) Gessain A, Cassar O.  Epidemiological Aspects and World Distribution of HTLV-1 Infection.  Frontiers in Microbiology. 2012;3.
  * (4) Matsuoka M, Jeang KT.  Human T-cell leukaemia virus type 1 (HTLV-1) infectivity and cellular transformation.  Nature Reviews Cancer. 2007;7(4):270–280.
  * (5) Bangham CRM.  HTLV-1 persistence and the oncogenesis of adult T cell leukemia/lymphoma.  Blood. 2023;.
  * (6) Yamagishi M, Fujikawa D, Watanabe T, Uchimaru K.  HTLV-1-Mediated Epigenetic Pathway to Adult T-Cell Leukemia-Lymphoma.  Front Microbiol. 2018;9:1686.
  * (7) Satou Y, Yasunaga JI, Yoshida M, Matsuoka M.  HTLV-I basic leucine zipper factor gene mRNA supports proliferation of adult T cell leukemia cells.  Proceedings of the National Academy of Sciences of the United States of America. 2006;103(3):720–725.
  * (8) Tanaka-Nakanishi A, Yasunaga JI, Takai K, Matsuoka M.  HTLV-1 bZIP factor suppresses apoptosis by attenuating the function of FoxO3a and altering its localization.  Cancer Research. 2014;74(1):188–200.
  * (9) Mitobe Y, Yasunaga JI, Furuta R, Matsuoka M.  HTLV-1 bZIP Factor RNA and Protein Impart Distinct Functions on T-cell Proliferation and Survival.  Cancer Research. 2015;75(19):4143–4152.
  * (10) Ma G, Yasunaga JI, Matsuoka M.  Multifaceted functions and roles of HBZ in HTLV-1 pathogenesis.  Retrovirology. 2016;13(1).
  * (11) Kataoka K, Nagata Y, Kitanaka A, Shiraishi Y, Shimamura T, Yasunaga JI, et al.  Integrated molecular analysis of adult T cell leukemia/lymphoma.  Nat Genet. 2015;47(11):1304–1315.
  * (12) Kataoka K, Ogawa S.  Genetic landscape of adult T-cell leukemia/lymphoma.  Rinsho Ketsueki. 2016;57(4):417–424.
  * (13) Kataoka K, Shiraishi Y, Takeda Y, Sakata S, Matsumoto M, Nagano S, et al.  Aberrant PD-L1 expression through 3’-UTR disruption in multiple cancers.  Nature. 2016;534(7607):402–406.
  * (14) Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ.  Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin , DNA-binding proteins and nucleosome position.  Nat Methods. 2013;10:1213–1218.
  * (15) Tanaka A, Ishitsuka Y, Ohta H, Fujimoto A, Yasunaga JI, Matsuoka M.  Systematic clustering algorithm for chromatin accessibility data and its application to hematopoietic cells.  PLoS Comput Biol. 2020;16(11):e1008422.
  * (16) Corces MR, Granja JM, Shams S, Louie BH, Seoane JA, Zhou W, et al.  The chromatin accessibility landscape of primary human cancers.  Science. 2018;362:1–58.
  * (17) Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al.  Model-based analysis of ChIP-Seq (MACS).  Genome Biology. 2008;9(9).
  * (18) Thurman RE, Rynes E, Humbert R, Vierstra J, Maurano MT, Haugen E, et al.  The accessible chromatin landscape of the human genome.  Nature. 2012;489(7414):75–82.  doi:10.1038/nature11232.
  * (19) Richardson JH, Edwards AJ, Cruickshank JK, Rudge P, Dalgleish AG.  In vivo cellular tropism of human T-cell leukemia virus type 1.  Journal of Virology. 1990;64:5682–5687.
  * (20) Henikoff JG, Belsky JA, Krassovsky K, MacAlpine DM, Henikoff S.  Epigenome characterization at single base-pair resolution.  Proceedings of the National Academy of Sciences of the United States of America. 2011;108(45):18318–18323.
  * (21) Mazumdar C, Shen Y, Xavy S, Zhao F, Reinisch A, Li R, et al.  Leukemia-Associated Cohesin Mutants Dominantly Enforce Stem Cell Programs and Impair Human Hematopoietic Progenitor Differentiation.  Cell Stem Cell. 2015;17(6):675–688.
  * (22) Robinson MD, Oshlack A.  A scaling normalization method for differential expression analysis of RNA-seq data.  Genome Biology. 2010;11(3):R25.  doi:10.1186/gb-2010-11-3-r25.
  * (23) Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al.  Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool.  BMC Bioinformatics. 2013;14(1):128.
  * (24) Robinson MD, McCarthy DJ, Smyth GK.  edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.  Bioinformatics. 2010;26(1):139–140.
  * (25) McCarthy DJ, Chen Y, Smyth GK.  Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation.  Nucleic Acids Res. 2012;40(10):4288–4297.
  * (26) Furuta R, Yasunaga JI, Miura M, Sugata K, Saito A, Akari H, et al.  Human T-cell leukemia virus type 1 infects multiple lineage hematopoietic cells in vivo.  PLoS Pathog. 2017;13(11):e1006722.
  * (27) Dong HY, Wilkes S, Yang H.  CD71 is selectively and ubiquitously expressed at high levels in erythroid precursors of all maturation stages: a comparative immunochemical study with glycophorin A and hemoglobin A.  Am J Surg Pathol. 2011;35(5):723–732.
  * (28) Marsee DK, Pinkus GS, Yu H.  CD71 (transferrin receptor): an effective marker for erythroid precursors in bone marrow biopsy specimens.  Am J Clin Pathol. 2010;134(3):429–435.
  * (29) Feinberg MW, Wara AK, Cao Z, Lebedeva MA, Rosenbauer F, Iwasaki H, et al.  The Kruppel-like factor KLF4 is a critical regulator of monocyte differentiation.  EMBO J. 2007;26(18):4138–4148.
  * (30) Takahashi R, Yamagishi M, Nakano K, Yamochi T, Yamochi T, Fujikawa D, et al.  Epigenetic deregulation of Ellis Van Creveld confers robust Hedgehog signaling in adult T-cell leukemia.  Cancer Sci. 2014;105(9):1160–1169.
  * (31) Sasaki H, Nishikata I, Shiraga T, Akamatsu E, Fukami T, Hidaka T, et al.  Overexpression of a cell adhesion molecule, TSLC1, as a possible molecular marker for acute-type adult T-cell leukemia.  Blood. 2005;105(3):1204–1213.
  * (32) Clark TG, Conway SJ, Scott IC, Labosky PA, Winnier G, Bundy J, et al.  The mammalian Tolloid-like 1 gene, Tll1, is necessary for normal septation and positioning of the heart.  Development. 1999;126(12):2631–2642.
  * (33) Matsuura K, Sawai H, Ikeo K, Ogawa S, Iio E, Isogawa M, et al.  Genome-Wide Association Study Identifies TLL1 Variant Associated With Development of Hepatocellular Carcinoma After Eradication of Hepatitis C Virus Infection.  Gastroenterology. 2017;152(6):1383–1394.
  * (34) Ge G, Greenspan DS.  BMP1 controls TGF$\beta$1 activation via cleavage of latent TGF$\beta$-binding protein.  J Cell Biol. 2006;175(1):111–120.
  * (35) Westerhausen DR, Hopkins WE, Billadello JJ.  Multiple transforming growth factor-beta-inducible elements regulate expression of the plasminogen activator inhibitor type-1 gene in Hep G2 cells.  J Biol Chem. 1991;266(2):1092–1100.
  * (36) Cabal-Hierro L, van Galen P, Prado MA, Higby KJ, Togami K, Mowery CT, et al.  Chromatin accessibility promotes hematopoietic and leukemia stem cell activity.  Nature Communications. 2020;11(1):1–18.
  * (37) Qu K, Zaba LC, Satpathy AT, Giresi PG, Li R, Jin Y, et al.  Chromatin Accessibility Landscape of Cutaneous T Cell Lymphoma and Dynamic Response to HDAC Inhibitors.  Cancer Cell. 2017;32(1):27–41.e4.
  * (38) Li Z, Schulz MH, Look T, Begemann M, Zenke M, Costa IG.  Identification of transcription factor binding sites using ATAC-seq.  Genome Biol. 2019;20(1):45.
  * (39) Yamazaki T, Yang XO, Chung Y, Fukunaga A, Nurieva R, Pappu B, et al.  CCR6 regulates the migration of inflammatory and regulatory T cells.  J Immunol. 2008;181(12):8391–8401.
  * (40) Yao R, Ma YL, Liang W, Li HH, Ma ZJ, Yu X, et al.  MicroRNA-155 modulates Treg and Th17 cells differentiation and Th17 cell function by targeting SOCS1.  PLoS One. 2012;7(10):e46082.
  * (41) Pichler K, Schneider G, Grassmann R.  MicroRNA miR-146a and further oncogenesis-related cellular microRNAs are dysregulated in HTLV-1-transformed T lymphocytes.  Retrovirology. 2008;5:100.
  * (42) Watanabe T.  Adult T-cell leukemia: molecular basis for clonal expansion and transformation of HTLV-1-infected T cells.  Blood. 2017;129(9):1071–1081.
  * (43) Staversky RJ, Byun DK, Georger MA, Zaffuto BJ, Goodman A, Becker MW, et al.  The Chemokine CCL3 Regulates Myeloid Differentiation and Hematopoietic Stem Cell Numbers.  Sci Rep. 2018;8(1):14691.
  * (44) Sidhu SS, Yuan S, Innes AL, Kerr S, Woodruff PG, Hou L, et al.  Roles of epithelial cell-derived periostin in TGF-beta activation, collagen production, and collagen gel elasticity in asthma.  Proc Natl Acad Sci U S A. 2010;107(32):14170–14175.
  * (45) Shimoyama M.  Diagnostic criteria and classification of clinical subtypes of adult T-cell leukaemia-lymphoma. A report from the Lymphoma Study Group (1984-87).  Br J Haematol. 1991;3:428–437.
  * (46) Corces MR, Buenrostro JD, Wu B, Greenside PG, Chan SM, Koenig JL, et al.  Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution.  Nat Genetics. 2016;48:1193–1203.
  * (47) Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al.  STAR: Ultrafast universal RNA-seq aligner.  Bioinformatics. 2013;29:15–21.
  * (48) Ponchel F, Toomes C, Bransfield K, Leong FT, Douglas SH, Field SL, et al.  Real-time PCR based on SYBR-Green I fluorescence: An alternative to the TaqMan assay for a relative quantification of gene rearrangements, gene amplifications and micro gene deletions.  BMC Biotechnol. 2003;3(1):18.
  * (49) Zhao T, Satou Y, Sugata K, Miyazato P, Green PL, Imamura T, et al.  HTLV-1 bZIP factor enhances TGF-$\beta$ signaling through p300 coactivator.  Blood. 2011;118(7):1865–1876.  doi:10.1182/blood-2010-12-326199.
  * (50) Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al.  Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities.  Mol Cell. 2010;38(4):576–589.